OilfieldToolsNet/opm-upscaling
4
0
Fork 0
mirror of https://github.com/OPM/opm-upscaling.git synced 2025-02-25 18:45:23 -06:00
Code Issues Packages Projects Releases Wiki Activity

Refactored linear solver interface.

Browse Source 
Now the ISTL solver is providing both AMG and non-AMG solvers.
This commit is contained in:
Atgeirr Flø Rasmussen
2010-11-02 21:10:21 +01:00
parent f189171155
commit 2e115bc386
6 changed files with 173 additions and 121 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
.hgsubstate
View File

@@ -1 +1 @@
0dfeeee61f230b304aee7b4069d3f2b96e1341cd dune/porsol/opmpressure
f26cf89a0b9e87a54d4a8166d27f71871f589e10 dune/porsol/opmpressure

21
dune/porsol/common/AbstractLinearSolver.hpp
View File

@@ -27,13 +27,6 @@
namespace Dune
{
struct LinearSolverResults
{
bool converged;
int iterations;
double reduction;
};
class AbstractLinearSolver
{
@@ -41,12 +34,20 @@ namespace Dune
virtual ~AbstractLinearSolver()
{
}
virtual void init(const parameter::ParameterGroup& param) = 0;
struct LinearSolverResults
{
bool converged;
int iterations;
double reduction;
};
virtual LinearSolverResults solve(int size, int nonzeros,
const int* ia, const int* ja, const double* sa,
const double* rhs, double* solution,
double relative_residual_tolerance,
int verbosity_level) = 0;
const double* rhs, double* solution) = 0;
};

145
dune/porsol/common/LinearSolverISTL.cpp
View File

@@ -29,42 +29,75 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/io.hh>
#include <dune/istl/overlappingschwarz.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <stdexcept>
namespace Dune
{
namespace {
typedef FieldVector<double, 1 > VectorBlockType;
typedef FieldMatrix<double, 1, 1> MatrixBlockType;
typedef BCRSMatrix <MatrixBlockType> Mat;
typedef BlockVector<VectorBlockType> Vector;
typedef MatrixAdapter<Mat,Vector,Vector> Operator;
LinearSolverISTL::LinearSolverResults
solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int verbosity);
LinearSolverISTL::LinearSolverResults
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int verbosity);
} // anonymous namespace
LinearSolverISTL::LinearSolverISTL()
: linsolver_residual_tolerance_(1e-8),
linsolver_verbosity_(0),
linsolver_type_(CG_AMG),
linsolver_save_system_(false)
{
}
LinearSolverISTL::~LinearSolverISTL()
{
}
void LinearSolverISTL::init(const parameter::ParameterGroup& param)
{
linsolver_residual_tolerance_ = param.getDefault("linsolver_residual_tolerance", linsolver_residual_tolerance_);
linsolver_verbosity_ = param.getDefault("linsolver_verbosity", linsolver_verbosity_);
linsolver_type_ = LinsolverType(param.getDefault("linsolver_type", int(linsolver_type_)));
linsolver_save_system_ = param.getDefault("linsolver_save_system", linsolver_save_system_);
if (linsolver_save_system_) {
linsolver_save_filename_ = param.getDefault("linsolver_save_filename", std::string("linsys"));
}
}
LinearSolverResults LinearSolverISTL::solve(int size, int nonzeros,
const int* ia, const int* ja, const double* sa,
const double* rhs, double* solution,
double relative_residual_tolerance,
int verbosity_level)
LinearSolverISTL::LinearSolverResults
LinearSolverISTL::solve(int size, int nonzeros,
const int* ia, const int* ja, const double* sa,
const double* rhs, double* solution)
{
// Build ISTL structures from input.
typedef FieldVector<double, 1 > VectorBlockType;
typedef FieldMatrix<double, 1, 1> MatrixBlockType;
typedef BCRSMatrix <MatrixBlockType> Matrix;
typedef BlockVector<VectorBlockType> Vector;
typedef MatrixAdapter<Matrix,Vector,Vector> Operator;
// System matrix
Matrix A(size, size, nonzeros, Matrix::row_wise);
for (Matrix::CreateIterator row = A.createbegin(); row != A.createend(); ++row) {
Mat A(size, size, nonzeros, Mat::row_wise);
for (Mat::CreateIterator row = A.createbegin(); row != A.createend(); ++row) {
int ri = row.index();
for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
row.insert(ja[i]);
@@ -82,10 +115,11 @@ namespace Dune
Vector x(size);
x = 0.0;
if (linsolver_save_system_)
{
// Save system
writeMatrixToMatlab(A, "ifshmat");
std::string rhsfile("ifshrhs");
// Save system to files.
writeMatrixToMatlab(A, linsolver_save_filename_ + "-mat");
std::string rhsfile(linsolver_save_filename_ + "-rhs");
std::ofstream rhsf(rhsfile.c_str());
rhsf.precision(15);
rhsf.setf(std::ios::scientific | std::ios::showpos);
@@ -93,7 +127,53 @@ namespace Dune
std::ostream_iterator<VectorBlockType>(rhsf, "\n"));
}
LinearSolverResults res;
switch (linsolver_type_) {
case CG_ILU0:
res = solveCG_ILU0(A, x, b, linsolver_residual_tolerance_, linsolver_verbosity_);
case CG_AMG:
res = solveCG_AMG(A, x, b, linsolver_residual_tolerance_, linsolver_verbosity_);
default:
std::cerr << "Unknown linsolver_type: " << int(linsolver_type_) << '\n';
throw std::runtime_error("Unknown linsolver_type");
}
std::copy(x.begin(), x.end(), solution);
return res;
}
namespace
{
LinearSolverISTL::LinearSolverResults
solveCG_ILU0(const Mat& A, Vector& x, Vector& b, double tolerance, int verbosity)
{
Operator opA(A);
// Construct preconditioner.
SeqILU0<Mat,Vector,Vector> precond(A, 1.0);
// Construct linear solver.
CGSolver<Vector> linsolve(opA, precond, tolerance, A.N(), verbosity);
// Solve system.
InverseOperatorResult result;
linsolve.apply(x, b, result);
// Output results.
LinearSolverISTL::LinearSolverResults res;
res.converged = result.converged;
res.iterations = result.iterations;
res.reduction = result.reduction;
return res;
}
LinearSolverISTL::LinearSolverResults
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int verbosity)
{
// Solve with AMG solver.
#define FIRST_DIAGONAL 1
#define SYMMETRIC 1
@@ -107,42 +187,49 @@ namespace Dune
#endif
#if SYMMETRIC
typedef Amg::SymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
typedef Amg::SymmetricCriterion<Mat,CouplingMetric> CriterionBase;
#else
typedef Amg::UnSymmetricCriterion<Matrix,CouplingMetric> CriterionBase;
typedef Amg::UnSymmetricCriterion<Mat,CouplingMetric> CriterionBase;
#endif
#if SMOOTHER_ILU
typedef SeqILU0<Matrix,Vector,Vector> Smoother;
typedef SeqILU0<Mat,Vector,Vector> Smoother;
#else
typedef SeqSSOR<Matrix,Vector,Vector> Smoother;
typedef SeqSSOR<Mat,Vector,Vector> Smoother;
#endif
typedef Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Amg::AMG<Operator,Vector,Smoother> Precond;
Operator opA(A);
// Construct preconditioner.
double relax = 1;
Precond::SmootherArgs smootherArgs;
smootherArgs.relaxationFactor = relax;
Criterion criterion;
criterion.setDebugLevel(verbosity_level);
criterion.setDebugLevel(verbosity);
#if ANISOTROPIC_3D
criterion.setDefaultValuesAnisotropic(3, 2);
#endif
Precond precond(opA, criterion, smootherArgs);
CGSolver<Vector> linsolve(opA, precond, relative_residual_tolerance, size, verbosity_level);
// Construct linear solver.
CGSolver<Vector> linsolve(opA, precond, tolerance, A.N(), verbosity);
// Solve system.
InverseOperatorResult result;
linsolve.apply(x, b, result);
std::copy(x.begin(), x.end(), solution);
LinearSolverResults res;
// Output results.
LinearSolverISTL::LinearSolverResults res;
res.converged = result.converged;
res.iterations = result.iterations;
res.reduction = result.reduction;
return res;
}
} // anonymous namespace
} // namespace Dune

27
dune/porsol/common/LinearSolverISTL.hpp
View File

@@ -38,21 +38,6 @@
#include <dune/porsol/common/AbstractLinearSolver.hpp>
// TODO: clean up includes.
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/io.hh>
#include <dune/istl/overlappingschwarz.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/pinfo.hh>
namespace Dune
{
@@ -60,13 +45,19 @@ namespace Dune
class LinearSolverISTL : public AbstractLinearSolver
{
public:
LinearSolverISTL();
virtual ~LinearSolverISTL();
virtual void init(const parameter::ParameterGroup& param);
virtual LinearSolverResults solve(int size, int nonzeros,
const int* ia, const int* ja, const double* sa,
const double* rhs, double* solution,
double relative_residual_tolerance,
int verbosity_level);
const double* rhs, double* solution);
private:
double linsolver_residual_tolerance_;
int linsolver_verbosity_;
enum LinsolverType { CG_ILU0 = 0, CG_AMG = 1 };
LinsolverType linsolver_type_;
bool linsolver_save_system_;
std::string linsolver_save_filename_;
};

44
dune/porsol/mimetic/IfshInterface.hpp
View File

@@ -176,7 +176,7 @@ namespace Dune
/// Type 0 selects a BiCGStab solver, type 1 selects AMG/CG.
///
template<class FluidInterface>
void solve(const FluidInterface& fl ,
void solve(const FluidInterface& fl ,
const std::vector<double>& sat,
const BCInterface& bc ,
const std::vector<double>& src,
@@ -185,13 +185,11 @@ namespace Dune
int linsolver_type = 1,
bool same_matrix = false)
{
if (linsolver_type != 1) {
MESSAGE("Requested linear solver type " << linsolver_type << ", but we only have AMG so far.");
}
if (same_matrix) {
MESSAGE("Requested reuse of preconditioner, not implemented so far.");
}
// Build totmob and omega.
int num_cells = sat.size();
std::vector<double> totmob(num_cells, 1.0);
std::vector<double> omega(num_cells, 0.0);
@@ -205,6 +203,8 @@ namespace Dune
double f_w = mob[0]/(mob[0] + mob[1]);
omega[cell] = rho[0]*f_w + rho[1]*(1.0 - f_w);
}
// Build bctypes and bcvalues.
int num_faces = pgrid_->numberOfFaces();
std::vector<HybridPressureSolver::FlowBCTypes> bctypes(num_faces, HybridPressureSolver::FBC_UNSET);
std::vector<double> bcvalues(num_faces, 0.0);
@@ -217,42 +217,26 @@ namespace Dune
} else if (face_bc.isNeumann()) {
bctypes[face] = HybridPressureSolver::FBC_FLUX;
bcvalues[face] = face_bc.outflux(); // TODO: may have to switch sign here depending on orientation.
if (bcvalues[face] != 0.0) {
THROW("Nonzero Neumann conditions not yet properly implemented (signs must be fixed)");
}
} else {
THROW("Unhandled boundary condition type.");
}
}
// typedef typename GridInterface::CellIterator CI;
// typedef typename CI::FaceIterator FI;
// for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
// for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
// if (f->boundary()) {
// int face = f->index();
// int bid = f->boundaryId();
// FlowBC face_bc = bc.flowCond(bid);
// if (face_bc.isDirichlet()) {
// bctypes[face] = HybridPressureSolver::FBC_PRESSURE;
// bcvalues[face] = face_bc.pressure();
// } else if (face_bc.isNeumann()) {
// bctypes[face] = HybridPressureSolver::FBC_FLUX;
// bcvalues[face] = face_bc.outflux(); // TODO: may have to switch sign here depending on orientation.
// } else {
// THROW("Unhandled boundary condition type.");
// }
// }
// }
// }
// Assemble system matrix and rhs.
ifsh_.assemble(src, totmob, omega, bctypes, bcvalues);
// static int count = 0;
// ++count;
// printSystem(std::string("linsys_mimetic-") + boost::lexical_cast<std::string>(count));
// Solve system.
HybridPressureSolver::LinearSystem s;
ifsh_.linearSystem(s);
LinearSolverResults res = linsolver_.solve(s.n, s.nnz, s.ia, s.ja, s.sa, s.b, s.x,
residual_tolerance, linsolver_verbosity);
parameter::ParameterGroup params;
params.insertParameter("linsolver_tolerance", boost::lexical_cast<std::string>(residual_tolerance));
params.insertParameter("linsolver_verbosity", boost::lexical_cast<std::string>(linsolver_verbosity));
params.insertParameter("linsolver_type", boost::lexical_cast<std::string>(linsolver_type));
linsolver_.init(params);
LinearSolverISTL::LinearSolverResults res = linsolver_.solve(s.n, s.nnz, s.ia, s.ja, s.sa, s.b, s.x);
if (!res.converged) {
THROW("Linear solver failed to converge in " << res.iterations << " iterations.\n"
<< "Residual reduction achieved is " << res.reduction << '\n');

55
dune/porsol/mimetic/TpfaInterface.hpp
View File

@@ -160,7 +160,7 @@ namespace Dune
/// Type 0 selects a BiCGStab solver, type 1 selects AMG/CG.
///
template<class FluidInterface>
void solve(const FluidInterface& r ,
void solve(const FluidInterface& fl ,
const std::vector<double>& sat,
const BCInterface& bc ,
const std::vector<double>& src,
@@ -169,21 +169,26 @@ namespace Dune
int linsolver_type = 1,
bool same_matrix = false)
{
if (linsolver_type != 1) {
MESSAGE("Requested linear solver type " << linsolver_type << ", but we only have AMG so far.");
}
if (same_matrix) {
MESSAGE("Requested reuse of preconditioner, not implemented so far.");
}
// Build totmob and omega.
int num_cells = sat.size();
std::vector<double> totmob(num_cells, 1.0);
std::vector<double> omega(num_cells, 0.0);
boost::array<double, FluidInterface::NumberOfPhases> mob ;
boost::array<double, FluidInterface::NumberOfPhases> rho ;
BOOST_STATIC_ASSERT(FluidInterface::NumberOfPhases == 2);
for (int cell = 0; cell < num_cells; ++cell) {
totmob[cell] = r.totalMobility(cell, sat[cell]);
double f_w = r.fractionalFlow(cell, sat[cell]);
omega[cell] = r.densityFirstPhase()*f_w + r.densitySecondPhase()*(1.0 - f_w);
fl.phaseMobilities(cell, sat[cell], mob);
fl.phaseDensities (cell, rho);
totmob[cell] = mob[0] + mob[1];
double f_w = mob[0]/(mob[0] + mob[1]);
omega[cell] = rho[0]*f_w + rho[1]*(1.0 - f_w);
}
// Build bctypes and bcvalues.
int num_faces = pgrid_->numberOfFaces();
std::vector<TPFAPressureSolver::FlowBCTypes> bctypes(num_faces, TPFAPressureSolver::FBC_UNSET);
std::vector<double> bcvalues(num_faces, 0.0);
@@ -196,42 +201,26 @@ namespace Dune
} else if (face_bc.isNeumann()) {
bctypes[face] = TPFAPressureSolver::FBC_FLUX;
bcvalues[face] = face_bc.outflux(); // TODO: may have to switch sign here depending on orientation.
if (bcvalues[face] != 0.0) {
THROW("Nonzero Neumann conditions not yet properly implemented (signs must be fixed)");
}
} else {
THROW("Unhandled boundary condition type.");
}
}
// typedef typename GridInterface::CellIterator CI;
// typedef typename CI::FaceIterator FI;
// for (CI c = pgrid_->cellbegin(); c != pgrid_->cellend(); ++c) {
// for (FI f = c->facebegin(); f != c-> faceend(); ++f) {
// if (f->boundary()) {
// int face = f->index();
// int bid = f->boundaryId();
// FlowBC face_bc = bc.flowCond(bid);
// if (face_bc.isDirichlet()) {
// bctypes[face] = TPFAPressureSolver::FBC_PRESSURE;
// bcvalues[face] = face_bc.pressure();
// } else if (face_bc.isNeumann()) {
// bctypes[face] = TPFAPressureSolver::FBC_FLUX;
// bcvalues[face] = face_bc.outflux(); // TODO: may have to switch sign here depending on orientation.
// } else {
// THROW("Unhandled boundary condition type.");
// }
// }
// }
// }
// Assemble system matrix and rhs.
psolver_.assemble(src, totmob, omega, bctypes, bcvalues);
// static int count = 0;
// ++count;
// printSystem(std::string("linsys_mimetic-") + boost::lexical_cast<std::string>(count));
// Solve system.
TPFAPressureSolver::LinearSystem s;
psolver_.linearSystem(s);
LinearSolverResults res = linsolver_.solve(s.n, s.nnz, s.ia, s.ja, s.sa, s.b, s.x,
residual_tolerance, linsolver_verbosity);
parameter::ParameterGroup params;
params.insertParameter("linsolver_tolerance", boost::lexical_cast<std::string>(residual_tolerance));
params.insertParameter("linsolver_verbosity", boost::lexical_cast<std::string>(linsolver_verbosity));
params.insertParameter("linsolver_type", boost::lexical_cast<std::string>(linsolver_type));
linsolver_.init(params);
LinearSolverISTL::LinearSolverResults res = linsolver_.solve(s.n, s.nnz, s.ia, s.ja, s.sa, s.b, s.x);
if (!res.converged) {
THROW("Linear solver failed to converge in " << res.iterations << " iterations.\n"
<< "Residual reduction achieved is " << res.reduction << '\n');