OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2024-12-21 23:13:27 -06:00
Code Issues Packages Projects Releases Wiki Activity

Moved ISTL right and left hand side construction for ISTL to solveSystem

Browse Source 
cherry-picked from add-scalarproduct and ported.
This commit is contained in:
Markus Blatt 2014-03-25 09:58:07 +01:00
parent bdc4573e2a
commit 3115b7868f
2 changed files with 28 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

47
opm/core/linalg/LinearSolverIstl.cpp
View File

@ -154,30 +154,12 @@ namespace Opm
A[ri][ja[i]] = sa[i]; A[ri][ja[i]] = sa[i];
} }
} }
// System RHS
Vector b(size);
std::copy(rhs, rhs + size, b.begin());
// System solution
Vector x(size);
x = 0.0;
if (linsolver_save_system_)
{
// 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);
std::copy(b.begin(), b.end(),
std::ostream_iterator<VectorBlockType>(rhsf, "\n"));
}
int maxit = linsolver_max_iterations_; int maxit = linsolver_max_iterations_;
if (maxit == 0) { if (maxit == 0) {
maxit = 5000; maxit = 5000;
} }
#if HAVE_MPI #if HAVE_MPI
if(comm.type()==typeid(ParallelISTLInformation)) if(comm.type()==typeid(ParallelISTLInformation))
{ {
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm; typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
@ -187,7 +169,7 @@ namespace Opm
Dune::OverlappingSchwarzOperator<Mat,Vector,Vector, Comm> Dune::OverlappingSchwarzOperator<Mat,Vector,Vector, Comm>
opA(A, istlComm); opA(A, istlComm);
Dune::OverlappingSchwarzScalarProduct<Vector,Comm> sp(istlComm); Dune::OverlappingSchwarzScalarProduct<Vector,Comm> sp(istlComm);
return solveSystem(opA, x, solution, b, sp, istlComm, maxit); return solveSystem(opA, solution, rhs, sp, istlComm, maxit);
} }
else else
#endif #endif
@ -195,15 +177,34 @@ namespace Opm
Dune::SeqScalarProduct<Vector> sp; Dune::SeqScalarProduct<Vector> sp;
Dune::Amg::SequentialInformation comm; Dune::Amg::SequentialInformation comm;
Operator opA(A); Operator opA(A);
return solveSystem(opA, x, solution, b, sp, comm, maxit); return solveSystem(opA, solution, rhs, sp, comm, maxit);
} }
} }
template<class O, class V, class S, class C> template<class O, class S, class C>
LinearSolverInterface::LinearSolverReport LinearSolverInterface::LinearSolverReport
LinearSolverIstl::solveSystem (O& opA, V& x, double* solution, V& b, LinearSolverIstl::solveSystem (O& opA, double* solution, const double* rhs,
S& sp, const C& comm, int maxit) const S& sp, const C& comm, int maxit) const
{ {
// System RHS
Vector b(opA.getmat().N());
std::copy(rhs, rhs + size, b.begin());
// System solution
Vector x(opA.getmat().M());
x = 0.0;
if (linsolver_save_system_)
{
// Save system to files.
writeMatrixToMatlab(opA.getmat(), 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);
std::copy(b.begin(), b.end(),
std::ostream_iterator<VectorBlockType>(rhsf, "\n"));
}
LinearSolverReport res; LinearSolverReport res;
switch (linsolver_type_) { switch (linsolver_type_) {
case CG_ILU0: case CG_ILU0:

10
opm/core/linalg/LinearSolverIstl.hpp
View File

@ -88,16 +88,14 @@ namespace Opm
private: private:
/// \brief Solve the linear system using ISTL /// \brief Solve the linear system using ISTL
/// \param[in] opA The linear operator of the system to solve. /// \param[in] opA The linear operator of the system to solve.
/// \param[in,out] x The vector with the initial guess that will store
/// the computed solution
/// \param[out] solution C array for storing the solution vector. /// \param[out] solution C array for storing the solution vector.
/// \param[in] b The vector containing the right hand side of the system. /// \param[in] rhs C array containing the right hand side.
/// \param[in] sp The scalar product to use. /// \param[in] sp The scalar product to use.
/// \param[in] comm The information about the parallel domain decomposition. /// \param[in] comm The information about the parallel domain decomposition.
/// \param[in] maxit The maximum number of iterations allowed. /// \param[in] maxit The maximum number of iterations allowed.
template<class O, class V, class S, class C> template<class O, class S, class C>
LinearSolverReport solveSystem(O& opA, V& x, double* solution, LinearSolverReport solveSystem(O& opA, double* solution, const double *rhs,
V& b, S& sp, const C& comm, int maxit) const; S& sp, const C& comm, int maxit) const;
double linsolver_residual_tolerance_; double linsolver_residual_tolerance_;
int linsolver_verbosity_; int linsolver_verbosity_;