Add simple gateway to Notay's AGMG solver package.
Hook up to build, but don't expose through the LinearSolverFactory. Only build tested. Won't be built unless the user explicitly enables AGMG support.
This commit is contained in:
Bård Skaflestad
parent
496be60220
commit
ca7b4a5da3
@ -276,7 +276,11 @@ endif
|
||||
if BUILD_AGMG
|
||||
libopmcore_la_SOURCES += \
|
||||
$(AGMG_SRCDIR)/dagmg.f90 \
|
||||
$(AGMG_SRCDIR)/dagmg_mumps.f90
|
||||
$(AGMG_SRCDIR)/dagmg_mumps.f90 \
|
||||
opm/core/linalg/LinearSolverAGMG.cpp
|
||||
|
||||
nobase_include_HEADERS += \
|
||||
opm/core/linalg/LinearSolverAGMG.hpp
|
||||
|
||||
libopmcore_la_LDFLAGS += \
|
||||
$(FCLIBS)
|
||||
|
||||
122
opm/core/linalg/LinearSolverAGMG.cpp
Normal file
122
opm/core/linalg/LinearSolverAGMG.cpp
Normal file
@ -0,0 +1,122 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
|
||||
#if HAVE_CONFIG_H
|
||||
#include <config.h>
|
||||
#endif
|
||||
|
||||
#include <opm/core/linalg/LinearSolverAGMG.hpp>
|
||||
#include <opm/core/utility/ErrorMacros.hpp>
|
||||
|
||||
#include <algorithm>
|
||||
#include <cassert>
|
||||
#include <iterator>
|
||||
#include <stdexcept>
|
||||
#include <vector>
|
||||
|
||||
#if HAVE_AGMG
|
||||
// Manual prototype of main gateway routine to DOUBLE PRECISION,
|
||||
// serial version of the AGMG software package.
|
||||
//
|
||||
// Note that both the matrix entries and column indices are writable.
|
||||
// The solver may permute the matrix entries within each row during
|
||||
// the setup phase.
|
||||
#define DAGMG_ F77_FUNC(dagmg, DAGMG)
|
||||
|
||||
extern "C"
|
||||
void
|
||||
DAGMG_(const int* N , // System size
|
||||
double* sa , // Non-zero elements
|
||||
int* ja , // Column indices
|
||||
const int* ia , // Row pointers
|
||||
const double* f , // Right-hand side
|
||||
double* x , // Solution
|
||||
int* ijob , // Main control parameter
|
||||
int* iprint, // Message output unit
|
||||
int* nrest , // Number of GCR restarts
|
||||
int* iter , // Maximum (and actual) number of iterations
|
||||
const double* tol ); // Residual reduction tolerance
|
||||
|
||||
#endif // HAVE_AGMG
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
LinearSolverAGMG::LinearSolverAGMG(const int max_it,
|
||||
const double rtol ,
|
||||
const bool is_spd)
|
||||
: max_it_(max_it),
|
||||
rtol_ (rtol) ,
|
||||
is_spd_(is_spd)
|
||||
{
|
||||
#if !HAVE_AGMG
|
||||
THROW("AGMG support is not enabled in this library");
|
||||
#endif // HAVE_AGMG
|
||||
}
|
||||
|
||||
LinearSolverInterface::LinearSolverReport
|
||||
LinearSolverAGMG::solve(const int size ,
|
||||
const int nonzeros,
|
||||
const int* ia ,
|
||||
const int* ja ,
|
||||
const double* sa ,
|
||||
const double* rhs ,
|
||||
double* solution) const
|
||||
{
|
||||
const std::vector<double>::size_type nnz = ia[size];
|
||||
|
||||
assert (nnz == std::vector<double>::size_type(nonzeros));
|
||||
|
||||
#if defined(NDEBUG)
|
||||
// Suppress warning about unused parameter.
|
||||
static_cast<void>(nonzeros);
|
||||
#endif
|
||||
|
||||
std::vector<double> a(sa, sa + nnz);
|
||||
|
||||
// Account for 1-based indexing.
|
||||
std::vector<int> i(ia, ia + std::vector<int>::size_type(size + 1));
|
||||
std::transform(i.begin(), i.end(), i.begin(),
|
||||
std::bind2nd(std::plus<int>(), 1));
|
||||
|
||||
std::vector<int> j(ja, ja + nnz);
|
||||
std::transform(j.begin(), j.end(), j.begin(),
|
||||
std::bind2nd(std::plus<int>(), 1));
|
||||
|
||||
LinearSolverInterface::LinearSolverReport rpt = {};
|
||||
rpt.iterations = max_it_;
|
||||
|
||||
int ijob = 0; // Setup + solution + cleanup, x0==0.
|
||||
|
||||
int nrest;
|
||||
if (is_spd_) {
|
||||
nrest = 1; // Use CG algorithm
|
||||
}
|
||||
else {
|
||||
nrest = 10; // Suggested default number of GCR restarts.
|
||||
}
|
||||
|
||||
int iprint = 0; // Suppress most output
|
||||
DAGMG_(& size, & a[0], & j[0], & i[0], rhs, solution,
|
||||
& ijob, & iprint, & nrest, & rpt.iterations, & rtol_);
|
||||
|
||||
rpt.converged = rpt.iterations <= max_it_;
|
||||
return rpt;
|
||||
}
|
||||
} // namespace Opm
|
||||
96
opm/core/linalg/LinearSolverAGMG.hpp
Normal file
96
opm/core/linalg/LinearSolverAGMG.hpp
Normal file
@ -0,0 +1,96 @@
|
||||
/*
|
||||
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/>.
|
||||
*/
|
||||
|
||||
#ifndef OPM_LINEARSOLVERAGMG_HEADER_INCLUDED
|
||||
#define OPM_LINEARSOLVERAGMG_HEADER_INCLUDED
|
||||
|
||||
/**
|
||||
* \file
|
||||
* Gateway to Notay's AGMG package implementing an aggregation-based
|
||||
* algebraic multigrid method.
|
||||
*
|
||||
* References:
|
||||
* * Y. Notay,
|
||||
* An aggregation-based algebraic multigrid method,
|
||||
* Electronic Transactions on Numerical Analysis, vol 37,
|
||||
* pp. 123-146, 2010.
|
||||
*
|
||||
* * A. Napov and Y. Notay,
|
||||
* An algebraic multigrid method with guaranteed convergence rate,
|
||||
* Report GANMN 10-03, Université Libre de Bruxelles, Brussels,
|
||||
* Belgium, 2010 (Revised 2011).
|
||||
*
|
||||
* * Y. Notay,
|
||||
* Aggregation-based algebraic multigrid for convection-diffusion
|
||||
* equations, Report GANMN 11-01, Université Libre de Bruxelles,
|
||||
* Brussels, Belgium, 2011.
|
||||
*/
|
||||
|
||||
#include <opm/core/linalg/LinearSolverInterface.hpp>
|
||||
|
||||
namespace Opm
|
||||
{
|
||||
/// Concrete class encapsulating Notay's AGMG package
|
||||
class LinearSolverAGMG : public LinearSolverInterface
|
||||
{
|
||||
public:
|
||||
/**
|
||||
* Constructor.
|
||||
* \param[in] max_it Maximum number of solver iterations.
|
||||
* \param[in] rtol Residual reduction tolerance.
|
||||
* \param[in] is_spd Whether or not the matrix is SPD. SPD
|
||||
* systems are solved using CG while other
|
||||
* systems are solved using GCR.
|
||||
*/
|
||||
LinearSolverAGMG(const int max_it = 100 ,
|
||||
const double rtol = 1.0e-6,
|
||||
const bool is_spd = false);
|
||||
|
||||
/// Destructor.
|
||||
virtual ~LinearSolverAGMG();
|
||||
|
||||
using LinearSolverInterface::solve;
|
||||
|
||||
/// Solve a linear system, with a matrix given in compressed
|
||||
/// sparse row format.
|
||||
/// \param[in] size Number of rows (and colums).
|
||||
/// \param[in] nonzeros Number of (structural) non-zeros.
|
||||
/// \param[in] ia Row pointers.
|
||||
/// \param[in] ja Column indices.
|
||||
/// \param[in] sa (structurally) non-zero elements.
|
||||
/// \param[in] rhs System right-hand side.
|
||||
/// \param[inout] solution System solution.
|
||||
/// \return Solver meta-data concerning most recent system solve.
|
||||
virtual LinearSolverInterface::LinearSolverReport
|
||||
solve(const int size, const int nonzeros,
|
||||
const int* ia, const int* ja, const double* sa,
|
||||
const double* rhs, double* solution) const;
|
||||
|
||||
private:
|
||||
int max_it_;
|
||||
double rtol_ ;
|
||||
bool is_spd_;
|
||||
};
|
||||
|
||||
|
||||
} // namespace Opm
|
||||
|
||||
|
||||
|
||||
#endif // OPM_LINEARSOLVERISTL_HEADER_INCLUDED
|
||||
Loading…
Reference in New Issue
Block a user