Using IncompTpfa.

This commit is contained in:
Atgeirr Flø Rasmussen 2012-02-20 13:59:34 +01:00
parent b84c957e2b
commit 97134499c5
2 changed files with 4 additions and 62 deletions

View File

@ -59,63 +59,6 @@ namespace Opm
{
class PressureSolver
{
public:
PressureSolver(const UnstructuredGrid* g,
const IncompPropertiesInterface& props)
: htrans_(g->cell_facepos[ g->number_of_cells ]),
trans_ (g->number_of_faces),
gpress_(g->cell_facepos[ g->number_of_cells ])
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_htrans_compute(gg, props.permeability(), &htrans_[0]);
h_ = ifs_tpfa_construct(gg);
}
~PressureSolver()
{
ifs_tpfa_destroy(h_);
}
template <class State>
void
solve(const UnstructuredGrid* g ,
const ::std::vector<double>& totmob,
const ::std::vector<double>& src ,
State& state )
{
UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(g);
tpfa_eff_trans_compute(gg, &totmob[0], &htrans_[0], &trans_[0]);
// No gravity
std::fill(gpress_.begin(), gpress_.end(), double(0.0));
ifs_tpfa_assemble(gg, &trans_[0], &src[0], &gpress_[0], h_);
using ImplicitTransportLinAlgSupport::CSRMatrixUmfpackSolver;
CSRMatrixUmfpackSolver linsolve;
linsolve.solve(h_->A, h_->b, h_->x);
ifs_tpfa_press_flux(gg, &trans_[0], h_,
&state.pressure()[0],
&state.faceflux()[0]);
}
private:
::std::vector<double> htrans_;
::std::vector<double> trans_ ;
::std::vector<double> gpress_;
struct ifs_tpfa_data* h_;
};
void
compute_porevolume(const UnstructuredGrid* g,
const Opm::IncompPropertiesInterface& props,

View File

@ -22,8 +22,7 @@
#include "Utilities.hpp"
#include <opm/core/pressure/tpfa/ifs_tpfa.h>
#include <opm/core/pressure/tpfa/trans_tpfa.h>
#include <opm/core/pressure/IncompTpfa.hpp>
#include <opm/core/grid.h>
#include <opm/core/GridManager.hpp>
@ -34,7 +33,6 @@
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/fluid/SimpleFluid2p.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
@ -307,12 +305,13 @@ main(int argc, char** argv)
double tot_porevol = std::accumulate(porevol.begin(), porevol.end(), 0.0);
// Solvers init.
Opm::PressureSolver psolver(grid->c_grid(), *props);
Opm::IncompTpfa psolver(*grid->c_grid(), props->permeability(), 0);
Opm::TransportModelPolymer tmodel(*grid->c_grid(), props->porosity(), &porevol[0], *props, polydata, method);
// State-related and source-related variables init.
std::vector<double> totmob;
std::vector<double> omega; // Empty dummy unless/until we include gravity here.
double init_sat = param.getDefault("init_sat", 0.0);
ReservoirState state(grid->c_grid(), props->numPhases(), init_sat);
// We need a separate reorder_sat, because the reorder
@ -365,7 +364,7 @@ main(int argc, char** argv)
compute_totmob(*props, state.saturation(), totmob);
pressure_timer.start();
psolver.solve(grid->c_grid(), totmob, src, state);
psolver.solve(totmob, omega, src, state.pressure(), state.faceflux());
pressure_timer.stop();
double pt = pressure_timer.secsSinceStart();
std::cout << "Pressure solver took: " << pt << " seconds." << std::endl;