OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

TransportModel* classes are now expected to have a custom solve() method. More:

Browse Source 
- Using new solve() method in spu_2p.
 - solve() implemented in terms of protected superclass method reorderAndTransport().
 - Removed unused code being replaced by solve().
This commit is contained in:
Atgeirr Flø Rasmussen 2012-02-10 10:48:18 +01:00
parent 11e5f76813
commit a48b261a3c
8 changed files with 116 additions and 114 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
Makefile.am
View File

@ -78,11 +78,11 @@ opm/core/pressure/mimetic/hybsys.c \
opm/core/transport/spu_explicit.c \
opm/core/transport/spu_implicit.c \
opm/core/transport/transport_source.c \
opm/core/transport/reorder/TransportModelInterface.cpp \
opm/core/transport/reorder/TransportModelTwophase.cpp \
opm/core/transport/reorder/reordersequence.c \
opm/core/transport/reorder/nlsolvers.c \
opm/core/transport/reorder/tarjan.c \
opm/core/transport/reorder/twophasetransport.cpp
opm/core/transport/reorder/tarjan.c
nobase_include_HEADERS = \
@ -188,8 +188,8 @@ opm/core/transport/reorder/TransportModelInterface.hpp \
opm/core/transport/reorder/TransportModelTwophase.hpp \
opm/core/transport/reorder/nlsolvers.h \
opm/core/transport/reorder/reordersequence.h \
opm/core/transport/reorder/tarjan.h \
opm/core/transport/reorder/twophasetransport.hpp
opm/core/transport/reorder/tarjan.h
if UMFPACK
libopmcore_la_SOURCES += opm/core/linalg/call_umfpack.c

24
examples/spu_2p.cpp
View File

@ -59,7 +59,7 @@
#include <opm/core/transport/CSRMatrixBlockAssembler.hpp>
#include <opm/core/transport/SinglePointUpwindTwoPhase.hpp>
#include <opm/core/transport/reorder/twophasetransport.hpp>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <boost/filesystem/convenience.hpp>
#include <boost/scoped_ptr.hpp>
@ -753,9 +753,14 @@ main(int argc, char** argv)
TwophaseFluid fluid(*props);
// Solvers init.
// Pressure solver.
PressureSolver psolver(grid->c_grid(), *props);
// Non-reordering solver.
TransportModel model (fluid, *grid->c_grid(), porevol, 0, guess_old_solution);
TransportSolver tsolver(model);
// Reordering solver.
Opm::TransportModelTwophase reorder_model(*grid->c_grid(), &porevol[0], *props);
// State-related and source-related variables init.
std::vector<double> totmob;
@ -763,10 +768,13 @@ main(int argc, char** argv)
// We need a separate reorder_sat, because the reorder
// code expects a scalar sw, not both sw and so.
std::vector<double> reorder_sat(grid->c_grid()->number_of_cells);
double flow_per_sec = 0.1*tot_porevol/Opm::unit::day;
// double flow_per_sec = 0.1*tot_porevol/Opm::unit::day;
double flow_per_sec = 0.1*porevol[0]/Opm::unit::day;
std::vector<double> src (grid->c_grid()->number_of_cells, 0.0);
src[0] = flow_per_sec;
src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
// src[0] = flow_per_sec;
// src[grid->c_grid()->number_of_cells - 1] = -flow_per_sec;
std::fill(src.begin(), src.begin() + src.size()/2, flow_per_sec);
std::fill(src.begin() + src.size()/2, src.end(), -flow_per_sec);
TransportSource* tsrc = create_transport_source(2, 2);
double ssrc[] = { 1.0, 0.0 };
double ssink[] = { 0.0, 1.0 };
@ -840,13 +848,7 @@ main(int argc, char** argv)
// boundary flows must be accumulated into
// source term following the same convention.
transport_timer.start();
Opm::reorderTransportTwophase(&porevol[0],
&reorder_src[0],
stepsize,
grid->c_grid(),
props.get(),
&state.faceflux()[0],
&reorder_sat[0]);
reorder_model.solve(&state.faceflux()[0], &reorder_src[0], stepsize, &reorder_sat[0]);
transport_timer.stop();
double tt = transport_timer.secsSinceStart();
std::cout << "Transport solver took: " << tt << " seconds." << std::endl;

49
opm/core/transport/reorder/TransportModelInterface.cpp Normal file
View File

@ -0,0 +1,49 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <opm/core/transport/reorder/TransportModelInterface.hpp>
#include <opm/core/transport/reorder/twophasetransport.hpp>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/grid.h>
#include <vector>
#include <cassert>
void Opm::TransportModelInterface::reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux)
{
// Compute sequence of single-cell problems
std::vector<int> sequence(grid.number_of_cells);
std::vector<int> components(grid.number_of_cells + 1);
int ncomponents;
compute_sequence(&grid, darcyflux, &sequence[0], &components[0], &ncomponents);
// Assume all strong components are single-cell domains.
assert(ncomponents == grid.number_of_cells);
for (int i = 0; i < grid.number_of_cells; ++i) {
#ifdef MATLAB_MEX_FILE
if (interrupt_signal) {
mexPrintf("Reorder loop interrupted by user: %d of %d "
"cells finished.\n", i, grid.number_of_cells);
break;
}
#endif
solveSingleCell(sequence[i]);
}
}

12
opm/core/transport/reorder/TransportModelInterface.hpp
View File

@ -20,16 +20,28 @@
#ifndef OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
#define OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED
class UnstructuredGrid;
namespace Opm
{
/// Interface for reordering transport models.
/// A transport model must provide the solveSingleCell()
/// method, and is expected to implement a solve() method
/// that will have an interface geared to the model's needs.
/// (The solve() method is therefore not virtual in this class).
/// The reorderAndTransport() method is provided as an
/// aid to implementing solve() in subclasses.
class TransportModelInterface
{
public:
virtual ~TransportModelInterface() {}
virtual void solveSingleCell(int cell) = 0;
protected:
void reorderAndTransport(const UnstructuredGrid& grid, const double* darcyflux);
};
} // namespace Opm
#endif // OPM_TRANSPORTMODELINTERFACE_HEADER_INCLUDED

47
opm/core/transport/reorder/TransportModelTwophase.cpp
View File

@ -27,28 +27,35 @@ namespace Opm
{
TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid* grid,
const Opm::IncompPropertiesInterface* props,
const double* darcyflux,
TransportModelTwophase::TransportModelTwophase(const UnstructuredGrid& grid,
const double* porevolume,
const double* source,
const double dt,
double* saturation)
const Opm::IncompPropertiesInterface& props)
: grid_(grid),
props_(props),
darcyflux_(darcyflux),
porevolume_(porevolume),
source_(source),
dt_(dt),
saturation_(saturation),
fractionalflow_(grid->number_of_cells, 0.0)
props_(props),
darcyflux_(0),
source_(0),
dt_(0.0),
saturation_(0),
fractionalflow_(grid.number_of_cells, -1.0)
{
if (props->numPhases() != 2) {
if (props.numPhases() != 2) {
THROW("Property object must have 2 phases");
}
visc_ = props->viscosity();
visc_ = props.viscosity();
}
void TransportModelTwophase::solve(const double* darcyflux,
const double* source,
const double dt,
double* saturation)
{
darcyflux_ = darcyflux;
source_ = source;
dt_ = dt;
saturation_ = saturation;
reorderAndTransport(grid_, darcyflux);
}
// Residual function r(s) for a single-cell implicit Euler transport
//
@ -73,17 +80,17 @@ namespace Opm
outflux = tm.source_[cell] <= 0 ? -tm.source_[cell] : 0.0;
dtpv = tm.dt_/tm.porevolume_[cell];
for (int i = tm.grid_->cell_facepos[cell]; i < tm.grid_->cell_facepos[cell+1]; ++i) {
int f = tm.grid_->cell_faces[i];
for (int i = tm.grid_.cell_facepos[cell]; i < tm.grid_.cell_facepos[cell+1]; ++i) {
int f = tm.grid_.cell_faces[i];
double flux;
int other;
// Compute cell flux
if (cell == tm.grid_->face_cells[2*f]) {
if (cell == tm.grid_.face_cells[2*f]) {
flux = tm.darcyflux_[f];
other = tm.grid_->face_cells[2*f+1];
other = tm.grid_.face_cells[2*f+1];
} else {
flux =-tm.darcyflux_[f];
other = tm.grid_->face_cells[2*f];
other = tm.grid_.face_cells[2*f];
}
// Add flux to influx or outflux, if interiour.
if (other != -1) {
@ -118,7 +125,7 @@ namespace Opm
{
double sat[2] = { s, 1.0 - s };
double mob[2];
props_->relperm(1, sat, &cell, mob, 0);
props_.relperm(1, sat, &cell, mob, 0);
mob[0] /= visc_[0];
mob[1] /= visc_[1];
return mob[0]/(mob[0] + mob[1]);

20
opm/core/transport/reorder/TransportModelTwophase.hpp
View File

@ -33,21 +33,23 @@ namespace Opm
class TransportModelTwophase : public TransportModelInterface
{
public:
TransportModelTwophase(const UnstructuredGrid* grid,
const Opm::IncompPropertiesInterface* props,
const double* darcyflux,
TransportModelTwophase(const UnstructuredGrid& grid,
const double* porevolume,
const double* source,
const double dt,
double* saturation);
const Opm::IncompPropertiesInterface& props);
void solve(const double* darcyflux,
const double* source,
const double dt,
double* saturation);
virtual void solveSingleCell(int cell);
private:
const UnstructuredGrid* grid_;
const IncompPropertiesInterface* props_;
const double* darcyflux_; // one flux per grid face
const UnstructuredGrid& grid_;
const double* porevolume_; // one volume per cell
const IncompPropertiesInterface& props_;
const double* darcyflux_; // one flux per grid face
const double* source_; // one source per cell
double dt_;
double* saturation_; // one per cell

46
opm/core/transport/reorder/twophasetransport.cpp
View File

@ -1,46 +0,0 @@
/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
#include <opm/core/transport/reorder/twophasetransport.hpp>
#include <opm/core/transport/reorder/reordersequence.h>
#include <opm/core/transport/reorder/TransportModelTwophase.hpp>
#include <opm/core/grid.h>
#include <cassert>
void Opm::reorderTransportTwophase(const double *porevolume,
const double *source,
const double dt,
const UnstructuredGrid *grid,
const IncompPropertiesInterface* props,
const double *darcyflux,
double *saturation)
{
// Set up transport model.
TransportModelTwophase tmodel(grid, props, darcyflux,
porevolume, source, dt, saturation);
// Compute sequence of single-cell problems
std::vector<int> sequence(grid->number_of_cells);
std::vector<int> components(grid->number_of_cells + 1);
int ncomponents;
compute_sequence(grid, darcyflux, &sequence[0], &components[0], &ncomponents);
// Assume all strong components are single-cell domains.
assert(ncomponents == grid->number_of_cells);
for (int i = 0; i < grid->number_of_cells; ++i) {
#ifdef MATLAB_MEX_FILE
if (interrupt_signal) {
mexPrintf("Reorder loop interrupted by user: %d of %d "
"cells finished.\n", i, grid->number_of_cells);
break;
}
#endif
tmodel.solveSingleCell(sequence[i]);
}
}
/* Local Variables: */
/* c-basic-offset:4 */
/* End: */

24
opm/core/transport/reorder/twophasetransport.hpp
View File

@ -1,24 +0,0 @@
/* Copyright 2011 (c) Jostein R. Natvig <Jostein.R.Natvig at sintef.no> */
#ifndef TWOPHASETRANSPORT_HPP_INCLUDED
#define TWOPHASETRANSPORT_HPP_INCLUDED
struct UnstructuredGrid;
namespace Opm
{
class IncompPropertiesInterface;
void reorderTransportTwophase(const double *porevolume,
const double *source,
const double dt,
const UnstructuredGrid *grid,
const IncompPropertiesInterface* props,
const double *darcyflux,
double *saturation);
} // namespace Opm
#endif /* TWOPHASETRANSPORT_HPP_INCLUDED */