Fix initialization of tracer heads, and output ordering.

This commit is contained in:
Atgeirr Flø Rasmussen 2014-06-25 09:56:37 +02:00
parent 3290c41f3a
commit fd606d907a

View File

@ -110,14 +110,15 @@ namespace Opm
darcyflux_ = darcyflux; darcyflux_ = darcyflux;
porevolume_ = porevolume; porevolume_ = porevolume;
source_ = source; source_ = source;
const int num_cells = grid_.number_of_cells;
#ifndef NDEBUG #ifndef NDEBUG
// Sanity check for sources. // Sanity check for sources.
const double cum_src = std::accumulate(source, source + grid_.number_of_cells, 0.0); const double cum_src = std::accumulate(source, source + num_cells, 0.0);
if (std::fabs(cum_src) > *std::max_element(source, source + grid_.number_of_cells)*1e-2) { if (std::fabs(cum_src) > *std::max_element(source, source + num_cells)*1e-2) {
OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src); OPM_THROW(std::runtime_error, "Sources do not sum to zero: " << cum_src);
} }
#endif #endif
tof.resize(grid_.number_of_cells); tof.resize(num_cells);
std::fill(tof.begin(), tof.end(), 0.0); std::fill(tof.begin(), tof.end(), 0.0);
tof_ = &tof[0]; tof_ = &tof[0];
@ -134,28 +135,36 @@ namespace Opm
// Find the tracer heads (injectors). // Find the tracer heads (injectors).
const int num_tracers = tracerheads.size(); const int num_tracers = tracerheads.size();
tracer.resize(grid_.number_of_cells*num_tracers); tracer.resize(num_cells*num_tracers);
std::fill(tracer.begin(), tracer.end(), 0.0); std::fill(tracer.begin(), tracer.end(), 0.0);
if (num_tracers > 0) { if (num_tracers > 0) {
tracerhead_by_cell_.clear(); tracerhead_by_cell_.clear();
tracerhead_by_cell_.resize(grid_.number_of_cells, NoTracerHead); tracerhead_by_cell_.resize(num_cells, NoTracerHead);
} }
for (int tr = 0; tr < num_tracers; ++tr) { for (int tr = 0; tr < num_tracers; ++tr) {
for (int i = 0; i < tracerheads[tr].size(); ++i) { for (int i = 0; i < tracerheads[tr].size(); ++i) {
const int cell = tracerheads[tr][i]; const int cell = tracerheads[tr][i];
tracer[cell*num_tracers + tr] = 1.0; tracer[num_cells * tr + cell] = 1.0;
tracerhead_by_cell_[cell] = tr; tracerhead_by_cell_[cell] = tr;
} }
} }
// Execute solve for tracers. // Execute solve for tracers.
std::vector<double> fake_pv(grid_.number_of_cells, 0.0); std::vector<double> fake_pv(num_cells, 0.0);
porevolume_ = fake_pv.data(); porevolume_ = fake_pv.data();
for (int tr = 0; tr < num_tracers; ++tr) { for (int tr = 0; tr < num_tracers; ++tr) {
tof_ = tracer.data() + tr * grid_.number_of_cells; tof_ = tracer.data() + tr * num_cells;
compute_tracer_ = true; compute_tracer_ = true;
executeSolve(); executeSolve();
} }
// Write output tracer data (transposing the computed data).
std::vector<double> computed = tracer;
for (int cell = 0; cell < num_cells; ++cell) {
for (int tr = 0; tr < num_tracers; ++tr) {
tracer[num_tracers * cell + tr] = computed[num_cells * tr + cell];
}
}
} }