20#ifndef OPM_COMPONENTTRANSPORT_HEADER_INCLUDED
21#define OPM_COMPONENTTRANSPORT_HEADER_INCLUDED
25#include <opm/core/utility/parameters/ParameterGroup.hpp>
42template <
class Gr
id,
class Rock,
class Flu
id,
class Wells>
49 : pgrid_(0), prock_(0), pfluid_(0), pwells_(0), ptrans_(0),
50 min_surfvol_threshold_(0.0),
51 single_step_only_(false),
56 void init(
const Opm::parameter::ParameterGroup& param)
58 min_surfvol_threshold_ = param.getDefault(
"min_surfvol_threshold", min_surfvol_threshold_);
59 single_step_only_ = param.getDefault(
"single_step_only", single_step_only_);
60 min_vtime_ = param.getDefault(
"min_vtime", min_vtime_);
67 const std::vector<double>& face_trans,
68 const typename Grid::Vector& gravity)
74 ptrans_ = &face_trans;
82 const CompVec& external_composition,
83 const std::vector<double>& face_flux,
84 const std::vector<PhaseVec>& cell_pressure,
85 const std::vector<PhaseVec>& face_pressure,
87 const double voldisclimit,
88 std::vector<CompVec>& cell_z)
90 int num_cells = pgrid_->numCells();
91 std::vector<CompVec> comp_change;
92 std::vector<double> cell_outflux;
93 std::vector<double> cell_max_ff_deriv;
94 std::vector<double> cell_gravflux;
95 std::vector<double> cell_max_mob_deriv;
96 double cur_time = 0.0;
97 updateFluidProperties(cell_pressure, face_pressure, cell_z,
98 external_pressure, external_composition);
102 std::vector<CompVec> cell_z_start;
103 std::cout <<
"Transport solver target time: " << dt << std::endl;
104 std::cout <<
" Step Stepsize Remaining time\n";
106 while (cur_time < dt) {
107 cell_z_start = cell_z;
108 computeChange(face_flux, comp_change, cell_outflux, cell_max_ff_deriv);
109 double min_time = 1e100;
110 for (
int cell = 0; cell < num_cells; ++cell) {
111 double pvol = prock_->
porosity(cell)*pgrid_->cellVolume(cell);
112 double vtime = pvol/(cell_outflux[cell]*cell_max_ff_deriv[cell]);
113 double gtime = 1e100;
114 double max_nonzero_time = 1e100;
116 if (comp_change[cell][comp] < 0.0) {
117 if (cell_z[cell][comp] > min_surfvol_threshold_) {
118 max_nonzero_time = std::min(max_nonzero_time,
119 -cell_z[cell][comp]*pvol/comp_change[cell][comp]);
121 comp_change[cell][comp] = 0.0;
122 cell_z[cell][comp] = 0.0;
126 vtime = std::max(vtime,min_vtime_);
127 double time = std::min(std::min(vtime, gtime), max_nonzero_time);
128 min_time = std::min(time, min_time);
131 double step_time = dt - cur_time;
132 if (min_time < step_time) {
133 step_time = min_time;
134 cur_time += min_time;
139 if ((dt - cur_time)/step_time > 10000) {
140 std::cout <<
"Collapsing transport stepsize detected." << std::endl;
143 std::cout.precision(10);
144 std::cout << std::setw(6) << count++
145 << std::setw(24) << step_time
146 << std::setw(24) << dt - cur_time << std::endl;
147 std::cout.precision(16);
148 for (
int cell = 0; cell < num_cells; ++cell) {
149 double pv = pgrid_->cellVolume(cell)*prock_->
porosity(cell);
150 comp_change[cell] *= (step_time/pv);
151 cell_z[cell] += comp_change[cell];
154 updateFluidProperties(cell_pressure, face_pressure, cell_z,
155 external_pressure, external_composition);
156 bool ok = volumeDiscrepancyAcceptable(voldisclimit);
159 cell_z = cell_z_start;
160 cur_time -= step_time;
163 if (single_step_only_) {
174 const Fluid* pfluid_;
175 const Wells* pwells_;
176 const std::vector<double>* ptrans_;
177 typename Grid::Vector gravity_;
181 std::vector<double> total_mobility;
182 std::vector<PhaseVec> fractional_flow;
184 template <
class G,
class R>
185 void computeNew(
const G& grid,
188 const typename Grid::Vector gravity,
189 const std::vector<PhaseVec>& cell_pressure,
190 const std::vector<PhaseVec>& face_pressure,
191 const std::vector<CompVec>& cell_z,
196 cell_pressure, face_pressure,
198 int num = grid.numCells();
199 total_mobility.resize(num);
200 fractional_flow.resize(num);
201#pragma omp parallel for
202 for (
int i = 0; i < num; ++i) {
203 total_mobility[i] = 0.0;
204 for (
int phase = 0; phase <
numPhases; ++phase) {
205 total_mobility[i] += cell_data.
mobility[i][phase];
207 fractional_flow[i] = cell_data.
mobility[i];
208 fractional_flow[i] *= (1.0/total_mobility[i]);
212 AllTransportFluidData fluid_data_;
213 struct TransportFluidData
222 TransportFluidData bdy_;
223 std::vector<int> perf_cells_;
224 std::vector<double> perf_flow_;
225 std::vector<TransportFluidData> perf_props_;
226 double min_surfvol_threshold_;
227 bool single_step_only_;
232 TransportFluidData computeProps(
const PhaseVec& pressure,
236 TransportFluidData data;
238 data.mobility = state.mobility_;
239 double total_mobility = 0.0;
240 for (
int phase = 0; phase <
numPhases; ++phase) {
241 total_mobility += state.mobility_[phase];
243 data.fractional_flow = state.mobility_;
244 data.fractional_flow /= total_mobility;
245 data.phase_to_comp = state.phase_to_comp_;
246 data.relperm = state.relperm_;
247 data.viscosity = state.viscosity_;
251 void updateFluidProperties(
const std::vector<PhaseVec>& cell_pressure,
252 const std::vector<PhaseVec>& face_pressure,
253 const std::vector<CompVec>& cell_z,
255 const CompVec& external_composition)
258 const double dummy_dt = 1.0;
259 fluid_data_.computeNew(*pgrid_, *prock_, *pfluid_, gravity_,
260 cell_pressure, face_pressure, cell_z, external_composition, dummy_dt);
263 bdy_ = computeProps(external_pressure, external_composition);
272 Wells::WellReport::report()->clearAll();
273 int num_cells = pgrid_->numCells();
274 for (
int cell = 0; cell < num_cells; ++cell) {
277 perf_cells_.push_back(cell);
278 perf_flow_.push_back(flow);
282 perf_props_.push_back(computeProps(well_pressure, well_mixture));
284 Wells::WellReport::report()->cellPressure.push_back(cell_pressure[cell][0]);
285 Wells::WellReport::report()->cellId.push_back(cell);
293 void cellData(
int cell, TransportFluidData& tfd)
const
295 tfd.saturation = fluid_data_.cell_data.saturation[cell];
296 tfd.mobility = fluid_data_.cell_data.mobility[cell];
297 tfd.fractional_flow = fluid_data_.fractional_flow[cell];
298 tfd.phase_to_comp = fluid_data_.cell_data.state_matrix[cell];
299 tfd.relperm = fluid_data_.cell_data.relperm[cell];
300 tfd.viscosity = fluid_data_.cell_data.viscosity[cell];
306 bool volumeDiscrepancyAcceptable(
const double voldisclimit)
const
308 double rel_voldiscr = *std::max_element(fluid_data_.relvoldiscr.begin(), fluid_data_.relvoldiscr.end());
309 if (rel_voldiscr > voldisclimit) {
310 std::cout <<
" Relative volume discrepancy too large: " << rel_voldiscr << std::endl;
319 void computeChange(
const std::vector<double>& face_flux,
320 std::vector<CompVec>& comp_change,
321 std::vector<double>& cell_outflux,
322 std::vector<double>& cell_max_ff_deriv)
324 const int num_cells = pgrid_->numCells();
327 comp_change.resize(num_cells,
zero);
328 cell_outflux.clear();
329 cell_outflux.resize(num_cells, 0.0);
330 cell_max_ff_deriv.clear();
331 cell_max_ff_deriv.resize(num_cells, 0.0);
332 CompVec surf_dens = pfluid_->surfaceDensities();
333 const int num_faces = pgrid_->numFaces();
334 std::vector<CompVec> face_change(num_faces);
335 std::vector<double> faces_max_ff_deriv(num_faces);
336#pragma omp parallel for
337 for (
int face = 0; face < num_faces; ++face) {
340 for (
int phase = 0; phase <
numPhases; ++phase) {
341 const double* At = &fluid_data_.face_data.state_matrix[face][0][0];
342 for (
int comp = 0; comp <
numPhases; ++comp) {
343 phase_dens[phase] += At[
numPhases*phase + comp]*surf_dens[comp];
348 TransportFluidData d[2];
349 for (
int ix = 0; ix < 2; ++ix) {
350 c[ix] = pgrid_->faceCell(face, ix);
352 cellData(c[ix], d[ix]);
362 typename Grid::Vector centroid_diff = c[0] >= 0 ? pgrid_->cellCentroid(c[0]) : pgrid_->faceCentroid(face);
363 centroid_diff -= c[1] >= 0 ? pgrid_->cellCentroid(c[1]) : pgrid_->faceCentroid(face);
364 double gravity_flux = gravity_*centroid_diff*ptrans_->operator[](face);
366 process_face(&(d[0].mobility[0]), &(d[1].mobility[0]),
367 &vstar[0], gravity_flux,
numPhases, &rho_star[0], upwind_dir);
371 double tot_mob = 0.0;
372 for (
int phase = 0; phase <
numPhases; ++phase) {
373 phase_mob[phase] = d[upwind_dir[phase] - 1].mobility[phase];
374 tot_mob += phase_mob[phase];
379 phase_flux *= face_flux[face];
381 if (gravity_flux != 0.0 && c[0] >= 0 && c[1] >= 0) {
383 double omega = ff*phase_dens;
384 for (
int phase = 0; phase <
numPhases; ++phase) {
385 double gf = (phase_dens[phase] - omega)*gravity_flux;
386 phase_flux[phase] -= phase_mob[phase]*gf;
391 double face_max_ff_deriv = 0.0;
394 int downwind_cell = c[upwind_dir[0]%2];
395 if (downwind_cell >= 0) {
398 PhaseVec upwind_viscosity = d[upwind_dir[0] - 1].viscosity;
399 PhaseVec upwind_sat = d[upwind_dir[0] - 1].saturation;
400 PhaseVec upwind_relperm = d[upwind_dir[0] - 1].relperm;
403 double downwind_totmob = 0.0;
404 double upwind_totmob = 0.0;
405 for (
int phase = 0; phase <
numPhases; ++phase) {
406 downwind_mob[phase] = fluid_data_.cell_data.relperm[downwind_cell][phase]/upwind_viscosity[phase];
407 downwind_totmob += downwind_mob[phase];
408 upwind_mob[phase] = upwind_relperm[phase]/upwind_viscosity[phase];
409 upwind_totmob += upwind_mob[phase];
411 PhaseVec downwind_ff = downwind_mob;
412 downwind_ff /= downwind_totmob;
414 upwind_ff /= upwind_totmob;
416 ff_diff -= downwind_ff;
417 for (
int phase = 0; phase <
numPhases; ++phase) {
418 if (std::fabs(ff_diff[phase]) > 1e-10) {
419 if (face_flux[face] != 0.0) {
420 double ff_deriv = ff_diff[phase]/(upwind_sat[phase] - fluid_data_.cell_data.saturation[downwind_cell][phase]);
422 face_max_ff_deriv = std::max(face_max_ff_deriv, std::fabs(ff_deriv));
427 faces_max_ff_deriv[face] = face_max_ff_deriv;
431 for (
int phase = 0; phase <
numPhases; ++phase) {
432 int upwind_ix = upwind_dir[phase] - 1;
433 CompVec z_in_phase = d[upwind_ix].phase_to_comp[phase];
434 z_in_phase *= phase_flux[phase];
435 change += z_in_phase;
437 face_change[face] = change;
441#pragma omp parallel for
442 for (
int cell = 0; cell < num_cells; ++cell) {
443 const int num_local_faces = pgrid_->numCellFaces(cell);
444 for (
int local = 0; local < num_local_faces; ++local) {
445 int face = pgrid_->cellFace(cell,local);
446 if (cell == pgrid_->faceCell(face, 0)) {
447 comp_change[cell] -= face_change[face];
448 if (face_flux[face] >= 0.0) {
449 cell_outflux[cell] += std::fabs(face_flux[face]);
451 }
else if (cell == pgrid_->faceCell(face, 1)) {
452 comp_change[cell] += face_change[face];
453 if (face_flux[face] < 0.0) {
454 cell_outflux[cell] += std::fabs(face_flux[face]);
457 cell_max_ff_deriv[cell] = std::max(cell_max_ff_deriv[cell],
458 faces_max_ff_deriv[face]);
463 int num_perf = perf_cells_.size();
464 for (
int perf = 0; perf < num_perf; ++perf) {
465 int cell = perf_cells_[perf];
466 double flow = perf_flow_[perf];
468 const TransportFluidData& fl = perf_props_[perf];
472 PhaseVec phase_flux = flow > 0.0 ? fl.fractional_flow : fluid_data_.fractional_flow[cell];
478 for (
int phase = 0; phase <
numPhases; ++phase) {
479 CompVec z_in_phase = fl.phase_to_comp[phase];
480 z_in_phase *= phase_flux[phase];
481 change += z_in_phase;
483 comp_change[cell] += change;
484 Wells::WellReport::report()->massRate.push_back(change);
497 process_face(
double *cmob1,
double *cmob2,
double *vstar,
double gf,
498 int np,
double *rho,
int *ix)
502 for (i=0; i<np; ++i) {
509 for(k=0; k<np; ++k) {
if (rho[k] > r) { r = rho[k]; j = k; } }
510 a = !(vstar[j]<0) && !(gf<0);
511 b = !(vstar[j]>0) && !(gf>0);
517 if (v < 0) { ix[j] = 2; m = cmob2[j]; }
518 else { ix[j] = 1; m = cmob1[j]; }
519 for (k=0; k<np; ++k) { vstar[k] -= m *(rho[k]-r)*gf;}
529 for(k=0; k<np; ++k) {
if (rho[k] < r) { r = rho[k]; j = k; } }
534 if (v < 0) { ix[j] = 2; m = cmob2[j]; }
535 else { ix[j] = 1; m = cmob1[j]; }
536 for (k=0; k<np; ++k) { vstar[k] -= m *(rho[k]-r)*gf;}
546 phase_upwind_directions(
int number_of_faces,
int *face_cells,
double *dflux,
double *gflux,
int np,
547 double *cmobility,
double *fdensity,
int *ix)
551 double *cmob1, *cmob2;
552 double *fden = (
double*)malloc(np *
sizeof *fden);
553 double *vstar = (
double*)malloc(np *
sizeof *vstar);
556 for (f=0; f<number_of_faces; ++f) {
559 for (k=0; k<np; ++k) {
561 fden [k] = fdensity[np*f+k];
563 c1 = face_cells[2*f + 0];
564 c2 = face_cells[2*f + 1];
565 if (c1 != -1) {cmob1 = cmobility+np*c1;}
566 else {cmob1 = cmobility+np*c2;}
567 if (c2 != -1) {cmob2 = cmobility+np*c2;}
568 else {cmob2 = cmobility+np*c1;}
569 process_face(cmob1, cmob2, vstar, gf, np, fden, ix+np*f);
Definition: BlackoilDefs.hpp:31
Dune::FieldVector< Scalar, numComponents > CompVec
Definition: BlackoilDefs.hpp:40
Dune::FieldVector< Scalar, numPhases > PhaseVec
Definition: BlackoilDefs.hpp:41
Dune::FieldMatrix< Scalar, numComponents, numPhases > PhaseToCompMatrix
Definition: BlackoilDefs.hpp:43
@ numComponents
Definition: BlackoilDefs.hpp:33
@ numPhases
Definition: BlackoilDefs.hpp:34
Definition: BlackoilFluid.hpp:42
FluidStateBlackoil FluidState
Definition: BlackoilFluid.hpp:44
Definition: ComponentTransport.hpp:44
ExplicitCompositionalTransport()
Default constructor. Does nothing.
Definition: ComponentTransport.hpp:48
void init(const Opm::parameter::ParameterGroup ¶m)
Definition: ComponentTransport.hpp:56
void setup(const Grid &grid, const Rock &rock, const Fluid &fluid, const Wells &wells, const std::vector< double > &face_trans, const typename Grid::Vector &gravity)
Definition: ComponentTransport.hpp:63
double transport(const PhaseVec &external_pressure, const CompVec &external_composition, const std::vector< double > &face_flux, const std::vector< PhaseVec > &cell_pressure, const std::vector< PhaseVec > &face_pressure, const double dt, const double voldisclimit, std::vector< CompVec > &cell_z)
Definition: ComponentTransport.hpp:81
A property class for porous media rock.
Definition: ImplicitTransportDefs.hpp:80
double porosity(int cell_index) const
Read-access to porosity.
Definition: Rock_impl.hpp:76
Dune::FieldVector< double, 3 > injectionMixture(int cell) const
Definition: Wells.hpp:141
double perforationPressure(int cell) const
Definition: Wells.hpp:130
double wellToReservoirFlux(int cell) const
Definition: Wells.hpp:136
Definition: BlackoilFluid.hpp:32
void zero(Matrix &A)
Zero-fill a.
Definition: Matrix.hpp:603
Definition: BlackoilFluid.hpp:406
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< PhaseVec > mobility
Definition: FluidStateBlackoil.hpp:84
PhaseVec saturation_
Definition: FluidStateBlackoil.hpp:42