Merge pull request #314 from rolk/314_agmg

Remove support for AGMG
This commit is contained in:
Bård Skaflestad committed 2013-08-14 01:40:39 -07:00
commit 9b313c2d2e
10 files changed
-529

No files matched your search

-11
View File
@@ -47,14 +47,6 @@ endmacro (prereqs_hook)
macro (sources_hook)
# these solvers are only compiled in if their dependency is found
if (NOT AGMG_FOUND)
list (REMOVE_ITEM opm-core_SOURCES
${PROJECT_SOURCE_DIR}/${opm-core_DIR}/core/linalg/LinearSolverAGMG.cpp
)
list (REMOVE_ITEM attic_SOURCES
${PROJECT_SOURCE_DIR}/attic/test_agmg.cpp
)
endif (NOT AGMG_FOUND)
if (NOT dune-istl_FOUND)
list (REMOVE_ITEM opm-core_SOURCES
${PROJECT_SOURCE_DIR}/${opm-core_DIR}/core/linalg/LinearSolverIstl.cpp
@@ -100,12 +92,9 @@ macro (sources_hook)
endmacro (sources_hook)
macro (fortran_hook)
# only include Fortran support if AGMG sources are available
set (${project}_FORTRAN_IF HAVE_AGMG)
endmacro (fortran_hook)
macro (tests_hook)
cond_disable_test ("AGMG")
cond_disable_test ("ERT")
endmacro (tests_hook)
-3
View File
@@ -42,7 +42,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/io/eclipse/writeECLData.cpp
opm/core/io/vag/vag.cpp
opm/core/io/vtk/writeVtkData.cpp
opm/core/linalg/LinearSolverAGMG.cpp
opm/core/linalg/LinearSolverFactory.cpp
opm/core/linalg/LinearSolverInterface.cpp
opm/core/linalg/LinearSolverIstl.cpp
@@ -187,7 +186,6 @@ list (APPEND ATTIC_FILES
attic/bo_resprop_test.cpp
attic/pvt_test.cpp
attic/relperm_test.cpp
attic/test_agmg.cpp
attic/test_cfs_tpfa.c
attic/test_ert.cpp
attic/test_jacsys.cpp
@@ -230,7 +228,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/io/eclipse/writeECLData.hpp
opm/core/io/vag/vag.hpp
opm/core/io/vtk/writeVtkData.hpp
opm/core/linalg/LinearSolverAGMG.hpp
opm/core/linalg/LinearSolverFactory.hpp
opm/core/linalg/LinearSolverInterface.hpp
opm/core/linalg/LinearSolverIstl.hpp
-125
View File
@@ -1,125 +0,0 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
Copyright 2012 Statoil ASA.
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 "config.h"
#include <cassert>
#include <cmath>
#include <cstddef>
#include <iomanip>
#include <iostream>
#include <vector>
#include <memory>
#include <opm/core/linalg/sparse_sys.h>
#include <opm/core/linalg/LinearSolverAGMG.hpp>
namespace {
std::size_t
compute_nnz(const std::size_t m)
{
assert (m > 0);
std::size_t nnz = m; // A(i,i)
nnz += (m > 1) ? 2 : 0; // A(0,1), A(m-1,m-2)
nnz += (m > 2) ? 2 * (m - 2) : 0; // A(i,i-1), A(i,i+1)
return nnz;
}
std::shared_ptr<CSRMatrix>
build_laplace_1d(const std::size_t m)
{
assert (m >= 2);
const std::size_t nnz = compute_nnz(m);
std::shared_ptr<CSRMatrix>
A(csrmatrix_new_known_nnz(m, nnz), csrmatrix_delete);
A->ia[ 0 ] = 0;
// First row
A->ia[ 0 + 1 ] = A->ia[ 0 ];
A->ja[ A->ia[0 + 1] ] = 0 + 0;
A->sa[ A->ia[0 + 1] ++ ] = 2.0;
A->ja[ A->ia[0 + 1] ] = 0 + 1;
A->sa[ A->ia[0 + 1] ++ ] = - 1.0;
// General rows
for (std::size_t i = 1; i < m - 1; ++i) {
A->ia[i + 1] = A->ia[i];
A->ja[ A->ia[i + 1] ] = int(i) - 1;
A->sa[ A->ia[i + 1] ++ ] = - 1.0;
A->ja[ A->ia[i + 1] ] = int(i) ;
A->sa[ A->ia[i + 1] ++ ] = 2.0;
A->ja[ A->ia[i + 1] ] = int(i) + 1;
A->sa[ A->ia[i + 1] ++ ] = - 1.0;
}
// Last row
A->ia[ (m - 1) + 1 ] = A->ia[ m - 1 ];
A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) - 1;
A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = - 1.0;
A->ja[ A->ia[ (m - 1) + 1 ] ] = int(m - 1) ;
A->sa[ A->ia[ (m - 1) + 1 ] ++ ] = 2.0;
return A;
}
}
int main()
{
const std::size_t m = 10;
std::shared_ptr<CSRMatrix> A = build_laplace_1d(m);
// Form right-hand side [1, 0, 0, ...., 0, 1]
std::vector<double> b(m, 0.0);
b[0] = 1.0; b.back() = 1.0;
// Allocate solution vector
std::vector<double> x(m);
// Create solver for SPD system.
Opm::LinearSolverAGMG linsolve(100, 1e-9, true);
linsolve.solve(A.get(), & b[0], & x[0]);
double e = 0.0;
for (std::size_t i = 0; i < m; ++i) {
const double d = x[i] - 1.0;
e += d * d;
}
std::cerr << "|| e ||_2 = "
<< std::scientific
<< std::setprecision(5)
<< std::sqrt(e) / double(m) << '\n';
return 0;
}
-94
View File
@@ -1,94 +0,0 @@
# - Find Notay's Algebraic Multigrid Solver
#
# Set the path to the source directory of AGMG in the cache variable
# AGMG_ROOT.
#
# Note the difference between AGMG_DIR and AGMG_ROOT. The former will
# cause find_package to switch to config mode and search for a file
# named agmg-config.cmake, thereby bypassing this module altogether,
# whereas the latter communicates the location of the library to this
# module.
#
# When found, add the contents of AGMG_SOURCES to your own list of
# sources to compile and link for the target.
#
# Use define_fc_func from UseFortranWrappers to write FC_FUNC to your
# own config.h, and declare the function dagmg using this macro.
# USE_MPI is an option that must be declared in the program
# if this is enabled, then we use the parallel version of MUMPS
# in package "libmumps-dev"; otherwise use serial version in
# "libmumps-seq-dev"
# use the same options when searching for MPI
if (AGMG_FIND_REQUIRED)
set (_required "REQUIRED")
endif ()
if (AGMG_FIND_QUIETLY)
set (_quiet "QUIET")
endif ()
# use different version of library depending on MPI or not
if (USE_MPI)
set (_seq "")
set (_par "_par")
find_package (MPI ${_required} ${_quiet})
set (_mpi_found "MPI_FOUND")
else ()
set (_seq "_seq")
set (_par "")
set (_mpi_found "")
endif ()
find_file (AGMG_SOURCES
dagmg${_par}.f90
PATHS ${AGMG_ROOT}
PATH_SUFFIXES SRC
DOC "Yvan Notay's Algebraic Multigrid Solver, Double Precision version"
NO_DEFAULT_PATH
)
# AGMG is dependent on having the MUMPS library present
find_path (MUMPS_INCLUDE_DIR
dmumps_struc.h
PATH_SUFFIXES include
)
find_library (MUMPS_LIBRARY
NAMES dmumps${_seq}
DOC "MUltifrontal Massive Parallel sparse direct Solver"
)
# make sure that we can compile Fortran code
if (AGMG_SOURCES)
enable_language (Fortran)
endif (AGMG_SOURCES)
# set value for config.h
if (AGMG_SOURCES)
set (HAVE_AGMG 1 CACHE INT "Is AGMG present?")
else (AGMG_SOURCES)
unset (HAVE_AGMG CACHE)
endif (AGMG_SOURCES)
# handle REQUIRED and QUIET standard options
include (FindPackageHandleStandardArgs)
find_package_handle_standard_args (AGMG
DEFAULT_MSG
AGMG_SOURCES
MUMPS_INCLUDE_DIR
MUMPS_LIBRARY
CMAKE_Fortran_COMPILER_SUPPORTS_F90
${_mpi_found}
)
# add our own compatibility routine to link with system MUMPS
if (AGMG_SOURCES)
list (APPEND AGMG_SOURCES "${CMAKE_CURRENT_SOURCE_DIR}/cmake/Templates/dagmg_mumps.f90")
list (APPEND AGMG_INCLUDE_DIRS "${MUMPS_INCLUDE_DIR}")
list (APPEND AGMG_LIBRARIES "${MUMPS_LIBRARY}")
endif ()
if (USE_MPI AND MPI_FOUND)
list (APPEND AGMG_INCLUDE_DIRS "${MPI_Fortran_INCLUDE_PATH}")
list (APPEND AGMG_LIBRARIES "${MPI_Fortran_LIBRARIES}")
endif ()
-3
View File
@@ -3,7 +3,6 @@
# defines that must be present in config.h for our headers
set (opm-core_CONFIG_VAR
HAVE_AGMG
HAVE_DUNE_ISTL
HAVE_DYNAMIC_BOOST_TEST
HAVE_ERT
@@ -28,8 +27,6 @@ set (opm-core_DEPS
"SuiteSparse COMPONENTS umfpack"
# solver
"SuperLU"
# Algebraic Multigrid
"AGMG"
# xml processing (for config parsing)
"TinyXML"
# Ensembles-based Reservoir Tools (ERT)
-6
View File
@@ -39,10 +39,6 @@ Optional Packages:
--with-opm=PATH specify parent of all OPM modules not specified
--with-opm-MODULE=PATH use given OPM module from a specified location
--with-superlu=PATH user defined path to SuperLU library
--with-agmg=PATH Include DOUBLE PRECISION version Notay's of AGMG
Algebraic Multigrid solver from specified source
location. Note: this option requires a complete,
working Fortran 90 environment.
--with-umfpack=PATH use UMFPACK/SuiteSparse from a specified location
--with-ert=PATH Use ERT libraries
--with-tinyxml=PATH use TinyXML library from a specified location
@@ -187,7 +183,6 @@ for OPT in "$@"; do
mpi-prefix)
rootvar="_MPI_PREFIX_PATH"
;;
agmg |\
ert |\
boost |\
superlu |\
@@ -280,7 +275,6 @@ for OPT in "$@"; do
use_ninja=""
pkgname=""
;;
agmg |\
ert |\
superlu)
pkgname="${pkgname^^}"
-11
View File
@@ -1,11 +0,0 @@
! AGMG expects an interface to MUMPS named like this
subroutine DAGMG_MUMPS (id)
! In here are the official defines for the parameter passed to
! MUMPS. These must match the corresponding ones in dagmg.f90
include 'dmumps_struc.h'
type (DMUMPS_STRUC) :: id
! Let the system library do all the heavy lifting
call DMUMPS (id)
end subroutine DAGMG_MUMPS
-150
View File
@@ -1,150 +0,0 @@
/*
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.
#ifdef HAVE_MPI
#define DAGMG_ FC_FUNC(dagmgpar, DAGMGPAR)
#else
#define DAGMG_ FC_FUNC(dagmg, DAGMG)
#endif
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)
{
}
LinearSolverAGMG::LinearSolverAGMG(const parameter::ParameterGroup& param)
: max_it_(100) ,
rtol_ (1.0e-8),
is_spd_(false),
linsolver_verbosity_(0)
{
max_it_ = param.getDefault("max_it", max_it_);
rtol_ = param.getDefault("rtol" , rtol_ );
is_spd_ = param.getDefault("is_spd", is_spd_);
linsolver_verbosity_ = param.getDefault("linsolver_verbosity",
linsolver_verbosity_);
}
LinearSolverAGMG::~LinearSolverAGMG() {}
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 = linsolver_verbosity_?6:-1; // 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;
}
void
LinearSolverAGMG::setTolerance(const double tol)
{
rtol_ = tol;
}
double
LinearSolverAGMG::getTolerance() const
{
return rtol_;
}
} // namespace Opm
-115
View File
@@ -1,115 +0,0 @@
/*
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>
#include <opm/core/utility/parameters/ParameterGroup.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-8,
const bool is_spd = false);
/**
* Constructor.
* \param[in] param ParameterGroup object containing the fields
* max_it,rtol,is_spd as used in the constructor.
*/
LinearSolverAGMG(const parameter::ParameterGroup& param);
/**
* 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 columns).
/// \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[in,out] 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;
/// Set tolerance for the linear solver.
/// \param[in] tol tolerance value
virtual void setTolerance(const double tol);
/// Get tolerance for the linear solver.
/// \return tolerance value
virtual double getTolerance() const;
private:
int max_it_;
double rtol_ ;
bool is_spd_;
int linsolver_verbosity_;
};
} // namespace Opm
#endif // OPM_LINEARSOLVERAGMG_HEADER_INCLUDED
-11
View File
@@ -31,11 +31,6 @@
#include <opm/core/linalg/LinearSolverIstl.hpp>
#endif
#if HAVE_AGMG
#include <opm/core/linalg/LinearSolverAGMG.hpp>
#endif
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <string>
@@ -75,12 +70,6 @@ namespace Opm
#endif
}
else if (ls == "agmg") {
#if HAVE_AGMG
solver_.reset(new LinearSolverAGMG(param));
#endif
}
else {
THROW("Linear solver " << ls << " is unknown.");
}