OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-11-29 20:43:49 -06:00
Code Issues Packages Projects Releases Wiki Activity

Implement solveMultiCell() properly.

Browse Source 
Interface change: solver now requires a linear solver (for the multi-cell blocks only).

Implementation uses new private method assembleSingleCell(), that is a modified copy
of solveSingleCell(). Should refactor.
This commit is contained in:
Atgeirr Flø Rasmussen 2013-04-17 12:58:15 +02:00
parent f8db9cad66
commit 25a2c3d00b
3 changed files with 120 additions and 3 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
examples/compute_tof_from_files.cpp
View File

@ -120,6 +120,9 @@ main(int argc, char** argv)
}
}
// Linear solver.
LinearSolverFactory linsolver(param);
// Choice of tof solver.
bool use_dg = param.getDefault("use_dg", false);
bool use_multidim_upwind = false;
@ -163,7 +166,7 @@ main(int argc, char** argv)
if (use_dg) {
dg_solver->solveTof(&flux[0], &porevol[0], &src[0], tof);
} else {
Opm::TofReorder tofsolver(grid, use_multidim_upwind);
Opm::TofReorder tofsolver(grid, linsolver, use_multidim_upwind);
if (compute_tracer) {
tofsolver.solveTofTracer(&flux[0], &porevol[0], &src[0], tof, tracer);
} else {

106
opm/core/tof/TofReorder.cpp
View File

@ -20,6 +20,7 @@
#include "config.h"
#include <opm/core/tof/TofReorder.hpp>
#include <opm/core/grid.h>
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <algorithm>
#include <numeric>
@ -30,17 +31,21 @@ namespace Opm
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] gri d A 2d or 3d grid.
/// \param[in] linsolver Linear solver used for multi-cell blocks.
/// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding.
TofReorder::TofReorder(const UnstructuredGrid& grid,
const LinearSolverInterface& linsolver,
const bool use_multidim_upwind)
: grid_(grid),
linsolver_(linsolver),
darcyflux_(0),
porevolume_(0),
source_(0),
tof_(0),
tracer_(0),
num_tracers_(0),
block_index_(grid.number_of_cells, OutsideBlock),
use_multidim_upwind_(use_multidim_upwind)
{
}
@ -195,6 +200,55 @@ namespace Opm
void TofReorder::assembleSingleCell(const int cell,
std::vector<int>& local_column,
std::vector<double>& local_coefficient,
double& rhs)
{
// Compute flux terms.
// Sources have zero tof, and therefore do not contribute
// to upwind_term. Sinks on the other hand, must be added
// to the downwind_flux (note sign change resulting from
// different sign conventions: pos. source is injection,
// pos. flux is outflow).
double downwind_flux = std::max(-source_[cell], 0.0);
double upwind_term = 0.0;
for (int i = grid_.cell_facepos[cell]; i < grid_.cell_facepos[cell+1]; ++i) {
int f = grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == grid_.face_cells[2*f]) {
flux = darcyflux_[f];
other = grid_.face_cells[2*f+1];
} else {
flux =-darcyflux_[f];
other = grid_.face_cells[2*f];
}
// Add flux to upwind_term or downwind_flux
if (flux < 0.0) {
// Using tof == 0 on inflow, so we only add a
// nonzero contribution if we are on an internal
// face.
if (other != -1) {
if (block_index_[other] == OutsideBlock) {
upwind_term += flux*tof_[other];
} else {
local_column.push_back(block_index_[other]);
local_coefficient.push_back(flux);
}
}
} else {
downwind_flux += flux;
}
}
local_column.push_back(block_index_[cell]);
local_coefficient.push_back(downwind_flux);
rhs = porevolume_[cell] - upwind_term;
}
void TofReorder::solveSingleCellMultidimUpwind(const int cell)
{
@ -247,10 +301,60 @@ namespace Opm
void TofReorder::solveMultiCell(const int num_cells, const int* cells)
{
#if 0
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
for (int i = 0; i < num_cells; ++i) {
solveSingleCell(cells[i]);
}
#else
std::cout << "Pretending to solve multi-cell dependent equation with " << num_cells << " cells." << std::endl;
// Give each cell a local index in this multi-cell block.
for (int ci = 0; ci < num_cells; ++ci) {
block_index_[cells[ci]] = ci;
}
// Create sparse matrix and rhs vector for block equation system.
std::vector<int> ia(num_cells+1);
std::vector<int> ja;
std::vector<double> sa;
ja.reserve(2*num_cells);
sa.reserve(2*num_cells);
std::vector<double> rhs(num_cells);
std::vector< std::pair<int, double> > rowdata;
for (int ci = 0; ci < num_cells; ++ci) {
ia[ci] = ja.size();
assembleSingleCell(cells[ci], ja, sa, rhs[ci]);
// We standardize sparse row format: must sort row content by column index.
int num_row_elem = ja.size() - ia[ci];
rowdata.resize(num_row_elem);
for (int i = 0; i < num_row_elem; ++i) {
rowdata[i].first = ja[ia[ci] + i];
rowdata[i].second = sa[ia[ci] + i];
}
std::sort(rowdata.begin(), rowdata.end());
for (int i = 0; i < num_row_elem; ++i) {
ja[ia[ci] + i] = rowdata[i].first;
sa[ia[ci] + i] = rowdata[i].second;
}
}
ia.back() = ja.size();
ASSERT(ja.size() == sa.size());
// Solve system.
std::vector<double> tof_block(num_cells);
LinearSolverInterface::LinearSolverReport rep
= linsolver_.solve(num_cells, ja.size(), &ia[0], &ja[0], &sa[0], &rhs[0], &tof_block[0]);
if (!rep.converged) {
THROW("Multicell system with " << num_cells << " failed to converge.");
}
// Write to global tof vector, and reset block indices to make
// it usable for next solveMultiCell() call.
for (int ci = 0; ci < num_cells; ++ci) {
tof_[cells[ci]] = tof_block[ci];
block_index_[cells[ci]] = OutsideBlock;
}
#endif
}

12
opm/core/tof/TofReorder.hpp
View File

@ -30,6 +30,7 @@ namespace Opm
{
class IncompPropertiesInterface;
class LinearSolverInterface;
/// Implements a first-order finite volume solver for
/// (single-phase) time-of-flight using reordering.
@ -42,9 +43,11 @@ namespace Opm
{
public:
/// Construct solver.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] grid A 2d or 3d grid.
/// \param[in] linsolver Linear solver used for multi-cell blocks.
/// \param[in] use_multidim_upwind If true, use multidimensional tof upwinding.
TofReorder(const UnstructuredGrid& grid,
const LinearSolverInterface& linsolver,
const bool use_multidim_upwind = false);
/// Solve for time-of-flight.
@ -79,6 +82,10 @@ namespace Opm
private:
virtual void solveSingleCell(const int cell);
void solveSingleCellMultidimUpwind(const int cell);
void assembleSingleCell(const int cell,
std::vector<int>& local_column,
std::vector<double>& local_coefficient,
double& rhs);
virtual void solveMultiCell(const int num_cells, const int* cells);
void multidimUpwindTerms(const int face, const int upwind_cell,
@ -86,12 +93,15 @@ namespace Opm
private:
const UnstructuredGrid& grid_;
const LinearSolverInterface& linsolver_;
const double* darcyflux_; // one flux per grid face
const double* porevolume_; // one volume per cell
const double* source_; // one volumetric source term per cell
double* tof_;
double* tracer_;
int num_tracers_;
enum { OutsideBlock = -1 };
std::vector<int> block_index_; // For solveMultiCell().
bool use_multidim_upwind_;
std::vector<double> face_tof_; // For multidim upwind face tofs.
mutable std::vector<int> adj_faces_; // For multidim upwind logic.