added: support for dune-istl linear solvers
This commit is contained in:
Arne Morten Kvarving
committed by
Knut Morten Okstad
Knut Morten Okstad
parent
0364126769
commit
efa6c99665
@@ -110,6 +110,12 @@ if(NOT PETSC_FOUND)
|
||||
${PROJECT_SOURCE_DIR}/src/LinAlg/PETScSolParams.C)
|
||||
endif()
|
||||
|
||||
if(NOT ISTL_FOUND)
|
||||
list(REMOVE_ITEM IFEM_SRCS
|
||||
${PROJECT_SOURCE_DIR}/src/LinAlg/ISTLMatrix.C
|
||||
${PROJECT_SOURCE_DIR}/src/LinAlg/ISTLSolParams.C)
|
||||
endif()
|
||||
|
||||
if(LRSPLINE_FOUND OR LRSpline_FOUND)
|
||||
file(GLOB LR_SRCS ${PROJECT_SOURCE_DIR}/src/ASM/LR/*.C)
|
||||
list(APPEND IFEM_SRCS ${LR_SRCS})
|
||||
|
||||
@@ -211,6 +211,16 @@ IF(IFEM_USE_OPENMP)
|
||||
ENDIF(OPENMP_FOUND)
|
||||
ENDIF(IFEM_USE_OPENMP)
|
||||
|
||||
# ISTL
|
||||
if(IFEM_USE_ISTL)
|
||||
find_package(ISTL)
|
||||
if(ISTL_FOUND)
|
||||
list(APPEND IFEM_DEPLIBS ${ISTL_LIBRARIES})
|
||||
list(APPEND IFEM_DEPINCLUDES ${ISTL_INCLUDE_DIRS})
|
||||
list(APPEND IFEM_DEFINITIONS -DHAS_ISTL=1 ${ISTL_DEFINITIONS})
|
||||
endif()
|
||||
endif()
|
||||
|
||||
# Portability issues
|
||||
include(CheckFunctionExists)
|
||||
set(CMAKE_REQUIRED_DEFINITIONS)
|
||||
|
||||
38
cmake/Modules/FindISTL.cmake
Normal file
38
cmake/Modules/FindISTL.cmake
Normal file
@@ -0,0 +1,38 @@
|
||||
find_package(PkgConfig)
|
||||
|
||||
set(OLD_PKG $ENV{PKG_CONFIG_PATH})
|
||||
set(ENV{PKG_CONFIG_PATH} $ENV{PETSC_DIR}/$ENV{PETSC_ARCH}/lib/pkgconfig)
|
||||
pkg_check_modules(ISTL dune-istl)
|
||||
set(ENV{PKG_CONFIG_PATH} ${OLD_PKG})
|
||||
|
||||
list(APPEND ISTL_DEFINITIONS -DHAVE_NULLPTR=1 -DISTL_VERSION="${ISTL_VERSION}")
|
||||
string(REPLACE "." ";" VERSION_LIST "${ISTL_VERSION}")
|
||||
if(VERSION_LIST)
|
||||
list(GET VERSION_LIST 0 DUNE_ISTL_VERSION_MAJOR)
|
||||
list(GET VERSION_LIST 1 DUNE_ISTL_VERSION_MINOR)
|
||||
list(GET VERSION_LIST 2 DUNE_ISTL_VERSION_PATCH)
|
||||
list(APPEND ISTL_DEFINITIONS -DDUNE_ISTL_VERSION_MAJOR=${DUNE_ISTL_VERSION_MAJOR}
|
||||
-DDUNE_ISTL_VERSION_MINOR=${DUNE_ISTL_VERSION_MINOR}
|
||||
-DDUNE_ISTL_VERSION_REVISION=${DUNE_ISTL_VERSION_PATCH})
|
||||
endif()
|
||||
|
||||
if(SUPERLU_FOUND)
|
||||
list(APPEND ISTL_DEFINITIONS -DHAVE_SUPERLU=1 -DSUPERLU_POST_2005_VERSION=1 -DSUPERLU_NTYPE=1 -DSUPERLU_MIN_VERSION_4_3=1)
|
||||
endif()
|
||||
|
||||
if(IFEM_USE_UMFPACK)
|
||||
find_package(SuiteSparse COMPONENTS umfpack)
|
||||
endif()
|
||||
|
||||
if(SuiteSparse_UMFPACK_FOUND)
|
||||
list(APPEND ISTL_DEFINITIONS -DHAVE_UMFPACK=1)
|
||||
list(APPEND ISTL_INCLUDE_DIRS ${UMFPACK_INCLUDE_DIR})
|
||||
list(APPEND ISTL_LIBRARIES ${UMFPACK_LIBRARY})
|
||||
endif()
|
||||
|
||||
|
||||
list(APPEND ISTL_INCLUDE_DIRS ${ISTL_INCLUDEDIR})
|
||||
|
||||
include(FindPackageHandleStandardArgs)
|
||||
find_package_handle_standard_args(ISTL DEFAULT_MSG
|
||||
ISTL_INCLUDE_DIRS ISTL_LIBRARIES)
|
||||
294
cmake/Modules/FindSuiteSparse.cmake
Normal file
294
cmake/Modules/FindSuiteSparse.cmake
Normal file
@@ -0,0 +1,294 @@
|
||||
# - Find Tim Davis' SuiteSparse collection of sparse matrix libraries
|
||||
#
|
||||
# Synopsis:
|
||||
# find_package (SuiteSparse COMPONENTS <list-of-components>)
|
||||
#
|
||||
# Components are:
|
||||
# amd Approximate Minimum Degree ordering
|
||||
# camd Constrained Approximate Minimum Degree ordering
|
||||
# colamd COLumn Approximate Minimum Degree ordering
|
||||
# ccolamd Constrained COLumn Approximate Minimum Degree ordering
|
||||
# cholmod Supernodal sparse Cholesky factorization and update
|
||||
# umfpack Unsymmetric MultiFrontal sparse LU factorization
|
||||
#
|
||||
# The following variables will be set:
|
||||
#
|
||||
# SuiteSparse_FOUND True if all dependencies are satisfied
|
||||
# SuiteSparse_Xxx_FOUND True if module Xxx is found
|
||||
# HAVE_SUITESPARSE_Xxx_H Binary value indicating presence of header
|
||||
# SuiteSparse_INCLUDE_DIRS Paths containing the SuiteSparse header files
|
||||
# SuiteSparse_LIBRARIES Name of the libraries which must be linked
|
||||
# SuiteSparse_DEFINITIONS Defines that must be passed to the compiler
|
||||
# SuiteSparse_LINKER_FLAGS Options that must be passed when linking
|
||||
#
|
||||
# The following options can be set to configure the module:
|
||||
#
|
||||
# SUITESPARSE_USE_STATIC Link with a static library, even if a
|
||||
# dynamic library is also present. Note that
|
||||
# setting this to OFF does not ensure that a
|
||||
# shared library will be used.
|
||||
#
|
||||
# See <http://www.cise.ufl.edu/research/sparse/SuiteSparse>.
|
||||
|
||||
# Copyright (C) 2012 Uni Research AS
|
||||
# This file is licensed under the GNU General Public License v3.0
|
||||
|
||||
function (try_compile_umfpack varname)
|
||||
include (CMakePushCheckState)
|
||||
include (CheckCSourceCompiles)
|
||||
cmake_push_check_state ()
|
||||
set (CMAKE_REQUIRED_INCLUDES ${UMFPACK_INCLUDE_DIRS})
|
||||
set (CMAKE_REQUIRED_LIBRARIES ${UMFPACK_LIBRARY} ${ARGN} ${SuiteSparse_EXTRA_LIBS})
|
||||
check_c_source_compiles (
|
||||
"#include <umfpack.h>
|
||||
int main (void) {
|
||||
void *Symbolic, *Numeric;
|
||||
double Info[UMFPACK_INFO], Control[UMFPACK_CONTROL];
|
||||
umfpack_dl_defaults(Control);
|
||||
umfpack_dl_symbolic(0, 0, 0, 0, 0,
|
||||
&Symbolic, Control, Info);
|
||||
umfpack_dl_numeric (0, 0, 0,
|
||||
Symbolic, &Numeric, Control, Info);
|
||||
umfpack_dl_free_symbolic(&Symbolic);
|
||||
umfpack_dl_solve(UMFPACK_A, 0, 0, 0, 0, 0,
|
||||
Numeric, Control, Info);
|
||||
umfpack_dl_free_numeric(&Numeric);
|
||||
umfpack_timer ();
|
||||
return 0;
|
||||
}" ${varname})
|
||||
cmake_pop_check_state ()
|
||||
set (${varname} "${${varname}}" PARENT_SCOPE)
|
||||
endfunction (try_compile_umfpack varname)
|
||||
|
||||
# variables to pass on to other packages
|
||||
if (FIND_QUIETLY)
|
||||
set (SuiteSparse_QUIET "QUIET")
|
||||
else (FIND_QUIETLY)
|
||||
set (SuiteSparse_QUIET "")
|
||||
endif (FIND_QUIETLY)
|
||||
|
||||
# we need to link to BLAS and LAPACK
|
||||
if (NOT BLAS_FOUND)
|
||||
find_package (BLAS ${SuiteSparse_QUIET} REQUIRED)
|
||||
endif (NOT BLAS_FOUND)
|
||||
if (NOT LAPACK_FOUND)
|
||||
find_package (LAPACK ${SuiteSparse_QUIET} REQUIRED)
|
||||
endif (NOT LAPACK_FOUND)
|
||||
|
||||
# we also need the math part of the runtime library
|
||||
find_library (MATH_LIBRARY NAMES "m")
|
||||
set (SuiteSparse_EXTRA_LIBS ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${MATH_LIBRARY})
|
||||
|
||||
# if we don't get any further clues about where to look, then start
|
||||
# roaming around the system
|
||||
set (_no_default_path "")
|
||||
|
||||
# search system directories by default
|
||||
set (SuiteSparse_SEARCH_PATH)
|
||||
|
||||
# pick up paths from the environment if specified there; these replace the
|
||||
# pre-defined paths so that we don't accidentially pick up old stuff
|
||||
if (NOT $ENV{SuiteSparse_DIR} STREQUAL "")
|
||||
set (SuiteSparse_SEARCH_PATH "$ENV{SuiteSparse_DIR}")
|
||||
endif (NOT $ENV{SuiteSparse_DIR} STREQUAL "")
|
||||
if (SuiteSparse_DIR)
|
||||
set (SuiteSparse_SEARCH_PATH "${SuiteSparse_DIR}")
|
||||
endif (SuiteSparse_DIR)
|
||||
# CMake uses _DIR suffix as default for config-mode files; it is unlikely
|
||||
# that we are building SuiteSparse ourselves; use _ROOT suffix to specify
|
||||
# location to pre-canned binaries
|
||||
if (NOT $ENV{SuiteSparse_ROOT} STREQUAL "")
|
||||
set (SuiteSparse_SEARCH_PATH "$ENV{SuiteSparse_ROOT}")
|
||||
endif (NOT $ENV{SuiteSparse_ROOT} STREQUAL "")
|
||||
if (SuiteSparse_ROOT)
|
||||
set (SuiteSparse_SEARCH_PATH "${SuiteSparse_ROOT}")
|
||||
endif (SuiteSparse_ROOT)
|
||||
# most commonly, we use the uppercase version of this variable
|
||||
if (SUITESPARSE_ROOT)
|
||||
set (SuiteSparse_SEARCH_PATH "${SUITESPARSE_ROOT}")
|
||||
endif (SUITESPARSE_ROOT)
|
||||
|
||||
# if we have specified a search path, then confine ourselves to that
|
||||
if (SuiteSparse_SEARCH_PATH)
|
||||
set (_no_default_path "NO_DEFAULT_PATH")
|
||||
endif (SuiteSparse_SEARCH_PATH)
|
||||
|
||||
# transitive closure of dependencies; after this SuiteSparse_MODULES is the
|
||||
# full list of modules that must be found to satisfy the user's link demands
|
||||
set (SuiteSparse_MODULES ${SuiteSparse_FIND_COMPONENTS})
|
||||
list (FIND SuiteSparse_MODULES "umfpack" UMFPACK_DESIRED)
|
||||
if (NOT UMFPACK_DESIRED EQUAL -1)
|
||||
list (APPEND SuiteSparse_MODULES amd cholmod)
|
||||
endif (NOT UMFPACK_DESIRED EQUAL -1)
|
||||
list (FIND SuiteSparse_MODULES "cholmod" CHOLMOD_DESIRED)
|
||||
if (NOT CHOLMOD_DESIRED EQUAL -1)
|
||||
list (APPEND SuiteSparse_MODULES amd camd colamd)
|
||||
endif (NOT CHOLMOD_DESIRED EQUAL -1)
|
||||
if (SuiteSparse_MODULES)
|
||||
list (REMOVE_DUPLICATES SuiteSparse_MODULES)
|
||||
endif (SuiteSparse_MODULES)
|
||||
|
||||
# if someone else already have found all the packages for us, then don't do anything
|
||||
set (SuiteSparse_EVERYTHING_FOUND TRUE)
|
||||
foreach (module IN LISTS SuiteSparse_MODULES)
|
||||
string (TOUPPER ${module} MODULE)
|
||||
if (NOT SuiteSparse_${MODULE}_FOUND)
|
||||
set (SuiteSparse_EVERYTHING_FOUND FALSE)
|
||||
break ()
|
||||
endif (NOT SuiteSparse_${MODULE}_FOUND)
|
||||
endforeach (module)
|
||||
if (SuiteSparse_EVERYTHING_FOUND)
|
||||
return ()
|
||||
endif (SuiteSparse_EVERYTHING_FOUND)
|
||||
|
||||
# only search in architecture-relevant directory
|
||||
if (CMAKE_SIZEOF_VOID_P)
|
||||
math (EXPR _BITS "8 * ${CMAKE_SIZEOF_VOID_P}")
|
||||
endif (CMAKE_SIZEOF_VOID_P)
|
||||
|
||||
# if we are told to link SuiteSparse statically, add these parts
|
||||
# to the name so we always match only that particular type of lib
|
||||
option (SUITESPARSE_USE_STATIC "Link SuiteSparse statically" OFF)
|
||||
mark_as_advanced (SUITESPARSE_USE_STATIC)
|
||||
if (SUITESPARSE_USE_STATIC)
|
||||
set (_pref_ "${CMAKE_STATIC_LIBRARY_PREFIX}")
|
||||
set (_suff_ "${CMAKE_STATIC_LIBRARY_SUFFIX}")
|
||||
else (SUITESPARSE_USE_STATIC)
|
||||
set (_pref_ "")
|
||||
set (_suff_ "")
|
||||
endif (SUITESPARSE_USE_STATIC)
|
||||
|
||||
# if SuiteSparse >= 4.0 we must also link with libsuitesparseconfig
|
||||
# assume that this is the case if we find the library; otherwise just
|
||||
# ignore it (older versions don't have a file named like this)
|
||||
find_library (config_LIBRARY
|
||||
NAMES "${_pref_}suitesparseconfig${_suff_}"
|
||||
PATHS ${SuiteSparse_SEARCH_PATH}
|
||||
PATH_SUFFIXES ".libs" "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}" "lib/ufsparse"
|
||||
${_no_default_path}
|
||||
)
|
||||
if (config_LIBRARY)
|
||||
list (APPEND SuiteSparse_EXTRA_LIBS ${config_LIBRARY})
|
||||
# POSIX.1-2001 REALTIME portion require us to link this library too for
|
||||
# clock_gettime() which is used by suitesparseconfig
|
||||
if ("${CMAKE_SYSTEM_NAME}" MATCHES "Linux")
|
||||
list (APPEND SuiteSparse_EXTRA_LIBS "-lrt")
|
||||
endif ("${CMAKE_SYSTEM_NAME}" MATCHES "Linux")
|
||||
endif (config_LIBRARY)
|
||||
|
||||
# search filesystem for each of the module individually
|
||||
foreach (module IN LISTS SuiteSparse_MODULES)
|
||||
string (TOUPPER ${module} MODULE)
|
||||
# search for files which implements this module
|
||||
find_path (${MODULE}_INCLUDE_DIR
|
||||
NAMES ${module}.h
|
||||
PATHS ${SuiteSparse_SEARCH_PATH}
|
||||
PATH_SUFFIXES "include" "include/suitesparse" "include/ufsparse" "${MODULE}/Include"
|
||||
${_no_default_path}
|
||||
)
|
||||
|
||||
find_library (${MODULE}_LIBRARY
|
||||
NAMES "${_pref_}${module}${_suff_}"
|
||||
PATHS ${SuiteSparse_SEARCH_PATH}
|
||||
PATH_SUFFIXES "lib/.libs" "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}" "lib/ufsparse" "${MODULE}/Lib"
|
||||
${_no_default_path}
|
||||
)
|
||||
# start out by including the module itself; other dependencies will be added later
|
||||
set (${MODULE}_INCLUDE_DIRS ${${MODULE}_INCLUDE_DIR})
|
||||
set (${MODULE}_LIBRARIES ${${MODULE}_LIBRARY})
|
||||
endforeach (module)
|
||||
|
||||
# insert any inter-modular dependencies here
|
||||
if (CHOLMOD_LIBRARY)
|
||||
list (APPEND CHOLMOD_LIBRARIES ${AMD_LIBRARIES} ${COLAMD_LIBRARIES})
|
||||
# optional libraries; don't insert any -NOT_FOUND paths
|
||||
if (CAMD_LIBRARY)
|
||||
list (APPEND CHOLMOD_LIBRARIES ${CAMD_LIBRARIES})
|
||||
endif (CAMD_LIBRARY)
|
||||
if (CCOLAMD_LIBRARY)
|
||||
list (APPEND CHOLMOD_LIBRARIES ${CCOLAMD_LIBRARIES})
|
||||
endif (CCOLAMD_LIBRARY)
|
||||
list (REVERSE CHOLMOD_LIBRARIES)
|
||||
# always remove the *first* library from the list
|
||||
list (REMOVE_DUPLICATES CHOLMOD_LIBRARIES)
|
||||
list (REVERSE CHOLMOD_LIBRARIES)
|
||||
endif (CHOLMOD_LIBRARY)
|
||||
if (UMFPACK_LIBRARY)
|
||||
set (UMFPACK_EXTRA_LIBS)
|
||||
# test if umfpack is usable with only amd and not cholmod
|
||||
try_compile_umfpack (HAVE_UMFPACK_WITHOUT_CHOLMOD ${AMD_LIBRARIES})
|
||||
if (HAVE_UMFPACK_WITHOUT_CHOLMOD)
|
||||
list (APPEND UMFPACK_EXTRA_LIBS ${AMD_LIBRARIES})
|
||||
else (HAVE_UMFPACK_WITHOUT_CHOLMOD)
|
||||
if (CHOLMOD_LIBRARIES)
|
||||
try_compile_umfpack (HAVE_UMFPACK_WITH_CHOLMOD ${CHOLMOD_LIBRARIES})
|
||||
if (HAVE_UMFPACK_WITH_CHOLMOD)
|
||||
list (APPEND UMFPACK_EXTRA_LIBS ${CHOLMOD_LIBRARIES})
|
||||
else (HAVE_UMFPACK_WITH_CHOLMOD)
|
||||
set (UMFPACK_EXTRA_LIBS "-NOTFOUND")
|
||||
endif (HAVE_UMFPACK_WITH_CHOLMOD)
|
||||
else (CHOLMOD_LIBRARIES)
|
||||
# if we don't have cholmod, then we certainly cannot have umfpack with cholmod
|
||||
set (UMFPACK_EXTRA_LIBS "-NOTFOUND")
|
||||
endif (CHOLMOD_LIBRARIES)
|
||||
endif (HAVE_UMFPACK_WITHOUT_CHOLMOD)
|
||||
list (APPEND UMFPACK_LIBRARIES ${UMFPACK_EXTRA_LIBS})
|
||||
list (REVERSE UMFPACK_LIBRARIES)
|
||||
list (REMOVE_DUPLICATES UMFPACK_LIBRARIES)
|
||||
list (REVERSE UMFPACK_LIBRARIES)
|
||||
endif (UMFPACK_LIBRARY)
|
||||
|
||||
# don't reset these sets; if two packages request SuiteSparse with
|
||||
# different modules, we want the sets to be merged
|
||||
#set (SuiteSparse_LIBRARIES "")
|
||||
#set (SuiteSparse_INCLUDE_DIRS "")
|
||||
|
||||
# determine which modules were found based on whether all dependencies
|
||||
# were satisfied; create a list of ALL modules (specified) that was found
|
||||
# (to be included in one swoop in CMakeLists.txt)
|
||||
set (SuiteSparse_FOUND TRUE)
|
||||
foreach (module IN LISTS SuiteSparse_FIND_COMPONENTS)
|
||||
string (TOUPPER ${module} MODULE)
|
||||
set (SuiteSparse_${MODULE}_FOUND TRUE)
|
||||
foreach (file IN LISTS ${MODULE}_INCLUDE_DIRS ${MODULE}_LIBRARIES)
|
||||
if (NOT EXISTS ${file})
|
||||
set (SuiteSparse_${MODULE}_FOUND FALSE)
|
||||
endif (NOT EXISTS ${file})
|
||||
endforeach (file)
|
||||
if (NOT SuiteSparse_${MODULE}_FOUND)
|
||||
set (SuiteSparse_FOUND FALSE)
|
||||
# use empty string instead of zero, so it can be tested with #ifdef
|
||||
# as well as #if in the source code
|
||||
set (HAVE_SUITESPARSE_${MODULE}_H "" CACHE INT "Is ${module} header present?")
|
||||
else (NOT SuiteSparse_${MODULE}_FOUND)
|
||||
set (HAVE_SUITESPARSE_${MODULE}_H 1 CACHE INT "Is ${module} header present?")
|
||||
list (APPEND SuiteSparse_LIBRARIES "${${MODULE}_LIBRARIES}")
|
||||
list (APPEND SuiteSparse_LINKER_FLAGS "${${MODULE}_LINKER_FLAGS}")
|
||||
list (APPEND SuiteSparse_INCLUDE_DIRS "${${MODULE}_INCLUDE_DIRS}")
|
||||
endif (NOT SuiteSparse_${MODULE}_FOUND)
|
||||
mark_as_advanced (HAVE_SUITESPARSE_${MODULE}_H)
|
||||
mark_as_advanced (${MODULE}_INCLUDE_DIR)
|
||||
mark_as_advanced (${MODULE}_LIBRARY)
|
||||
endforeach (module)
|
||||
|
||||
if (SuiteSparse_INCLUDE_DIRS)
|
||||
list (REMOVE_DUPLICATES SuiteSparse_INCLUDE_DIRS)
|
||||
endif (SuiteSparse_INCLUDE_DIRS)
|
||||
if (SuiteSparse_LIBRARIES)
|
||||
list (REVERSE SuiteSparse_LIBRARIES)
|
||||
list (REMOVE_DUPLICATES SuiteSparse_LIBRARIES)
|
||||
list (REVERSE SuiteSparse_LIBRARIES)
|
||||
endif (SuiteSparse_LIBRARIES)
|
||||
|
||||
# print a message to indicate status of this package
|
||||
include (FindPackageHandleStandardArgs)
|
||||
find_package_handle_standard_args (SuiteSparse
|
||||
DEFAULT_MSG
|
||||
SuiteSparse_LIBRARIES
|
||||
SuiteSparse_INCLUDE_DIRS
|
||||
)
|
||||
|
||||
# add these after checking to not pollute the message output (checking for
|
||||
# BLAS and LAPACK is REQUIRED so if they are not found, we'll have failed
|
||||
# already; suitesparseconfig is "optional" anyway)
|
||||
list (APPEND SuiteSparse_LIBRARIES ${SuiteSparse_EXTRA_LIBS})
|
||||
@@ -4,11 +4,13 @@ OPTION(IFEM_USE_SUPERLU "Compile with SuperLU support?" ON)
|
||||
OPTION(IFEM_USE_SUPERLU_MT "Compile with SuperLU_MT support?" OFF)
|
||||
OPTION(IFEM_USE_PETSC "Compile with PETSc support?" ON)
|
||||
OPTION(IFEM_USE_MPI "Compile with MPI support?" OFF)
|
||||
OPTION(IFEM_USE_ISTL "Compile with dune-istl support?" ON)
|
||||
OPTION(IFEM_USE_SLEPC "Compile with SLEPc support?" OFF)
|
||||
OPTION(IFEM_USE_SPR "Compile with SPR support?" OFF)
|
||||
OPTION(IFEM_USE_SAMG "Compile with SAMG support?" OFF)
|
||||
OPTION(IFEM_USE_HDF5 "Compile with HDF5 support?" ON)
|
||||
OPTION(IFEM_USE_VTFWRITER "Compile with VTFWriter support?" ON)
|
||||
OPTION(IFEM_USE_UMFPACK "Compile with UMFPACK support?" ON)
|
||||
OPTION(IFEM_AS_SUBMODULE "Compile IFEM as a submodule of apps?" OFF)
|
||||
OPTION(IFEM_WHOLE_PROG_OPTIM "Compile IFEM with link-time optimizations?" OFF)
|
||||
OPTION(IFEM_TEST_MEMCHECK "Run tests through valgrind?" OFF)
|
||||
|
||||
@@ -53,6 +53,12 @@ macro(IFEM_add_unittests IFEM_PATH)
|
||||
file(GLOB TEST_SOURCES ${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C;${IFEM_PATH}/src/LinAlg/Test/*.C;${IFEM_PATH}/src/SIM/Test/*.C)
|
||||
if(NOT PETSC_FOUND)
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestPETScMatrix.C)
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestISTLPETScMatrix.C)
|
||||
endif()
|
||||
|
||||
if(NOT ISTL_FOUND)
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestISTLMatrix.C)
|
||||
list(REMOVE_ITEM TEST_SOURCES ${IFEM_PATH}/src/LinAlg/Test/TestISTLPETScMatrix.C)
|
||||
endif()
|
||||
|
||||
IFEM_add_test_app("${TEST_SOURCES}"
|
||||
@@ -66,6 +72,9 @@ macro(IFEM_add_unittests IFEM_PATH)
|
||||
if(PETSC_FOUND)
|
||||
list(APPEND TEST_SRCS_MPI ${IFEM_PATH}/src/LinAlg/Test/MPI/TestPETScMatrix.C)
|
||||
endif()
|
||||
if(ISTL_FOUND)
|
||||
list(APPEND TEST_SRCS_MPI ${IFEM_PATH}/src/LinAlg/Test/MPI/TestISTLMatrix.C)
|
||||
endif()
|
||||
add_executable(IFEM-MPI-test EXCLUDE_FROM_ALL
|
||||
${IFEM_PATH}/src/IFEM-test.C ${TEST_SRCS_MPI})
|
||||
target_link_libraries(IFEM-MPI-test ${IFEM_LIBRARIES} ${IFEM_DEPLIBS} gtest)
|
||||
|
||||
@@ -52,6 +52,8 @@ bool AlgEqSystem::init (SystemMatrix::Type mtype, const LinSolParams* spar,
|
||||
SystemVector::Type vtype = SystemVector::STD;
|
||||
if (mtype == SystemMatrix::PETSC)
|
||||
vtype = SystemVector::PETSC;
|
||||
if (mtype == SystemMatrix::ISTL)
|
||||
vtype = SystemVector::ISTL;
|
||||
|
||||
for (i = 0; i < b.size(); i++)
|
||||
if (!b[i])
|
||||
|
||||
@@ -114,6 +114,13 @@ int IFEM::Init (int arg_c, char** arg_v, const char* title)
|
||||
std::cout <<"disabled";
|
||||
#endif
|
||||
|
||||
std::cout <<"\n ISTL support: ";
|
||||
#if HAS_ISTL
|
||||
std::cout <<"enabled (v"<< ISTL_VERSION <<")";
|
||||
#else
|
||||
std::cout <<"disabled";
|
||||
#endif
|
||||
|
||||
std::cout <<"\n VTF support: ";
|
||||
#if HAS_VTFAPI == 2
|
||||
std::cout <<"enabled (v2)";
|
||||
|
||||
253
src/LinAlg/ISTLMatrix.C
Normal file
253
src/LinAlg/ISTLMatrix.C
Normal file
@@ -0,0 +1,253 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file ISTLMatrix.C
|
||||
//!
|
||||
//! \date Mar 2 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Representation of the system matrix in ISTL format.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "ISTLMatrix.h"
|
||||
#include "LinSolParams.h"
|
||||
#include "LinAlgInit.h"
|
||||
#include "ProcessAdm.h"
|
||||
#include "SIMenums.h"
|
||||
#include "ASMstruct.h"
|
||||
#include "Profiler.h"
|
||||
#include "SAMpatch.h"
|
||||
#include "DomainDecomposition.h"
|
||||
|
||||
|
||||
ISTLVector::ISTLVector(const ProcessAdm& padm) : adm(padm)
|
||||
{
|
||||
LinAlgInit::increfs();
|
||||
}
|
||||
|
||||
|
||||
ISTLVector::ISTLVector(const ProcessAdm& padm, size_t n) : adm(padm)
|
||||
{
|
||||
x.resize(n);
|
||||
LinAlgInit::increfs();
|
||||
}
|
||||
|
||||
|
||||
ISTLVector::ISTLVector(const ProcessAdm& padm, const Real* values, size_t n) : adm(padm)
|
||||
{
|
||||
x.resize(n);
|
||||
this->restore(values);
|
||||
LinAlgInit::increfs();
|
||||
}
|
||||
|
||||
|
||||
ISTLVector::ISTLVector(const ISTLVector& vec) : adm(vec.adm)
|
||||
{
|
||||
x = vec.x;
|
||||
LinAlgInit::increfs();
|
||||
}
|
||||
|
||||
|
||||
ISTLVector::~ISTLVector()
|
||||
{
|
||||
LinAlgInit::decrefs();
|
||||
}
|
||||
|
||||
|
||||
void ISTLVector::init(Real value)
|
||||
{
|
||||
StdVector::init(value);
|
||||
x = value;
|
||||
}
|
||||
|
||||
size_t ISTLVector::dim() const
|
||||
{
|
||||
return x.size();
|
||||
}
|
||||
|
||||
|
||||
void ISTLVector::redim(size_t n)
|
||||
{
|
||||
x.resize(n);
|
||||
StdVector::redim(n);
|
||||
}
|
||||
|
||||
|
||||
bool ISTLVector::beginAssembly()
|
||||
{
|
||||
for (size_t i = 0; i < size(); ++i)
|
||||
x[i] = (*this)[i];
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool ISTLVector::endAssembly()
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
Real ISTLVector::L1norm() const
|
||||
{
|
||||
return x.one_norm();
|
||||
}
|
||||
|
||||
|
||||
Real ISTLVector::L2norm() const
|
||||
{
|
||||
return x.two_norm();
|
||||
}
|
||||
|
||||
|
||||
Real ISTLVector::Linfnorm() const
|
||||
{
|
||||
return x.infinity_norm();
|
||||
}
|
||||
|
||||
|
||||
ISTLMatrix::ISTLMatrix (const ProcessAdm& padm, const LinSolParams& spar,
|
||||
LinAlg::LinearSystemType ltype) :
|
||||
SparseMatrix(SUPERLU, 1), adm(padm),
|
||||
solParams(spar, adm), linsysType(ltype)
|
||||
{
|
||||
LinAlgInit::increfs();
|
||||
|
||||
setParams = true;
|
||||
nLinSolves = 0;
|
||||
}
|
||||
|
||||
|
||||
ISTLMatrix::ISTLMatrix (const ISTLMatrix& B) :
|
||||
adm(B.adm), solParams(B.solParams.get(), B.adm)
|
||||
{
|
||||
A = B.A;
|
||||
|
||||
LinAlgInit::increfs();
|
||||
|
||||
setParams = true;
|
||||
}
|
||||
|
||||
|
||||
ISTLMatrix::~ISTLMatrix ()
|
||||
{
|
||||
LinAlgInit::decrefs();
|
||||
}
|
||||
|
||||
|
||||
void ISTLMatrix::initAssembly (const SAM& sam, bool b)
|
||||
{
|
||||
SparseMatrix::initAssembly(sam, b);
|
||||
SparseMatrix::preAssemble(sam, b);
|
||||
|
||||
std::vector<std::set<int>> dofc;
|
||||
sam.getDofCouplings(dofc);
|
||||
|
||||
// Set correct number of rows and columns for matrix.
|
||||
size_t sum = 0;
|
||||
for (const auto& it : dofc)
|
||||
sum += it.size();
|
||||
|
||||
A.setSize(rows(), cols(), sum);
|
||||
A.setBuildMode(ISTL::Mat::random);
|
||||
|
||||
for (size_t i = 0; i < dofc.size(); ++i)
|
||||
A.setrowsize(i,dofc[i].size());
|
||||
A.endrowsizes();
|
||||
|
||||
for (size_t i = 0; i < dofc.size(); ++i)
|
||||
for (const auto& it : dofc[i])
|
||||
A.addindex(i, it-1);
|
||||
|
||||
A.endindices();
|
||||
|
||||
A = 0;
|
||||
}
|
||||
|
||||
bool ISTLMatrix::beginAssembly()
|
||||
{
|
||||
for (size_t j = 0; j < cols(); ++j)
|
||||
for (int i = IA[j]; i < IA[j+1]; ++i)
|
||||
A[JA[i]][j] = SparseMatrix::A[i];
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool ISTLMatrix::endAssembly()
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
void ISTLMatrix::init ()
|
||||
{
|
||||
SparseMatrix::init();
|
||||
|
||||
// Set all matrix elements to zero
|
||||
A = 0;
|
||||
}
|
||||
|
||||
|
||||
|
||||
bool ISTLMatrix::solve (SystemVector& B, bool newLHS, Real*)
|
||||
{
|
||||
if (!pre)
|
||||
std::tie(solver, pre, op) = solParams.setupPC(A);
|
||||
|
||||
ISTLVector* Bptr = dynamic_cast<ISTLVector*>(&B);
|
||||
if (!Bptr || !solver || !pre)
|
||||
return false;
|
||||
|
||||
try {
|
||||
Dune::InverseOperatorResult r;
|
||||
ISTL::Vec b(Bptr->getVector());
|
||||
Bptr->getVector() = 0;
|
||||
solver->apply(Bptr->getVector(), b, r);
|
||||
} catch (Dune::ISTLError& e) {
|
||||
std::cerr << "ISTL exception " << e << std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < rows(); ++i)
|
||||
(*Bptr)(i+1) = Bptr->getVector()[i];
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
bool ISTLMatrix::solve (const SystemVector& b, SystemVector& x, bool newLHS)
|
||||
{
|
||||
if (!pre)
|
||||
std::tie(solver, pre, op) = solParams.setupPC(A);
|
||||
|
||||
const ISTLVector* Bptr = dynamic_cast<const ISTLVector*>(&b);
|
||||
if (!Bptr || ! solver || !pre)
|
||||
return false;
|
||||
|
||||
ISTLVector* Xptr = dynamic_cast<ISTLVector*>(&x);
|
||||
if (!Xptr)
|
||||
return false;
|
||||
|
||||
try {
|
||||
Dune::InverseOperatorResult r;
|
||||
solver->apply(Xptr->getVector(),
|
||||
const_cast<ISTL::Vec&>(Bptr->getVector()), r);
|
||||
} catch (Dune::ISTLError& e) {
|
||||
std::cerr << "ISTL exception " << e << std::endl;
|
||||
return false;
|
||||
}
|
||||
|
||||
for (size_t i = 0; i < rows(); ++i)
|
||||
(*Xptr)(i+1) = Xptr->getVector()[i];
|
||||
|
||||
return true;
|
||||
}
|
||||
|
||||
|
||||
Real ISTLMatrix::Linfnorm () const
|
||||
{
|
||||
return A.infinity_norm();
|
||||
}
|
||||
163
src/LinAlg/ISTLMatrix.h
Normal file
163
src/LinAlg/ISTLMatrix.h
Normal file
@@ -0,0 +1,163 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file ISTLMatrix.h
|
||||
//!
|
||||
//! \date Mar 2 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Representation of the system matrix in ISTL format with interface
|
||||
//! to ISTL routines for solving linear equation systems.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#ifndef _ISTL_MATRIX_H
|
||||
#define _ISTL_MATRIX_H
|
||||
|
||||
#include "SystemMatrix.h"
|
||||
#include "SparseMatrix.h"
|
||||
#include "ISTLSolParams.h"
|
||||
#include "LinAlgenums.h"
|
||||
#include <memory>
|
||||
|
||||
class LinSolParams;
|
||||
class SAMpatch;
|
||||
|
||||
|
||||
/*!
|
||||
\brief Class for representing the system vector in ISTL format.
|
||||
*/
|
||||
|
||||
class ISTLVector : public StdVector
|
||||
{
|
||||
public:
|
||||
//! \brief Constructor creating an empty vector.
|
||||
ISTLVector(const ProcessAdm& padm);
|
||||
//! \brief Constructor creating a vector of length \a n.
|
||||
ISTLVector(const ProcessAdm& padm, size_t n);
|
||||
//! \brief Constructor creating a vector from an array.
|
||||
ISTLVector(const ProcessAdm& padm, const Real* values, size_t n);
|
||||
//! \brief Copy constructor.
|
||||
ISTLVector(const ISTLVector& vec);
|
||||
//! \brief Destructor.
|
||||
virtual ~ISTLVector();
|
||||
#endif
|
||||
|
||||
//! \brief Returns the vector type.
|
||||
virtual Type getType() const { return ISTL; }
|
||||
|
||||
//! \brief Initializes the vector to a given scalar value.
|
||||
virtual void init(Real value = Real(0));
|
||||
|
||||
//! \brief Returns the dimension of the system vector.
|
||||
virtual size_t dim() const;
|
||||
|
||||
//! \brief Sets the dimension of the system vector.
|
||||
virtual void redim(size_t n);
|
||||
|
||||
//! \brief Begins communication step needed in parallel vector assembly.
|
||||
//! \details Must be called together with endAssembly after vector assembly
|
||||
//! is completed on each processor and before the linear system is solved.
|
||||
virtual bool beginAssembly();
|
||||
|
||||
//! \brief Ends communication step needed in parallel vector assembly.
|
||||
//! \details Must be called together with beginAssembly after vector assembly
|
||||
//! is completed on each processor and before the linear system is solved.
|
||||
virtual bool endAssembly();
|
||||
|
||||
//! \brief L1-norm of vector.
|
||||
virtual Real L1norm() const;
|
||||
|
||||
//! \brief L2-norm of vector.
|
||||
virtual Real L2norm() const;
|
||||
|
||||
//! \brief Linfinity-norm of vector.
|
||||
virtual Real Linfnorm() const;
|
||||
|
||||
//! \brief Return associated process administrator
|
||||
const ProcessAdm& getAdm() const { return adm; }
|
||||
typedef Dune::BlockVector<Dune::FieldVector<double,1>> Vec;
|
||||
|
||||
Vec& getVector() { return x; }
|
||||
const Vec& getVector() const { return x; }
|
||||
|
||||
protected:
|
||||
ISTL::Vec x; //!< ISTL vector
|
||||
const ProcessAdm& adm; //!< Process administrator
|
||||
};
|
||||
|
||||
|
||||
/*!
|
||||
\brief Class for representing the system matrix in ISTL format.
|
||||
*/
|
||||
|
||||
class ISTLMatrix : public SparseMatrix
|
||||
{
|
||||
public:
|
||||
//! \brief Constructor.
|
||||
ISTLMatrix(const ProcessAdm& padm, const LinSolParams& spar,
|
||||
LinAlg::LinearSystemType ltype);
|
||||
//! \brief Copy constructor.
|
||||
ISTLMatrix(const ISTLMatrix& A);
|
||||
//! \brief The destructor frees the dynamically allocated arrays.
|
||||
virtual ~ISTLMatrix();
|
||||
|
||||
//! \brief Returns the matrix type.
|
||||
virtual Type getType() const { return ISTL; }
|
||||
|
||||
//! \brief Returns the dimension of the system matrix.
|
||||
virtual size_t dim(int = 1) const { return 0; }
|
||||
|
||||
//! \brief Creates a copy of the system matrix and returns a pointer to it.
|
||||
virtual SystemMatrix* copy() const { return new ISTLMatrix(*this); }
|
||||
|
||||
//! \brief Initializes the element assembly process.
|
||||
//! \details Must be called once before the element assembly loop.
|
||||
//! The PETSc data structures are initialized and the all symbolic operations
|
||||
//! that are needed before the actual assembly can start are performed.
|
||||
//! \param[in] sam Auxiliary data describing the FE model topology, etc.
|
||||
virtual void initAssembly(const SAM& sam, bool);
|
||||
|
||||
//! \brief Initializes the matrix to zero assuming it is properly dimensioned.
|
||||
virtual void init();
|
||||
|
||||
//! \brief Begins communication step needed in parallel matrix assembly.
|
||||
//! \details Must be called together with endAssembly after matrix assembly
|
||||
//! is completed on each processor and before the linear system is solved.
|
||||
virtual bool beginAssembly();
|
||||
//! \brief Ends communication step needed in parallel matrix assembly.
|
||||
//! \details Must be called together with beginAssembly after matrix assembly
|
||||
//! is completed on each processor and before the linear system is solved.
|
||||
virtual bool endAssembly();
|
||||
|
||||
//! \brief Solves the linear system of equations for a given right-hand-side.
|
||||
//! \param B Right-hand-side vector on input, solution vector on output
|
||||
//! \param[in] newLHS \e true if the left-hand-side matrix has been updated
|
||||
virtual bool solve(SystemVector& B, bool newLHS, Real*);
|
||||
|
||||
//! \brief Solves the linear system of equations for a given right-hand-side.
|
||||
//! \param[in] B Right-hand-side vector
|
||||
//! \param[out] x Solution vector
|
||||
//! \param[in] newLHS \e true if the left-hand-side matrix has been updated
|
||||
virtual bool solve(const SystemVector& B, SystemVector& x, bool newLHS);
|
||||
|
||||
//! \brief Returns the L-infinity norm of the matrix.
|
||||
virtual Real Linfnorm() const;
|
||||
|
||||
//! \brief Returns the ISTL matrix (for assignment).
|
||||
virtual ISTL::Mat& getMatrix() { return A; }
|
||||
//! \brief Returns the ISTL matrix (for read access).
|
||||
virtual const ISTL::Mat& getMatrix() const { return A; }
|
||||
|
||||
protected:
|
||||
ISTL::Mat A; //!< The actual ISTL matrix
|
||||
std::unique_ptr<ISTL::Operator> op; //!< The matrix adapter
|
||||
std::unique_ptr<ISTL::InverseOperator> solver; //!< Solver to use
|
||||
std::unique_ptr<ISTL::Preconditioner> pre; //!< Preconditioner to use
|
||||
const ProcessAdm& adm; //!< Process administrator
|
||||
ISTLSolParams solParams; //!< Linear solver parameters
|
||||
bool setParams; //!< If linear solver parameters are set
|
||||
int nLinSolves; //!< Number of linear solves
|
||||
LinAlg::LinearSystemType linsysType; //!< Linear system type
|
||||
};
|
||||
510
src/LinAlg/ISTLSolParams.C
Normal file
510
src/LinAlg/ISTLSolParams.C
Normal file
@@ -0,0 +1,510 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file ISTLSolParams.C
|
||||
//!
|
||||
//! \date Mar 10 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Linear solver parameters for ISTL.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "ISTLSolParams.h"
|
||||
#include "ASMstruct.h"
|
||||
#include "LinSolParams.h"
|
||||
#include "ProcessAdm.h"
|
||||
#include "SAMpatch.h"
|
||||
#include <dune/istl/overlappingschwarz.hh>
|
||||
#include <dune/istl/paamg/amg.hh>
|
||||
#include <dune/istl/matrixmatrix.hh>
|
||||
#include <dune/common/version.hh>
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
|
||||
#include <dune/istl/paamg/fastamg.hh>
|
||||
#include <dune/istl/paamg/twolevelmethod.hh>
|
||||
#endif
|
||||
|
||||
|
||||
namespace ISTL {
|
||||
|
||||
BlockPreconditioner::BlockPreconditioner(const ISTL::Mat& A,
|
||||
const DomainDecomposition& dd_,
|
||||
const std::string& schurType) :
|
||||
dd(dd_)
|
||||
{
|
||||
blocks.resize(dd.getNoBlocks()*dd.getNoBlocks()+1);
|
||||
size_t k = 0;
|
||||
for (size_t i = 0; i < dd.getNoBlocks(); ++i)
|
||||
for (size_t j = 0; j < dd.getNoBlocks(); ++j)
|
||||
extractBlock(blocks[k++], A, dd.getBlockEqs(i), dd.getBlockEqs(j), dd.getG2LEQ(j));
|
||||
|
||||
// Scale blocks[1] with inverse diagonal of blocks[0]
|
||||
blocks.back() = blocks[1];
|
||||
for (auto row = blocks.back().begin(); row != blocks.back().end(); ++row)
|
||||
for (auto it = row->begin(); it != row->end(); ++it)
|
||||
*it /= blocks[0][row.index()][row.index()];
|
||||
|
||||
// blocks[1] = D - C*diag(A)^-1*B
|
||||
ISTL::Mat SP;
|
||||
Dune::matMultMat(SP, blocks[2], blocks.back());
|
||||
blocks[2] = blocks[1];
|
||||
subtractMatrices(blocks[1], blocks[3], SP);
|
||||
|
||||
// clear memory for unused blocks
|
||||
blocks.resize(schurType == "diag" ? 2 : 3);
|
||||
|
||||
blockPre.resize(dd.getNoBlocks());
|
||||
blockOp.resize(dd.getNoBlocks());
|
||||
}
|
||||
|
||||
|
||||
void BlockPreconditioner::pre(ISTL::Vec& x, ISTL::Vec& b)
|
||||
{
|
||||
ISTL::Vec tempx, tempb;
|
||||
for (size_t block = 0; block < dd.getNoBlocks(); ++block) {
|
||||
ISTL::Vec tempx, tempb;
|
||||
tempx.resize(dd.getBlockEqs(block).size());
|
||||
tempb.resize(dd.getBlockEqs(block).size());
|
||||
size_t i = 0;
|
||||
for (auto& it : dd.getBlockEqs(block)) {
|
||||
tempx[i] = x[it-1];
|
||||
tempb[i++] = b[it-1];
|
||||
}
|
||||
|
||||
blockPre[block]->pre(tempx, tempb);
|
||||
|
||||
i = 0;
|
||||
for (auto& it : dd.getBlockEqs(block)) {
|
||||
x[it-1] = tempx[i];
|
||||
b[it-1] = tempb[i++];
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void BlockPreconditioner::apply(ISTL::Vec& v, const ISTL::Vec& d)
|
||||
{
|
||||
// backwards due to upper schur complement
|
||||
ISTL::Vec tempx, tempb;
|
||||
for (int block = dd.getNoBlocks()-1; block >= 0; --block) {
|
||||
tempb.resize(dd.getBlockEqs(block).size());
|
||||
size_t i = 0;
|
||||
for (auto& it : dd.getBlockEqs(block))
|
||||
tempb[i++] = d[it-1];
|
||||
|
||||
if (block == 0 && blocks.size() > 2)
|
||||
blocks.back().mmv(tempx, tempb);
|
||||
|
||||
tempx.resize(dd.getBlockEqs(block).size());
|
||||
tempx = 0;
|
||||
blockPre[block]->apply(tempx, tempb);
|
||||
|
||||
i = 0;
|
||||
for (auto& it : dd.getBlockEqs(block))
|
||||
v[it-1] = tempx[i++];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void BlockPreconditioner::post(ISTL::Vec& x)
|
||||
{
|
||||
ISTL::Vec tempx, tempb;
|
||||
for (size_t block = 0; block < dd.getNoBlocks(); ++block) {
|
||||
ISTL::Vec tempx;
|
||||
tempx.resize(dd.getBlockEqs(block).size());
|
||||
size_t i = 0;
|
||||
for (auto& it : dd.getBlockEqs(block))
|
||||
tempx[i++] = x[it-1];
|
||||
|
||||
blockPre[block]->post(tempx);
|
||||
|
||||
i = 0;
|
||||
for (auto& it : dd.getBlockEqs(block))
|
||||
x[it-1] = tempx[i++];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void BlockPreconditioner::extractBlock(ISTL::Mat& B, const ISTL::Mat& A,
|
||||
const std::set<int>& eqs_row1,
|
||||
const std::set<int>& eqs_col1,
|
||||
const std::map<int,int>& eqs_col_g2l)
|
||||
{
|
||||
std::vector<int> eqs_row(eqs_row1.begin(), eqs_row1.end());
|
||||
std::vector<int> eqs_col(eqs_col1.begin(), eqs_col1.end());
|
||||
size_t sum=0;
|
||||
std::vector<std::set<int>> adj;
|
||||
adj.resize(eqs_row.size());
|
||||
size_t i = 0;
|
||||
|
||||
for (auto& it3 : eqs_row) {
|
||||
auto it = A.begin()+it3-1;
|
||||
for (auto it2 = it->begin(); it2 != it->end(); ++it2) {
|
||||
auto itc = eqs_col_g2l.find(it2.index()+1);
|
||||
if (itc != eqs_col_g2l.end()) {
|
||||
adj[i].insert(itc->second-1);
|
||||
++sum;
|
||||
}
|
||||
}
|
||||
++i;
|
||||
}
|
||||
|
||||
B.setSize(eqs_row.size(), eqs_col.size(), sum);
|
||||
B.setBuildMode(ISTL::Mat::random);
|
||||
|
||||
for (size_t i = 0; i < adj.size(); i++)
|
||||
B.setrowsize(i,adj[i].size());
|
||||
B.endrowsizes();
|
||||
|
||||
for (size_t i = 0; i < adj.size(); i++)
|
||||
for (auto& it : adj[i])
|
||||
B.addindex(i, it);
|
||||
|
||||
B.endindices();
|
||||
B = 0;
|
||||
for (auto it = B.begin(); it != B.end(); ++it)
|
||||
for (auto it2 = it->begin(); it2 != it->end(); ++it2)
|
||||
*it2 = A[eqs_row[it.index()]-1][eqs_col[it2.index()]-1];
|
||||
}
|
||||
|
||||
|
||||
void BlockPreconditioner::subtractMatrices(ISTL::Mat& A, const ISTL::Mat& B,
|
||||
const ISTL::Mat& C)
|
||||
{
|
||||
size_t sum=0;
|
||||
std::vector<std::set<int>> adj;
|
||||
adj.resize(B.N());
|
||||
for (auto row = B.begin(); row != B.end(); ++row) {
|
||||
for (auto it = row->begin(); it != row->end(); ++it)
|
||||
adj[row.index()].insert(it.index());
|
||||
}
|
||||
for (auto row = C.begin(); row != C.end(); ++row) {
|
||||
for (auto it = row->begin(); it != row->end(); ++it)
|
||||
adj[row.index()].insert(it.index());
|
||||
sum += adj[row.index()].size();
|
||||
}
|
||||
|
||||
A.setSize(B.N(), B.M(), sum);
|
||||
A.setBuildMode(ISTL::Mat::random);
|
||||
|
||||
for (size_t i = 0; i < adj.size(); i++)
|
||||
A.setrowsize(i,adj[i].size());
|
||||
A.endrowsizes();
|
||||
|
||||
for (size_t i = 0; i < adj.size(); i++)
|
||||
for (auto& it : adj[i])
|
||||
A.addindex(i, it);
|
||||
|
||||
A.endindices();
|
||||
A = 0;
|
||||
for (auto row = B.begin(); row != B.end(); ++row)
|
||||
for (auto it = row->begin(); it != row->end(); ++it)
|
||||
A[row.index()][it.index()] = *it;
|
||||
|
||||
for (auto row = C.begin(); row != C.end(); ++row)
|
||||
for (auto it = row->begin(); it != row->end(); ++it)
|
||||
A[row.index()][it.index()] -= *it;
|
||||
}
|
||||
|
||||
|
||||
} // namespace ISTL
|
||||
|
||||
|
||||
/*! \brief Helper template for setting up a solver with the appropriate
|
||||
preconditioner type.
|
||||
\details We cannot instance using dynamic polymorphism, the solver need
|
||||
access to the real type for the preconditioner. We can however
|
||||
call the solver in the interface class scope afterwards.
|
||||
*/
|
||||
template<class Prec>
|
||||
static ISTL::InverseOperator* setupWithPreType(const LinSolParams& solParams,
|
||||
ISTL::Operator& op,
|
||||
Prec& pre)
|
||||
{
|
||||
std::string type = solParams.getStringValue("type");
|
||||
double rtol = solParams.getDoubleValue("rtol");
|
||||
int maxits = solParams.getIntValue("maxits");
|
||||
int verbosity = solParams.getIntValue("verbosity");
|
||||
int restart = solParams.getIntValue("gmres_restart_iterations");
|
||||
if (type == "bcgs")
|
||||
return new Dune::BiCGSTABSolver<ISTL::Vec>(op, pre, rtol, maxits, verbosity);
|
||||
else if (type == "cg")
|
||||
return new Dune::CGSolver<ISTL::Vec>(op, pre, rtol, maxits, verbosity);
|
||||
else if (type == "minres")
|
||||
return new Dune::MINRESSolver<ISTL::Vec>(op, pre, rtol, maxits, verbosity);
|
||||
else if (type == "gmres")
|
||||
return new Dune::RestartedGMResSolver<ISTL::Vec>(op, pre, rtol, restart,
|
||||
maxits, verbosity);
|
||||
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
|
||||
/*! \brief Helper template for setting up an AMG preconditioner with a given smoother */
|
||||
|
||||
template<class Smoother>
|
||||
static ISTL::Preconditioner* setupAMG(const LinSolParams& params,
|
||||
size_t block, ISTL::Operator& op,
|
||||
std::unique_ptr<ISTL::InverseOperator>* solver)
|
||||
{
|
||||
// The coupling metric used in the AMG
|
||||
typedef Dune::Amg::FirstDiagonal CouplingMetric;
|
||||
// The coupling criterion used in the AMG
|
||||
typedef Dune::Amg::SymmetricCriterion<ISTL::Mat, CouplingMetric> CritBase;
|
||||
// The coarsening criterion used in the AMG
|
||||
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
|
||||
Criterion crit;
|
||||
typedef typename Dune::Amg::AMG<ISTL::Operator, ISTL::Vec, Smoother> AMG;
|
||||
typename AMG::SmootherArgs args;
|
||||
args.relaxationFactor = 1.0;
|
||||
args.iterations = std::max(1, params.getBlock(block).getIntValue("multigrid_no_smooth"));
|
||||
|
||||
if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
|
||||
crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
|
||||
|
||||
if (params.getBlock(block).hasValue("multigrid_no_smooth")) {
|
||||
int val = params.getBlock(block).getIntValue("multigrid_no_smooth");
|
||||
crit.setNoPreSmoothSteps(val);
|
||||
crit.setNoPostSmoothSteps(val);
|
||||
}
|
||||
|
||||
if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
|
||||
crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
|
||||
|
||||
auto result = new AMG(op, crit, args);
|
||||
if (solver)
|
||||
solver->reset(setupWithPreType(params, op, *result));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
|
||||
/*! \brief Helper template for setting up a AMG preconditioner
|
||||
*! with fine smoother differing from the other levels */
|
||||
|
||||
template<class Smoother, class FineSmoother>
|
||||
static ISTL::Preconditioner* setupAMG2_full(const LinSolParams& params, size_t block,
|
||||
ISTL::Operator& op,
|
||||
std::unique_ptr<ISTL::InverseOperator>* solver,
|
||||
FineSmoother* fsmooth)
|
||||
{
|
||||
typedef Dune::Amg::FirstDiagonal CouplingMetric;
|
||||
typedef Dune::Amg::SymmetricCriterion<ISTL::Mat,CouplingMetric> CritBase;
|
||||
typedef Dune::Amg::CoarsenCriterion<CritBase> Criterion;
|
||||
typedef Dune::Amg::AggregationLevelTransferPolicy<ISTL::Operator,Criterion> TransferPolicy;
|
||||
typedef Dune::Amg::OneStepAMGCoarseSolverPolicy<ISTL::Operator,Smoother,Criterion> CoarsePolicy;
|
||||
typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
|
||||
typedef typename Dune::Amg::TwoLevelMethod<ISTL::Operator, CoarsePolicy, FineSmoother> AMG2;
|
||||
|
||||
SmootherArgs args;
|
||||
args.relaxationFactor = 1.0;
|
||||
|
||||
Criterion crit;
|
||||
if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
|
||||
crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
|
||||
|
||||
if (params.getBlock(block).hasValue("multigrid_no_smooth")) {
|
||||
int val = params.getBlock(block).getIntValue("multigrid_no_smooth");
|
||||
crit.setNoPreSmoothSteps(val);
|
||||
crit.setNoPostSmoothSteps(val);
|
||||
}
|
||||
|
||||
if (params.getBlock(block).hasValue("multigrid_max_coarse_size"))
|
||||
crit.setCoarsenTarget(params.getBlock(block).getIntValue("multigrid_max_coarse_size"));
|
||||
|
||||
CoarsePolicy coarsePolicy(args, crit);
|
||||
TransferPolicy policy(crit);
|
||||
Dune::shared_ptr<FineSmoother> fsp(fsmooth);
|
||||
auto result = new AMG2(op, fsp, policy, coarsePolicy);
|
||||
if (solver)
|
||||
solver->reset(setupWithPreType(params, op, *result));
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
template<class FineSmoother>
|
||||
static ISTL::Preconditioner* setupAMG2_smoother(const LinSolParams& params, size_t block,
|
||||
ISTL::Operator& op,
|
||||
std::unique_ptr<ISTL::InverseOperator>* solver,
|
||||
FineSmoother* fsmooth)
|
||||
{
|
||||
std::string smoother = params.getBlock(block).getStringValue("multigrid_smoother");
|
||||
if (smoother == "ilu")
|
||||
return setupAMG2_full<ISTL::ILU0>(params, block, op, solver, fsmooth);
|
||||
else if (smoother == "sor")
|
||||
return setupAMG2_full<ISTL::SOR>(params, block, op, solver, fsmooth);
|
||||
else if (smoother == "ssor")
|
||||
return setupAMG2_full<ISTL::SSOR>(params, block, op, solver, fsmooth);
|
||||
else if (smoother == "jacobi")
|
||||
return setupAMG2_full<ISTL::GJ>(params, block, op, solver, fsmooth);
|
||||
else {
|
||||
std::cerr << "** ISTLSolParams ** Invalid smoother " << smoother << "." << std::endl;
|
||||
return nullptr;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
static ISTL::Preconditioner* setupAMG2(const LinSolParams& params, size_t block,
|
||||
ISTL::Operator& op, std::unique_ptr<ISTL::InverseOperator>* solver)
|
||||
{
|
||||
std::string smoother = params.getBlock(block).getStringValue("multigrid_finesmoother");
|
||||
if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "ilu") {
|
||||
auto fsmooth = new ISTL::ILU0(op.getmat(), 1.0);
|
||||
return setupAMG2_smoother(params, block, op, solver, fsmooth);
|
||||
} else if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "sor") {
|
||||
auto fsmooth = new ISTL::SOR(op.getmat(), 1, 1.0);
|
||||
return setupAMG2_smoother(params, block, op, solver, fsmooth);
|
||||
} else if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "ssor") {
|
||||
auto fsmooth = new ISTL::SSOR(op.getmat(), 1, 1.0);
|
||||
return setupAMG2_smoother(params, block, op, solver, fsmooth);
|
||||
} else if (params.getBlock(block).getStringValue("multigrid_finesmoother") == "jacobi") {
|
||||
auto fsmooth = new ISTL::GJ(op.getmat(), 1, 1.0);
|
||||
return setupAMG2_smoother(params, block, op, solver, fsmooth);
|
||||
} else {
|
||||
std::cerr << "** ISTLSolParams ** Invalid fine smoother " << smoother << "." << std::endl;
|
||||
return nullptr;
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
|
||||
/*! \brief Conditionally setup a KSP */
|
||||
template<class Type>
|
||||
static ISTL::Preconditioner*
|
||||
setupSolver(Type* pre,
|
||||
ISTL::Operator& op,
|
||||
const LinSolParams& solParams,
|
||||
std::unique_ptr<ISTL::InverseOperator>* solver)
|
||||
{
|
||||
if (solver)
|
||||
solver->reset(setupWithPreType(solParams, op, *pre));
|
||||
|
||||
return pre;
|
||||
}
|
||||
|
||||
ISTL::Preconditioner* ISTLSolParams::setupPCInternal(ISTL::Mat& A,
|
||||
ISTL::Operator& op,
|
||||
size_t block,
|
||||
std::unique_ptr<ISTL::InverseOperator>* solver)
|
||||
{
|
||||
std::string prec = solParams.getBlock(block).getStringValue("pc");
|
||||
if (prec == "ilu") {
|
||||
int fill_level = solParams.getBlock(block).getIntValue("ilu_fill_level");
|
||||
if (fill_level == 0)
|
||||
return setupSolver(new ISTL::ILU0(A, 1.0), op, solParams, solver);
|
||||
else
|
||||
return setupSolver(new ISTL::ILUn(A, fill_level, 1.0),
|
||||
op, solParams, solver);
|
||||
} else if (prec == "lu")
|
||||
return setupSolver(new ISTL::LU(new ISTL::LUType(A)),
|
||||
op, solParams, solver);
|
||||
else if (prec == "sor")
|
||||
return setupSolver(new ISTL::SOR(A, 1, 1.0), op, solParams, solver);
|
||||
else if (prec == "ssor")
|
||||
return setupSolver(new ISTL::SSOR(A, 1, 1.0), op, solParams, solver);
|
||||
else if (prec == "jacobi")
|
||||
return setupSolver(new ISTL::GJ(A, 1, 1.0), op, solParams, solver);
|
||||
else if (prec == "gs")
|
||||
return setupSolver(new ISTL::GS(A, 1, 1.0), op, solParams, solver);
|
||||
else if (prec == "asm" || prec == "asmlu") {
|
||||
size_t nx = solParams.getBlock(block).getIntValue("asm_nx");
|
||||
nx = std::max(1ul, nx);
|
||||
size_t ny = solParams.getBlock(block).getIntValue("asm_ny");
|
||||
ny = std::max(1ul, ny);
|
||||
size_t nz = solParams.getBlock(block).getIntValue("asm_nz");
|
||||
nz = std::max(1ul, nz);
|
||||
int overlap = solParams.getBlock(block).getIntValue("asm_overlap");
|
||||
auto locSubdDofs = adm.dd.getSubdomains(nx, ny, nz, overlap, block);
|
||||
|
||||
if (prec == "asmlu") {
|
||||
ISTL::ASMLU::subdomain_vector ddofs(locSubdDofs.size());
|
||||
for (size_t i = 0; i < locSubdDofs.size(); ++i)
|
||||
for (const auto& it : locSubdDofs[i])
|
||||
ddofs[i].insert(it-1);
|
||||
|
||||
return setupSolver(new ISTL::ASMLU(A, ddofs), op, solParams, solver);
|
||||
} else {
|
||||
ISTL::ASM::subdomain_vector ddofs(locSubdDofs.size());
|
||||
for (size_t i = 0; i < locSubdDofs.size(); ++i) {
|
||||
for (size_t i = 0; i < locSubdDofs.size(); ++i)
|
||||
for (const auto& it : locSubdDofs[i])
|
||||
ddofs[i].insert(it-1);
|
||||
}
|
||||
|
||||
return setupSolver(new ISTL::ASM(A, ddofs), op, solParams, solver);
|
||||
}
|
||||
} else if (prec == "amg" || prec == "gamg") {
|
||||
if (solParams.getBlock(block).getStringValue("multigrid_smoother") !=
|
||||
solParams.getBlock(block).getStringValue("multigrid_finesmoother")) {
|
||||
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
|
||||
return setupAMG2(solParams, block, op, solver);
|
||||
#else
|
||||
std::cerr << "** ISTLSolParams ** Separate fine smoother not implemented." << std::endl;
|
||||
return nullptr;
|
||||
#endif
|
||||
} else {
|
||||
std::string smoother = solParams.getBlock(block).getStringValue("multigrid_smoother");
|
||||
if (smoother.empty()) {
|
||||
std::cerr << "** ISTLSolParams ** No smoother defined for AMG, defaulting to ILU" << std::endl;
|
||||
smoother = "ilu";
|
||||
}
|
||||
if (smoother == "ilu")
|
||||
return setupAMG<ISTL::ILU0>(solParams, block, op, solver);
|
||||
else if (smoother == "sor")
|
||||
return setupAMG<ISTL::SOR>(solParams, block, op, solver);
|
||||
else if (smoother == "ssor")
|
||||
return setupAMG<ISTL::SSOR>(solParams, block, op, solver);
|
||||
else if (smoother == "jacobi")
|
||||
return setupAMG<ISTL::GJ>(solParams, block, op, solver);
|
||||
}
|
||||
}
|
||||
|
||||
return nullptr;
|
||||
}
|
||||
|
||||
|
||||
std::tuple<std::unique_ptr<ISTL::InverseOperator>,
|
||||
std::unique_ptr<ISTL::Preconditioner>,
|
||||
std::unique_ptr<ISTL::Operator>> ISTLSolParams::setupPC(ISTL::Mat& A)
|
||||
{
|
||||
std::unique_ptr<ISTL::InverseOperator> solver;
|
||||
std::unique_ptr<ISTL::Preconditioner> pre;
|
||||
std::unique_ptr<ISTL::Operator> op;
|
||||
|
||||
if (solParams.getNoBlocks() > 2) {
|
||||
std::cerr << "*** ISTLSolParams ** More than two blocks are not implemented." << std::endl;
|
||||
return std::make_tuple(nullptr, nullptr, nullptr);
|
||||
}
|
||||
|
||||
if (solParams.getNoBlocks() > 1) {
|
||||
op.reset(new ISTL::Operator(A));
|
||||
ISTL::BlockPreconditioner* bpre = new ISTL::BlockPreconditioner(A, adm.dd, solParams.getStringValue("schur"));
|
||||
pre.reset(bpre);
|
||||
for (size_t i = 0; i < adm.dd.getNoBlocks(); ++i)
|
||||
bpre->getBlockPre(i).reset(setupPCInternal(bpre->getBlock(i),
|
||||
bpre->getBlockOp(i), i, nullptr));
|
||||
solver.reset(setupWithPreType(solParams, *op, *bpre));
|
||||
} else {
|
||||
#ifdef HAVE_MPI
|
||||
if (adm.isParallel()) {
|
||||
Dune::IndexInfoFromGrid<int, int> index;
|
||||
for (size_t i = 0; i < adm.dd.getMLGEQ().size(); ++i)
|
||||
index.addLocalIndex(std::make_tuple(adm.dd.getGlobalEq(i+1), i, 0));
|
||||
Dune::OwnerOverlapCopyCommunication<int,int> comm(index, *adm.getCommunicator());
|
||||
//ISTL::ParMatrixAdapter* nop = new ISTL::ParMatrixAdapter(A, comm);
|
||||
//op.reset(nop);
|
||||
//pre.reset(setupPCInternal(A, *nop, 0, &solver));
|
||||
} else
|
||||
#endif
|
||||
{
|
||||
ISTL::Operator* nop = new ISTL::Operator(A);
|
||||
op.reset(nop);
|
||||
pre.reset(setupPCInternal(A, *nop, 0, &solver));
|
||||
}
|
||||
}
|
||||
|
||||
return std::make_tuple(std::move(solver), std::move(pre), std::move(op));
|
||||
}
|
||||
258
src/LinAlg/ISTLSolParams.h
Normal file
258
src/LinAlg/ISTLSolParams.h
Normal file
@@ -0,0 +1,258 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file ISTLSolParams.h
|
||||
//!
|
||||
//! \date Mar 10 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Linear solver parameters for ISTL matrices.
|
||||
//! \details Includes linear solver method, preconditioner
|
||||
//! and convergence criteria.
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#ifndef _ISTL_SOLPARAMS_H
|
||||
#define _ISTL_SOLPARAMS_H
|
||||
|
||||
#include "ISTLSupport.h"
|
||||
#include <dune/istl/operators.hh>
|
||||
#include <dune/istl/solvercategory.hh>
|
||||
|
||||
#include <iostream>
|
||||
#include <set>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
|
||||
class DomainDecomposition;
|
||||
class LinSolParams;
|
||||
class ProcessAdm;
|
||||
|
||||
|
||||
/*! This implements a Schur-decomposition based preconditioner for the
|
||||
* block system
|
||||
* [A B]
|
||||
* [C D]
|
||||
*
|
||||
* The preconditioner is
|
||||
* [Apre B]
|
||||
* [ P]
|
||||
* Here Apre is some preconditioner for A and P some preconditioner for
|
||||
* S = D - C*diag(A)^-1*B. The B block may be dropped.
|
||||
!*/
|
||||
|
||||
namespace ISTL {
|
||||
|
||||
class BlockPreconditioner : public Preconditioner {
|
||||
public:
|
||||
// define the category
|
||||
enum {
|
||||
//! \brief The category the preconditioner is part of.
|
||||
category=Dune::SolverCategory::sequential
|
||||
};
|
||||
|
||||
//! \brief Constructor
|
||||
//! \param[in] A The system matrix
|
||||
//! \param[in] dd Domain decomposition
|
||||
BlockPreconditioner(const ISTL::Mat& A, const DomainDecomposition& dd_,
|
||||
const std::string& schurType);
|
||||
|
||||
//! \brief Destructor
|
||||
virtual ~BlockPreconditioner()
|
||||
{}
|
||||
|
||||
//! \brief Preprocess preconditioner
|
||||
virtual void pre(ISTL::Vec& x, ISTL::Vec& b);
|
||||
|
||||
//! \brief Applies the preconditioner
|
||||
//! \param[out] v The resulting vector
|
||||
//! \param[in] d The vector to apply the preconditioner to
|
||||
virtual void apply(ISTL::Vec& v, const ISTL::Vec& d);
|
||||
|
||||
//! \brief Post-process function
|
||||
virtual void post(ISTL::Vec& x);
|
||||
|
||||
//! \brief Obtain reference to a block preconditioner
|
||||
std::unique_ptr<ISTL::Preconditioner>& getBlockPre(size_t block)
|
||||
{
|
||||
return blockPre[block];
|
||||
}
|
||||
|
||||
//! \brief Obtain reference to a block matrix
|
||||
ISTL::Mat& getBlock(size_t block)
|
||||
{
|
||||
return blocks[block];
|
||||
}
|
||||
|
||||
//! \brief Obtain matrix adaptor for a block matrix
|
||||
ISTL::Operator& getBlockOp(size_t block)
|
||||
{
|
||||
if (!blockOp[block])
|
||||
blockOp[block].reset(new ISTL::Operator(blocks[block]));
|
||||
|
||||
return *blockOp[block];
|
||||
}
|
||||
|
||||
protected:
|
||||
//! \brief Build block from block equations.
|
||||
static void extractBlock(ISTL::Mat& B, const ISTL::Mat& A,
|
||||
const std::set<int>& eqs_row,
|
||||
const std::set<int>& eqs_col,
|
||||
const std::map<int,int>& eqs_col_g2l);
|
||||
|
||||
//! \brief Find A = B - C
|
||||
static void subtractMatrices(ISTL::Mat& A,
|
||||
const ISTL::Mat& B,
|
||||
const ISTL::Mat& C);
|
||||
|
||||
std::vector<std::unique_ptr<ISTL::Preconditioner>> blockPre; //!< The preconditioners
|
||||
std::vector<std::unique_ptr<ISTL::Operator>> blockOp; //!< The matrix adaptors for the blocks
|
||||
std::vector<ISTL::Mat> blocks; //!< Matrix blocks
|
||||
const DomainDecomposition& dd; //!< Domain decomposition
|
||||
};
|
||||
|
||||
#ifdef HAVE_MPI
|
||||
/**
|
||||
* \brief An overlapping schwarz operator.
|
||||
*
|
||||
* This operator represents a parallel matrix product using
|
||||
* sequential data structures together with a parallel index set
|
||||
* describing an overlapping domain decomposition and the communication.
|
||||
* \tparam M The type of the sequential matrix to use,
|
||||
* e.g. BCRSMatrix or another matrix type fulfilling the
|
||||
* matrix interface of ISTL.
|
||||
* \tparam X The type of the sequential vector to use for the left hand side,
|
||||
* e.g. BlockVector or another type fulfilling the ISTL
|
||||
* vector interface.
|
||||
* \tparam Y The type of the sequential vector to use for the right hand side,
|
||||
* e..g. BlockVector or another type fulfilling the ISTL
|
||||
* vector interface.
|
||||
* \tparam C The type of the communication object.
|
||||
* This must either be OwnerOverlapCopyCommunication or a type
|
||||
* implementing the same interface.
|
||||
*
|
||||
* This is a modified version which does summation of the vector after evaluation.
|
||||
*/
|
||||
template<class M, class X, class Y, class C>
|
||||
class OverlappingSchwarzOperator : public Dune::AssembledLinearOperator<M,X,Y>
|
||||
{
|
||||
public:
|
||||
//! \brief The type of the matrix we operate on.
|
||||
//!
|
||||
//! E.g. BCRSMatrix or another matrix type fulfilling the
|
||||
//! matrix interface of ISTL
|
||||
typedef M matrix_type;
|
||||
//! \brief The type of the domain.
|
||||
//!
|
||||
//! E.g. BlockVector or another type fulfilling the ISTL
|
||||
//! vector interface.
|
||||
typedef X domain_type;
|
||||
//! \brief The type of the range.
|
||||
//!
|
||||
//! E.g. BlockVector or another type fulfilling the ISTL
|
||||
//! vector interface.
|
||||
typedef Y range_type;
|
||||
//! \brief The field type of the range
|
||||
typedef typename X::field_type field_type;
|
||||
//! \brief The type of the communication object.
|
||||
//!
|
||||
//! This must either be OwnerOverlapCopyCommunication or a type
|
||||
//! implementing the same interface.
|
||||
typedef C communication_type;
|
||||
|
||||
enum {
|
||||
//! \brief The solver category.
|
||||
category=Dune::SolverCategory::overlapping
|
||||
};
|
||||
|
||||
/**
|
||||
* @brief constructor: just store a reference to a matrix.
|
||||
*
|
||||
* @param A The assembled matrix.
|
||||
* @param com The communication object for syncing overlap and copy
|
||||
* data points. (E.~g. OwnerOverlapCopyCommunication )
|
||||
*/
|
||||
OverlappingSchwarzOperator (const matrix_type& A, const communication_type& com)
|
||||
: _A_(A), communication(com)
|
||||
{}
|
||||
|
||||
//! apply operator to x: \f$ y = A(x) \f$
|
||||
virtual void apply (const X& x, Y& y) const
|
||||
{
|
||||
// Y y2(y.size());
|
||||
_A_.umv(x,y);
|
||||
communication.addOwnerOverlapToAll(y,y);
|
||||
}
|
||||
|
||||
//! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$
|
||||
virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
|
||||
{
|
||||
Y y2(y.size());
|
||||
_A_.usmv(alpha, x, y2);
|
||||
y = 0;
|
||||
communication.addOwnerOverlapToAll(y2,y);
|
||||
}
|
||||
|
||||
//! get the sequential assembled linear operator.
|
||||
virtual const matrix_type& getmat () const
|
||||
{
|
||||
return _A_;
|
||||
}
|
||||
|
||||
private:
|
||||
const matrix_type& _A_;
|
||||
const communication_type& communication;
|
||||
};
|
||||
|
||||
typedef OverlappingSchwarzOperator<ISTL::Mat,ISTL::Vec,ISTL::Vec,Dune::OwnerOverlapCopyCommunication<int,int>> ParMatrixAdapter;
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
|
||||
/*!
|
||||
\brief Class for ISTL solver parameters.
|
||||
\details It contains information about solver method, preconditioner
|
||||
and convergence criteria.
|
||||
*/
|
||||
|
||||
class ISTLSolParams
|
||||
{
|
||||
public:
|
||||
//! \brief Constructor.
|
||||
//! \param spar Linear solver parameters to use
|
||||
//! \param padm Process administrator to use
|
||||
ISTLSolParams(const LinSolParams& spar, const ProcessAdm& padm) :
|
||||
solParams(spar), adm(padm)
|
||||
{}
|
||||
|
||||
|
||||
//! \brief Setup solver and preconditioner.
|
||||
//! \param A Matrix to use
|
||||
//! \return tuple with (solver, preconditioner, op)
|
||||
std::tuple<std::unique_ptr<ISTL::InverseOperator>,
|
||||
std::unique_ptr<ISTL::Preconditioner>,
|
||||
std::unique_ptr<ISTL::Operator>>
|
||||
setupPC(ISTL::Mat& A);
|
||||
|
||||
//! \brief Obtain linear solver parameters.
|
||||
const LinSolParams& get() const { return solParams; }
|
||||
|
||||
protected:
|
||||
//! \brief Internal helper function for setting up a preconditioner.
|
||||
//! \param A Matrix to construct preconditioner for
|
||||
//! \param op Matrix adaptor to use (used with AMG)
|
||||
//! \param block Block to read settings from
|
||||
//! \param solver Solver object to instance. nullptr to do nothing.
|
||||
ISTL::Preconditioner* setupPCInternal(ISTL::Mat& A,
|
||||
ISTL::Operator& op,
|
||||
size_t block,
|
||||
std::unique_ptr<ISTL::InverseOperator>* solver);
|
||||
|
||||
const LinSolParams& solParams; //!< Reference to linear solver parameters.
|
||||
const ProcessAdm& adm; //!< Reference to process administrator.
|
||||
};
|
||||
|
||||
#endif
|
||||
76
src/LinAlg/ISTLSupport.h
Normal file
76
src/LinAlg/ISTLSupport.h
Normal file
@@ -0,0 +1,76 @@
|
||||
// $Id$
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file ISTLSupport.h
|
||||
//!
|
||||
//! \date May 12 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief IFEM ISTL support
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#ifndef _ISTL_SUPPORT_H_
|
||||
#define _ISTL_SUPPORT_H_
|
||||
|
||||
#include <vector>
|
||||
|
||||
#ifdef HAS_ISTL
|
||||
#include <dune/istl/bcrsmatrix.hh>
|
||||
#include <dune/istl/overlappingschwarz.hh>
|
||||
#include <dune/istl/preconditioners.hh>
|
||||
#include <dune/istl/solvers.hh>
|
||||
#include <dune/istl/schwarz.hh>
|
||||
#include <dune/istl/owneroverlapcopy.hh>
|
||||
#ifdef HAVE_SUPERLU
|
||||
#include <dune/istl/superlu.hh>
|
||||
#endif
|
||||
#ifdef HAVE_UMFPACK
|
||||
#include <dune/istl/umfpack.hh>
|
||||
#endif
|
||||
|
||||
namespace ISTL
|
||||
{
|
||||
typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>> Mat; //!< A sparse system matrix
|
||||
typedef Dune::BlockVector<Dune::FieldVector<double,1>> Vec; //!< A vector
|
||||
typedef Dune::MatrixAdapter<Mat,Vec,Vec> Operator; //!< A serial matrix operator
|
||||
typedef Dune::InverseOperator<Vec, Vec> InverseOperator; //!< Linear system inversion abstraction
|
||||
typedef Dune::Preconditioner<Vec,Vec> Preconditioner; //!< Preconditioner abstraction
|
||||
|
||||
/*! \brief Wrapper template to avoid memory leaks */
|
||||
template<template<class M> class Pre>
|
||||
class IOp2Pre : public Dune::InverseOperator2Preconditioner<Pre<ISTL::Mat>, Dune::SolverCategory::sequential> {
|
||||
typedef Dune::InverseOperator2Preconditioner<Pre<ISTL::Mat>, Dune::SolverCategory::sequential> SolverType;
|
||||
public:
|
||||
IOp2Pre(Pre<ISTL::Mat>* iop) : SolverType(*iop)
|
||||
{
|
||||
m_op.reset(iop);
|
||||
}
|
||||
protected:
|
||||
std::unique_ptr<Pre<ISTL::Mat>> m_op;
|
||||
};
|
||||
|
||||
#if defined(HAVE_UMFPACK)
|
||||
typedef Dune::UMFPack<ISTL::Mat> LUType;
|
||||
typedef IOp2Pre<Dune::UMFPack> LU;
|
||||
#elif defined(HAVE_SUPERLU)
|
||||
typedef Dune::SuperLU<ISTL::Mat> LUType;
|
||||
typedef IOp2Pre<Dune::SuperLU> LU;
|
||||
#endif
|
||||
|
||||
// Preconditioner types
|
||||
typedef Dune::SeqILU0<Mat, Vec, Vec> ILU0; //!< Sequential ILU(0)
|
||||
typedef Dune::SeqILUn<Mat, Vec, Vec> ILUn; //!< Sequential ILU(n)
|
||||
typedef Dune::SeqSOR<Mat, Vec, Vec> SOR; //!< Sequential SOR
|
||||
typedef Dune::SeqSSOR<Mat, Vec,Vec> SSOR; //!< Sequential SSOR
|
||||
typedef Dune::SeqJac<Mat, Vec, Vec> GJ; //!< Sequential Gauss-Jacobi
|
||||
typedef Dune::SeqGS<Mat, Vec, Vec> GS; //!< Sequential Gauss-Seidel
|
||||
typedef Dune::SeqOverlappingSchwarz<Mat, Vec,
|
||||
Dune::AdditiveSchwarzMode,
|
||||
LUType> ASMLU; //!< Sequential additive Schwarz w/ LU subdomain solvers
|
||||
typedef Dune::SeqOverlappingSchwarz<Mat, Vec> ASM; //!< Sequential additive Schwarz w/ ILU subdomain solvers
|
||||
}
|
||||
#endif
|
||||
|
||||
#endif
|
||||
@@ -13,12 +13,15 @@
|
||||
|
||||
#include "SystemMatrix.h"
|
||||
#include "DenseMatrix.h"
|
||||
#ifdef HAS_ISTL
|
||||
#include "ISTLMatrix.h"
|
||||
#endif
|
||||
#include "SPRMatrix.h"
|
||||
#include "SparseMatrix.h"
|
||||
#ifdef HAS_PETSC
|
||||
#include "LinSolParams.h"
|
||||
#include "PETScMatrix.h"
|
||||
#endif
|
||||
#include "LinSolParams.h"
|
||||
|
||||
|
||||
SystemVector* SystemVector::create (const ProcessAdm& adm, Type vectorType)
|
||||
@@ -26,6 +29,9 @@ SystemVector* SystemVector::create (const ProcessAdm& adm, Type vectorType)
|
||||
switch (vectorType)
|
||||
{
|
||||
case STD : return new StdVector();
|
||||
#ifdef HAS_ISTL
|
||||
case ISTL : return new ISTLVector(adm);
|
||||
#endif
|
||||
#ifdef HAS_PETSC
|
||||
case PETSC : return new PETScVector(adm);
|
||||
#endif
|
||||
@@ -73,6 +79,10 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType,
|
||||
if (matrixType == PETSC)
|
||||
return new PETScMatrix(padm,spar,ltype);
|
||||
#endif
|
||||
#ifdef HAS_ISTL
|
||||
if (matrixType == ISTL)
|
||||
return new ISTLMatrix(padm,spar,ltype);
|
||||
#endif
|
||||
|
||||
return SystemMatrix::create(padm,matrixType,ltype);
|
||||
}
|
||||
@@ -88,7 +98,16 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType,
|
||||
<< std::endl;
|
||||
exit(1);
|
||||
}
|
||||
#else
|
||||
#endif
|
||||
#ifndef HAS_ISTL
|
||||
if (matrixType == ISTL) {
|
||||
std::cerr <<"SystemMatrix::create: ISTL not compiled in, bailing out..."
|
||||
<< std::endl;
|
||||
exit(1);
|
||||
}
|
||||
#endif
|
||||
|
||||
#if defined(HAS_PETSC) || defined(HAS_ISTL)
|
||||
// Use default settings when no parameters are provided by user
|
||||
static LinSolParams defaultPar;
|
||||
#endif
|
||||
@@ -99,6 +118,9 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType,
|
||||
case SPR : return new SPRMatrix();
|
||||
case SPARSE: return new SparseMatrix(SparseMatrix::SUPERLU,num_thread_SLU);
|
||||
case SAMG : return new SparseMatrix(SparseMatrix::S_A_M_G);
|
||||
#ifdef HAS_ISTL
|
||||
case ISTL : return new ISTLMatrix(padm,defaultPar,ltype);
|
||||
#endif
|
||||
#ifdef HAS_PETSC
|
||||
case PETSC : return new PETScMatrix(padm,defaultPar,ltype);
|
||||
#endif
|
||||
|
||||
@@ -30,7 +30,7 @@ class SystemVector
|
||||
{
|
||||
public:
|
||||
//! \brief The available system vector formats.
|
||||
enum Type { STD = 0, PETSC = 1 };
|
||||
enum Type { STD = 0, PETSC = 1, ISTL = 2 };
|
||||
|
||||
//! \brief Static method creating a vector of the given type.
|
||||
static SystemVector* create(const ProcessAdm& padm, Type vectorType = STD);
|
||||
@@ -203,7 +203,7 @@ class SystemMatrix
|
||||
public:
|
||||
//! \brief The available system matrix formats.
|
||||
enum Type { DENSE = 0, SPR = 1, SPARSE = 2, SAMG = 3,
|
||||
PETSC = 4 };
|
||||
PETSC = 4, ISTL = 5 };
|
||||
|
||||
//! \brief Static method creating a matrix of the given type.
|
||||
static SystemMatrix* create(const ProcessAdm& padm, Type matrixType,
|
||||
|
||||
110
src/LinAlg/Test/MPI/TestISTLMatrix.C
Normal file
110
src/LinAlg/Test/MPI/TestISTLMatrix.C
Normal file
@@ -0,0 +1,110 @@
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file TestISTLMatrix.C
|
||||
//!
|
||||
//! \date Apr 27 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Unit tests for parallel ISTL matrices
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "AlgEqSystem.h"
|
||||
#include "ISTLMatrix.h"
|
||||
#include "ISTLSolParams.h"
|
||||
#include "SIM2D.h"
|
||||
#include "ASMs2D.h"
|
||||
#include "IntegrandBase.h"
|
||||
#include "SAM.h"
|
||||
#include "IFEM.h"
|
||||
|
||||
#include "gtest/gtest.h"
|
||||
|
||||
#include <fstream>
|
||||
|
||||
|
||||
typedef std::vector<int> IntVec;
|
||||
|
||||
static IntVec readIntVector(const std::string& file)
|
||||
{
|
||||
std::vector<int> result;
|
||||
std::ifstream f(file);
|
||||
size_t size;
|
||||
f >> size;
|
||||
result.resize(size);
|
||||
for (size_t j=0;j<size;++j)
|
||||
f >> result[j];
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
class DummyIntegrandForDummies : public IntegrandBase {
|
||||
public:
|
||||
DummyIntegrandForDummies(unsigned short int n = 0) : IntegrandBase(n) {}
|
||||
};
|
||||
|
||||
|
||||
class InspectMatrixSIM : public SIM2D {
|
||||
public:
|
||||
InspectMatrixSIM(unsigned char n1 = 2, bool check = false) :
|
||||
SIM2D(n1, check) { myProblem = new DummyIntegrandForDummies; }
|
||||
|
||||
InspectMatrixSIM(const std::vector<unsigned char>& n, bool check = false) :
|
||||
SIM2D(n, check) { myProblem = new DummyIntegrandForDummies; }
|
||||
|
||||
SystemMatrix* getMatrix() { return myEqSys->getMatrix(0); }
|
||||
SystemVector* getVector() { return myEqSys->getVector(0); }
|
||||
};
|
||||
|
||||
|
||||
TEST(TestISTLMatrix, AssembleMPI)
|
||||
{
|
||||
InspectMatrixSIM sim(1);
|
||||
sim.read("src/LinAlg/Test/refdata/petsc_test.xinp");
|
||||
sim.opt.solver = SystemMatrix::ISTL;
|
||||
sim.preprocess();
|
||||
sim.initSystem(SystemMatrix::ISTL);
|
||||
|
||||
Matrix stencil(4,4);
|
||||
stencil(1,1) = stencil(2,2) = stencil(3,3) = stencil(4,4) = 1.0;
|
||||
|
||||
for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
|
||||
sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
|
||||
|
||||
sim.getMatrix()->beginAssembly();
|
||||
sim.getMatrix()->endAssembly();
|
||||
|
||||
// now inspect the matrix
|
||||
const ProcessAdm& adm = sim.getProcessAdm();
|
||||
ISTL::Mat& mat = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
|
||||
ISTL::Vec b(mat.N()), b2(mat.N());
|
||||
|
||||
try {
|
||||
Dune::OwnerOverlapCopyCommunication<int,int> comm(*adm.getCommunicator());
|
||||
comm.indexSet().beginResize();
|
||||
typedef Dune::ParallelLocalIndex<Dune::OwnerOverlapCopyAttributeSet::AttributeSet> LI;
|
||||
for (size_t i = 0; i < adm.dd.getMLGEQ().size(); ++i) {
|
||||
int gid = adm.dd.getGlobalEq(i+1);
|
||||
comm.indexSet().add(gid-1, LI(i, gid >= adm.dd.getMinEq() ?
|
||||
Dune::OwnerOverlapCopyAttributeSet::owner :
|
||||
Dune::OwnerOverlapCopyAttributeSet::overlap));
|
||||
}
|
||||
comm.indexSet().endResize();
|
||||
comm.remoteIndices().setIncludeSelf(true);
|
||||
comm.remoteIndices().template rebuild<false>();
|
||||
|
||||
ISTL::ParMatrixAdapter op(mat, comm);
|
||||
|
||||
b = 1.0;
|
||||
op.apply(b, b2);
|
||||
} catch (Dune::ISTLError e) {
|
||||
std::cerr << e << std::endl;
|
||||
ASSERT_TRUE(false);
|
||||
}
|
||||
|
||||
IntVec v = readIntVector("src/LinAlg/Test/refdata/petsc_matrix_diagonal.ref");
|
||||
for (size_t i = 1; i <= adm.dd.getMLGEQ().size(); ++i)
|
||||
ASSERT_FLOAT_EQ(v[adm.dd.getGlobalEq(i)-1], b2[i-1]);
|
||||
}
|
||||
205
src/LinAlg/Test/TestISTLMatrix.C
Normal file
205
src/LinAlg/Test/TestISTLMatrix.C
Normal file
@@ -0,0 +1,205 @@
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file TestISTLMatrix.C
|
||||
//!
|
||||
//! \date Apr 27 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Unit tests for ISTL matrices
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "AlgEqSystem.h"
|
||||
#include "ISTLMatrix.h"
|
||||
#include "SIM2D.h"
|
||||
#include "ASMs2D.h"
|
||||
#include "IntegrandBase.h"
|
||||
#include "SAM.h"
|
||||
#include <dune/istl/io.hh>
|
||||
|
||||
#include "gtest/gtest.h"
|
||||
|
||||
#include <fstream>
|
||||
|
||||
|
||||
typedef std::vector<int> IntVec;
|
||||
|
||||
static IntVec readIntVector(const std::string& file)
|
||||
{
|
||||
std::vector<int> result;
|
||||
std::ifstream f(file);
|
||||
size_t size;
|
||||
f >> size;
|
||||
result.resize(size);
|
||||
for (size_t j=0;j<size;++j)
|
||||
f >> result[j];
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
class DummyIntegrandForDummies : public IntegrandBase {
|
||||
public:
|
||||
DummyIntegrandForDummies(unsigned short int n = 0) : IntegrandBase(n) {}
|
||||
};
|
||||
|
||||
|
||||
class InspectMatrixSIM : public SIM2D {
|
||||
public:
|
||||
InspectMatrixSIM(unsigned char n1 = 2, bool check = false) :
|
||||
SIM2D(n1, check) { myProblem = new DummyIntegrandForDummies; }
|
||||
|
||||
InspectMatrixSIM(const std::vector<unsigned char>& n, bool check = false) :
|
||||
SIM2D(n, check) { myProblem = new DummyIntegrandForDummies; }
|
||||
|
||||
SystemMatrix* getMatrix() { return myEqSys->getMatrix(0); }
|
||||
};
|
||||
|
||||
|
||||
class InspectBlockPreconditioner : public ISTL::BlockPreconditioner {
|
||||
public:
|
||||
InspectBlockPreconditioner(const ISTL::Mat& A,
|
||||
const DomainDecomposition& dd) :
|
||||
ISTL::BlockPreconditioner(A, dd, "upper") {}
|
||||
|
||||
void extractBlock(ISTL::Mat& B, const ISTL::Mat& A,
|
||||
const std::set<int>& eqs_row,
|
||||
const std::set<int>& eqs_col,
|
||||
const std::map<int,int>& eqs_col_g2l)
|
||||
{
|
||||
ISTL::BlockPreconditioner::extractBlock(B, A, eqs_row, eqs_col, eqs_col_g2l);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
TEST(TestISTLMatrix, Assemble)
|
||||
{
|
||||
InspectMatrixSIM sim(1);
|
||||
sim.read("src/LinAlg/Test/refdata/petsc_test.xinp");
|
||||
sim.opt.solver = SystemMatrix::ISTL;
|
||||
sim.preprocess();
|
||||
sim.initSystem(SystemMatrix::ISTL);
|
||||
|
||||
Matrix stencil(4,4);
|
||||
stencil(1,1) = stencil(2,2) = stencil(3,3) = stencil(4,4) = 1.0;
|
||||
|
||||
for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
|
||||
sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
|
||||
|
||||
sim.getMatrix()->beginAssembly();
|
||||
sim.getMatrix()->endAssembly();
|
||||
|
||||
// now inspect the matrix
|
||||
ISTL::Mat& mat = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
|
||||
IntVec v = readIntVector("src/LinAlg/Test/refdata/petsc_matrix_diagonal.ref");
|
||||
|
||||
ASSERT_EQ(mat.N(), v.size());
|
||||
|
||||
for (size_t i = 0; i < v.size(); ++i)
|
||||
ASSERT_FLOAT_EQ(double(v[i]), mat[i][i]);
|
||||
}
|
||||
|
||||
|
||||
TEST(TestISTLMatrix, AssembleBasisBlocks)
|
||||
{
|
||||
InspectMatrixSIM sim({1,1});
|
||||
sim.read("src/LinAlg/Test/refdata/petsc_test_blocks_basis.xinp");
|
||||
sim.opt.solver = SystemMatrix::ISTL;
|
||||
sim.preprocess();
|
||||
sim.initSystem(SystemMatrix::ISTL);
|
||||
|
||||
Matrix stencil(13,13);
|
||||
for (size_t i = 1; i<= 9; ++i)
|
||||
stencil(i,i) = 1.0;
|
||||
for (size_t i = 10; i<= 13; ++i)
|
||||
stencil(i,i) = 2.0;
|
||||
|
||||
for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
|
||||
sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
|
||||
|
||||
sim.getMatrix()->beginAssembly();
|
||||
sim.getMatrix()->endAssembly();
|
||||
|
||||
const ProcessAdm& adm = sim.getProcessAdm();
|
||||
ISTL::Mat& A = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
|
||||
InspectBlockPreconditioner block(A, adm.dd);
|
||||
|
||||
// now inspect the matrix blocks
|
||||
std::vector<ISTL::Mat> mat(5);
|
||||
block.extractBlock(mat[0], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
|
||||
block.extractBlock(mat[1], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
|
||||
block.extractBlock(mat[2], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
|
||||
block.extractBlock(mat[3], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
|
||||
mat[4] = block.getBlock(1); // S = D
|
||||
|
||||
for (size_t b = 0; b < mat.size(); ++b) {
|
||||
if (b == 0 || b >= 3) { // diagonal blocks
|
||||
std::stringstream str;
|
||||
str << "src/LinAlg/Test/refdata/petsc_matrix_diagonal_basis_block" << std::min(b/2 + 1, 2ul) << ".ref";
|
||||
IntVec v = readIntVector(str.str());
|
||||
ASSERT_EQ(v.size(), mat[b].N());
|
||||
for (size_t i = 0; i < v.size(); ++i)
|
||||
ASSERT_FLOAT_EQ(double(v[i]), mat[b][i][i]);
|
||||
}
|
||||
|
||||
// check that no values outside the diagonal are != 0
|
||||
for (auto r = mat[b].begin(); r != mat[b].end(); ++r) {
|
||||
for (auto it = r->begin(); it != r->end(); ++it)
|
||||
if (r.index() != it.index())
|
||||
ASSERT_FLOAT_EQ(*it, 0.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
TEST(TestISTLMatrix, AssembleComponentBlocks)
|
||||
{
|
||||
InspectMatrixSIM sim(2);
|
||||
sim.read("src/LinAlg/Test/refdata/petsc_test_blocks_components.xinp");
|
||||
sim.opt.solver = SystemMatrix::ISTL;
|
||||
sim.preprocess();
|
||||
sim.initSystem(SystemMatrix::ISTL);
|
||||
|
||||
Matrix stencil(4*2,4*2);
|
||||
for (size_t i = 1; i<= 4; ++i) {
|
||||
stencil(2*i-1,2*i-1) = 1.0;
|
||||
stencil(2*i,2*i) = 2.0;
|
||||
}
|
||||
|
||||
for (int iel = 1; iel <= sim.getSAM()->getNoElms(); ++iel)
|
||||
sim.getMatrix()->assemble(stencil, *sim.getSAM(), iel);
|
||||
|
||||
sim.getMatrix()->beginAssembly();
|
||||
sim.getMatrix()->endAssembly();
|
||||
|
||||
const ProcessAdm& adm = sim.getProcessAdm();
|
||||
ISTL::Mat& A = static_cast<ISTLMatrix*>(sim.getMatrix())->getMatrix();
|
||||
InspectBlockPreconditioner block(A, adm.dd);
|
||||
|
||||
// now inspect the matrix blocks
|
||||
std::vector<ISTL::Mat> mat(5);
|
||||
block.extractBlock(mat[0], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
|
||||
block.extractBlock(mat[1], A, adm.dd.getBlockEqs(0), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
|
||||
block.extractBlock(mat[2], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(0), adm.dd.getG2LEQ(0));
|
||||
block.extractBlock(mat[3], A, adm.dd.getBlockEqs(1), adm.dd.getBlockEqs(1), adm.dd.getG2LEQ(1));
|
||||
mat[4] = block.getBlock(1); // S = D
|
||||
|
||||
for (size_t b = 0; b < mat.size(); ++b) {
|
||||
if (b == 0 || b >= 3) { // diagonal blocks
|
||||
std::stringstream str;
|
||||
str << "src/LinAlg/Test/refdata/petsc_matrix_diagonal_components_block" << std::min(b/2 + 1, 2ul) << ".ref";
|
||||
IntVec v = readIntVector(str.str());
|
||||
ASSERT_EQ(v.size(), mat[b].N());
|
||||
for (size_t i = 0; i < v.size(); ++i)
|
||||
ASSERT_FLOAT_EQ(double(v[i]), mat[b][i][i]);
|
||||
}
|
||||
|
||||
// check that no values outside the diagonal are != 0
|
||||
for (auto r = mat[b].begin(); r != mat[b].end(); ++r) {
|
||||
for (auto it = r->begin(); it != r->end(); ++it)
|
||||
if (r.index() != it.index())
|
||||
ASSERT_FLOAT_EQ(*it, 0.0);
|
||||
}
|
||||
}
|
||||
}
|
||||
82
src/LinAlg/Test/TestISTLPETScMatrix.C
Normal file
82
src/LinAlg/Test/TestISTLPETScMatrix.C
Normal file
@@ -0,0 +1,82 @@
|
||||
//==============================================================================
|
||||
//!
|
||||
//! \file TestISTLPETScMatrix.C
|
||||
//!
|
||||
//! \date Apr 27 2016
|
||||
//!
|
||||
//! \author Arne Morten Kvarving / SINTEF
|
||||
//!
|
||||
//! \brief Unit tests comparing ISTL and PETSc matrices
|
||||
//!
|
||||
//==============================================================================
|
||||
|
||||
#include "AlgEqSystem.h"
|
||||
#include "ISTLMatrix.h"
|
||||
#include "PETScMatrix.h"
|
||||
#include "SIM2D.h"
|
||||
#include "ASMs2D.h"
|
||||
#include "IntegrandBase.h"
|
||||
#include "SAM.h"
|
||||
#include <dune/istl/io.hh>
|
||||
|
||||
#include "gtest/gtest.h"
|
||||
|
||||
#include <fstream>
|
||||
|
||||
|
||||
class DummyIntegrandForDummies : public IntegrandBase {
|
||||
public:
|
||||
DummyIntegrandForDummies(unsigned short int n = 0) : IntegrandBase(n) {}
|
||||
};
|
||||
|
||||
|
||||
class InspectMatrixSIM : public SIM2D {
|
||||
public:
|
||||
InspectMatrixSIM() :
|
||||
SIM2D({1,1}, false) { myProblem = new DummyIntegrandForDummies; }
|
||||
|
||||
SystemMatrix* getMatrix() { return myEqSys->getMatrix(0); }
|
||||
};
|
||||
|
||||
|
||||
|
||||
TEST(TestISTLPETScMatrix, SchurComplement)
|
||||
{
|
||||
Matrix stencil(13,13);
|
||||
for (size_t i = 1; i<= 13; ++i)
|
||||
for (size_t j = 1; j <= 13; ++j)
|
||||
stencil(i,j) = 1.0;
|
||||
|
||||
std::array<InspectMatrixSIM,2> sim;
|
||||
for (size_t i = 0; i < 2; ++i) {
|
||||
sim[i].read("src/LinAlg/Test/refdata/petsc_test_blocks_basis.xinp");
|
||||
sim[i].opt.solver = i == 0 ? SystemMatrix::PETSC : SystemMatrix::ISTL;
|
||||
sim[i].preprocess();
|
||||
sim[i].initSystem(i == 0 ? SystemMatrix::PETSC : SystemMatrix::ISTL);
|
||||
|
||||
for (int iel = 1; iel <= sim[i].getSAM()->getNoElms(); ++iel)
|
||||
sim[i].getMatrix()->assemble(stencil, *sim[i].getSAM(), iel);
|
||||
|
||||
sim[i].getMatrix()->beginAssembly();
|
||||
sim[i].getMatrix()->endAssembly();
|
||||
}
|
||||
|
||||
const ProcessAdm& adm = sim[1].getProcessAdm();
|
||||
ISTL::Mat& A = static_cast<ISTLMatrix*>(sim[1].getMatrix())->getMatrix();
|
||||
ISTL::BlockPreconditioner block(A, adm.dd, "upper");
|
||||
|
||||
ISTL::Mat& S = block.getBlock(1);
|
||||
PETScSolParams params(LinSolParams(), adm);
|
||||
params.setupSchurComplement(static_cast<PETScMatrix*>(sim[0].getMatrix())->getBlockMatrices());
|
||||
|
||||
// check that matrices are the same
|
||||
for (size_t r = 0; r < S.N(); ++r) {
|
||||
const PetscInt* cols;
|
||||
PetscInt ncols;
|
||||
const PetscScalar* vals;
|
||||
MatGetRow(params.getSchurComplement(), r, &ncols, &cols, &vals);
|
||||
for (PetscInt i = 0; i < ncols; ++i)
|
||||
ASSERT_FLOAT_EQ(vals[i], S[r][cols[i]]);
|
||||
MatRestoreRow(params.getSchurComplement(), r, &ncols, &cols, &vals);
|
||||
}
|
||||
}
|
||||
@@ -68,6 +68,8 @@ void SIMoptions::setLinearSolver (const std::string& eqsolver)
|
||||
solver = SystemMatrix::SAMG;
|
||||
else if (eqsolver == "petsc")
|
||||
solver = SystemMatrix::PETSC;
|
||||
else if (eqsolver == "istl")
|
||||
solver = SystemMatrix::ISTL;
|
||||
}
|
||||
|
||||
|
||||
@@ -223,6 +225,8 @@ bool SIMoptions::parseOldOptions (int argc, char** argv, int& i)
|
||||
solver = SystemMatrix::SAMG;
|
||||
else if (!strcmp(argv[i],"-petsc"))
|
||||
solver = SystemMatrix::PETSC;
|
||||
else if (!strcmp(argv[i],"-istl"))
|
||||
solver = SystemMatrix::ISTL;
|
||||
else if (!strncmp(argv[i],"-lag",4))
|
||||
discretization = ASM::Lagrange;
|
||||
else if (!strncmp(argv[i],"-spec",5))
|
||||
|
||||
Reference in New Issue
Block a user