From f2044db3a43d73504943aa9904d105e7e739f124 Mon Sep 17 00:00:00 2001
From: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue, 6 May 2014 15:21:57 +0200
Subject: [PATCH] Squashed commit of the following:
commit 920fe22f2dd8e4aa5e1e1a012474e254387cfe51
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Tue May 6 15:12:54 2014 +0200
Complete PETSc solver wrapper
commit ac4505aa7890d944c75b28ed8154a56ca59b9f1e
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon May 5 18:13:11 2014 +0200
New PETSc solver test
commit c88efb737e512f7c9360a7ce3be1257cf7dc8d9f
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon May 5 11:41:07 2014 +0200
Classify KSP and PC options for PETSc
commit a1575c81a736e8a2a26dd9b7cb10e3ef2a951720
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri May 2 17:18:42 2014 +0200
Split files of PETSc options
commit 726e93a9db2350339738ee2bf903a4ab01263642
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri May 2 17:01:04 2014 +0200
LU configuration (compiled)
commit 3d9c7f3fd3ab2285c76dbff928c1b65febe462f1
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Fri May 2 16:52:03 2014 +0200
General template for PETSc solver configuration (compiled)
commit ac64c755b575ee448dd214f9e9fdc4da8db21147
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Wed Apr 30 18:03:38 2014 +0200
Organize the option data
commit 5abf2fd64be950fc675901f06e03da52efa0e97c
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon Apr 28 18:13:19 2014 +0200
Add a template funtion to set PC options.
commit 6d17b57645084ba39b0d513153058efd0d7aceb2
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Mon Apr 28 17:01:45 2014 +0200
Add ILU and SOR options
commit ee5f50187e6952c043741d38905bf424841d8345
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu Apr 24 18:03:48 2014 +0200
Add more solver option
commit bba799fc6b284230cf759c83c9bfcd98683eae78
Author: Wenqing Wang <wenqing.wang@ufz.de>
Date: Thu Apr 17 17:17:30 2014 +0200
Modify the PETSc linear solver class
---
FileIO/RapidXmlIO/RapidVtuInterface.cpp | 4 +-
FileIO/RapidXmlIO/RapidVtuInterface.h | 2 +-
MathLib/CMakeLists.txt | 4 +
.../PETScPC_KSP_Chebyshev_Option.cpp | 34 ++
.../KSP_Option/PETScPC_KSP_Chebyshev_Option.h | 52 ++++
.../KSP_Option/PETScPC_KSP_GMRES_Option.cpp | 65 ++++
.../KSP_Option/PETScPC_KSP_GMRES_Option.h | 61 ++++
.../PETScPC_KSP_Richards_Option.cpp | 31 ++
.../KSP_Option/PETScPC_KSP_Richards_Option.h | 50 +++
.../PETSc/PC_Option/PETScPC_AMG_Option.cpp | 42 +++
.../PETSc/PC_Option/PETScPC_AMG_Option.h | 61 ++++
.../PETSc/PC_Option/PETScPC_ASM_Option.cpp | 46 +++
.../PETSc/PC_Option/PETScPC_ASM_Option.h | 54 ++++
.../PETSc/PC_Option/PETScPC_ILU_Option.cpp | 69 +++++
.../PETSc/PC_Option/PETScPC_ILU_Option.h | 65 ++++
.../PETSc/PC_Option/PETScPC_LU_Option.cpp | 54 ++++
.../PETSc/PC_Option/PETScPC_LU_Option.h | 50 +++
.../PETSc/PC_Option/PETScPC_SOR_Option.cpp | 67 ++++
.../PETSc/PC_Option/PETScPC_SOR_Option.h | 66 ++++
MathLib/LinAlg/PETSc/PETScLinearSolver.cpp | 183 +++++++++++
MathLib/LinAlg/PETSc/PETScLinearSolver.h | 83 +++++
.../LinAlg/PETSc/PETScLinearSolverOption.cpp | 114 +++++++
.../LinAlg/PETSc/PETScLinearSolverOption.h | 161 ++++++++++
MathLib/LinAlg/PETSc/PETScMatrix.cpp | 4 +-
MathLib/LinAlg/PETSc/PETScTools.cpp | 41 +++
MathLib/LinAlg/PETSc/PETScTools.h | 45 +++
Tests/MathLib/TestLinearSolver.cpp | 292 +++++++++++++++++-
27 files changed, 1786 insertions(+), 14 deletions(-)
create mode 100644 MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.h
create mode 100644 MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.h
create mode 100644 MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.h
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.h
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.h
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.h
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.h
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.cpp
create mode 100644 MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.h
create mode 100644 MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
create mode 100644 MathLib/LinAlg/PETSc/PETScLinearSolver.h
create mode 100644 MathLib/LinAlg/PETSc/PETScLinearSolverOption.cpp
create mode 100644 MathLib/LinAlg/PETSc/PETScLinearSolverOption.h
create mode 100644 MathLib/LinAlg/PETSc/PETScTools.cpp
create mode 100644 MathLib/LinAlg/PETSc/PETScTools.h
diff --git a/FileIO/RapidXmlIO/RapidVtuInterface.cpp b/FileIO/RapidXmlIO/RapidVtuInterface.cpp
index a4c472cb31a..647a97cf8f6 100644
--- a/FileIO/RapidXmlIO/RapidVtuInterface.cpp
+++ b/FileIO/RapidXmlIO/RapidVtuInterface.cpp
@@ -316,14 +316,14 @@ bool RapidVtuInterface::isVTKUnstructuredGrid(const rapidxml::xml_node<>* vtk_ro
return false;
}
-unsigned char* RapidVtuInterface::uncompressData(const rapidxml::xml_node<>* node)
+char* RapidVtuInterface::uncompressData(const rapidxml::xml_node<>* node)
{
rapidxml::xml_node<>* data_node = node->first_node("AppendedData");
char* compressed_data (NULL);
if (data_node)
compressed_data = data_node->value();
- return NULL; // Does here return compressed_data?
+ return compressed_data; //? return NULL; // Does here return compressed_data?
}
diff --git a/FileIO/RapidXmlIO/RapidVtuInterface.h b/FileIO/RapidXmlIO/RapidVtuInterface.h
index 7a38b776fd8..1c95e616365 100644
--- a/FileIO/RapidXmlIO/RapidVtuInterface.h
+++ b/FileIO/RapidXmlIO/RapidVtuInterface.h
@@ -77,7 +77,7 @@ private:
/// Construct an Element-object from the data given to the method and the data at the current stream position.
static MeshLib::Element* readElement(std::stringstream &iss, const std::vector<MeshLib::Node*> &nodes, unsigned material, unsigned type);
- static unsigned char* uncompressData(const rapidxml::xml_node<>* node);
+ static char* uncompressData(const rapidxml::xml_node<>* node);
bool _use_compressor;
};
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index 4e01547a2ef..2bb0ed21eea 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -31,6 +31,10 @@ ENDIF ()
IF (OGS_USE_PETSC)
GET_SOURCE_FILES(SOURCES_LINALG_PETSC LinAlg/PETSc)
SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_PETSC})
+ GET_SOURCE_FILES(SOURCES_LINALG_PETSC_KSP LinAlg/PETSc/KSP_Option)
+ SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_PETSC_KSP})
+ GET_SOURCE_FILES(SOURCES_LINALG_PETSC_PC LinAlg/PETSc/PC_Option)
+ SET ( SOURCES ${SOURCES} ${SOURCES_LINALG_PETSC_PC})
ENDIF ()
IF (METIS_FOUND)
diff --git a/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.cpp b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.cpp
new file mode 100644
index 00000000000..7e06d701d0e
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.cpp
@@ -0,0 +1,34 @@
+/*!
+ \file PETScPC_KSP_Chebyshev_Option.cpp
+ \brief Define the configuration data for the PETSc Chebyshev linear solver.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_KSP_Chebyshev_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_KSP_Chebyshev_Option::
+PETScPC_KSP_Chebyshev_Option(const boost::property_tree::ptree &option)
+ : emin_chebyshev(0.01), emax_chebyshev(100.0)
+{
+ auto val = option.get_optional<double>("smallest_eignvalue");
+ emin_chebyshev = *val;
+
+ val = option.get_optional<double>("maximum_eignvalue");
+ emax_chebyshev = *val;
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.h b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.h
new file mode 100644
index 00000000000..6433846f6af
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Chebyshev_Option.h
@@ -0,0 +1,52 @@
+/*!
+ \file PETScPC_KSP_Chebyshev_Option.h
+ \brief Define the configuration data for the PETSc Chebyshev linear solver.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCPC_KSP_CHEBYSHEV_OPTION_H_
+#define PETSCPC_KSP_CHEBYSHEV_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+/*!
+ PETSc KSP Chebyshev options
+
+ */
+struct PETScPC_KSP_Chebyshev_Option
+{
+ PETScPC_KSP_Chebyshev_Option(const boost::property_tree::ptree &option);
+
+ /// Set Chebyshev option
+ void setOption(KSP &solver)
+ {
+ KSPChebyshevSetEigenvalues(solver, emax_chebyshev, emin_chebyshev);
+ }
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_KSP_Chebyshev_Option& opt)
+ {
+ emin_chebyshev = opt.emin_chebyshev;
+ emax_chebyshev = opt.emax_chebyshev;
+ }
+
+ /// Smallest eignvalue for Chebyshev.
+ PetscReal emin_chebyshev;
+ /// maximum eignvalue for Chebyshev.
+ PetscReal emax_chebyshev;
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.cpp b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.cpp
new file mode 100644
index 00000000000..5dfe3499f72
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.cpp
@@ -0,0 +1,65 @@
+/*!
+ \file PETScPC_KSP_GMRES_Option.cpp
+ \brief Define the configuration data for the PETSc GMRES linear solver.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_KSP_GMRES_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_KSP_GMRES_Option::
+PETScPC_KSP_GMRES_Option(const boost::property_tree::ptree &option)
+ : restart_number_gmres(30), is_modified_gram_schmidt_gmres(false),
+ refine_type_gmres(KSP_GMRES_CGS_REFINE_NEVER)
+{
+ auto val = option.get_optional<double>("restart_number");
+ restart_number_gmres = *val;
+
+ boost::optional<bool> bool_vals = option.get_optional<bool>("is_modified_ram_schmidt_orthog");
+ is_modified_gram_schmidt_gmres = *bool_vals;
+
+ auto refine_type = option.get_optional<int>("refine_type");
+ switch(*refine_type)
+ {
+ case 0:
+ refine_type_gmres = KSP_GMRES_CGS_REFINE_NEVER;
+ break;
+ case 1:
+ refine_type_gmres = KSP_GMRES_CGS_REFINE_IFNEEDED;
+ break;
+ case 2:
+ refine_type_gmres = KSP_GMRES_CGS_REFINE_ALWAYS;
+ break;
+ default:
+ refine_type_gmres = KSP_GMRES_CGS_REFINE_NEVER;
+ break;
+ }
+}
+
+/// Set Chebyshev option
+void PETScPC_KSP_GMRES_Option::setOption(KSP &solver)
+{
+ KSPGMRESSetRestart(solver, restart_number_gmres);
+
+ if(is_modified_gram_schmidt_gmres)
+ {
+ KSPGMRESSetOrthogonalization(solver, KSPGMRESClassicalGramSchmidtOrthogonalization);
+ }
+
+ KSPGMRESSetCGSRefinementType(solver, refine_type_gmres);
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.h b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.h
new file mode 100644
index 00000000000..7aacabf8030
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_GMRES_Option.h
@@ -0,0 +1,61 @@
+/*!
+ \file PETScPC_KSP_GMRES_Option.h
+ \brief Define the configuration data for the PETSc GMRES linear solver.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCPC_KSP_GMRES_OPTION_H_
+#define PETSCPC_KSP_GMRES_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+/*!
+ PETSc KSP GMRES options
+
+ */
+struct PETScPC_KSP_GMRES_Option
+{
+ PETScPC_KSP_GMRES_Option(const boost::property_tree::ptree &option);
+
+ /// Set GMRES option
+ void setOption(KSP &solver);
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_KSP_GMRES_Option& opt)
+ {
+ restart_number_gmres = opt.restart_number_gmres;
+ is_modified_gram_schmidt_gmres = opt.is_modified_gram_schmidt_gmres;
+ refine_type_gmres = opt.refine_type_gmres;
+ }
+
+ /// Restart number of GMRES.
+ PetscInt restart_number_gmres;
+
+ /// Flag for the modified Gram-Schmidt orthogonalization.
+ bool is_modified_gram_schmidt_gmres;
+
+ /*!
+ \brief Refinement type for GMRES.
+ This iterative refinement is used to improve the stability
+ of orthogonalization.
+ KSPGMRESCGSRefinementType is a enum type of PETSc defined as
+ typedef enum {KSP_GMRES_CGS_REFINE_NEVER, KSP_GMRES_CGS_REFINE_IFNEEDED,
+ KSP_GMRES_CGS_REFINE_ALWAYS} KSPGMRESCGSRefinementType;
+ */
+ KSPGMRESCGSRefinementType refine_type_gmres;
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.cpp b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.cpp
new file mode 100644
index 00000000000..378f2984316
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.cpp
@@ -0,0 +1,31 @@
+/*!
+ \file PETScPC_KSP_Richards_Option.cpp
+ \brief Define the configuration data for the PETSc Richards linear solver.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_KSP_Richards_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_KSP_Richards_Option::
+PETScPC_KSP_Richards_Option(const boost::property_tree::ptree &option)
+ : damping_factor_richards(1.0)
+{
+ auto damping_factor = option.get_optional<double>("damping_factor");
+ damping_factor_richards = *damping_factor;
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.h b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.h
new file mode 100644
index 00000000000..2610e3f25c8
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/KSP_Option/PETScPC_KSP_Richards_Option.h
@@ -0,0 +1,50 @@
+/*!
+ \file PETScPC_KSP_Richards_Option.h
+ \brief Define the configuration data for the PETSc Richards linear solver.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCPC_KSP_RICHARDS_OPTION_H_
+#define PETSCPC_KSP_RICHARDS_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+
+/*!
+ PETSc KSP Richards options
+
+ */
+struct PETScPC_KSP_Richards_Option
+{
+ PETScPC_KSP_Richards_Option(const boost::property_tree::ptree &option);
+
+ /// Set Richards option
+ void setOption(KSP &solver)
+ {
+ KSPRichardsonSetScale(solver, damping_factor_richards);
+ }
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_KSP_Richards_Option& opt)
+ {
+ damping_factor_richards = opt.damping_factor_richards;
+ }
+
+ /// Damping factor for Richards.
+ PetscReal damping_factor_richards;
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.cpp b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.cpp
new file mode 100644
index 00000000000..2b703134924
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.cpp
@@ -0,0 +1,42 @@
+/*!
+ \file PETScPC_AMG_Option.cpp
+ \brief Define the configuration data for PETSc AMG preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_AMG_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_AMG_Option::
+PETScPC_AMG_Option(const boost::property_tree::ptree &option)
+ : nsmooths(1), type(PCGAMGAGG), is_agg(true)
+{
+ auto val = option.get_optional<double>("agg_nsmooths");
+ nsmooths = *val;
+
+ auto type_name = option.get_optional<std::string>("type");
+ if(type_name->find("agg") != std::string::npos)
+ {
+ type = PCGAMGAGG;
+ }
+ if(type_name->find("geo") != std::string::npos)
+ {
+ type = PCGAMGGEO;
+ is_agg = false;
+ }
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.h b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.h
new file mode 100644
index 00000000000..cc57b78e8f4
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_AMG_Option.h
@@ -0,0 +1,61 @@
+/*!
+ \file PETScPC_AMG_Option.h
+ \brief Define the configuration data for PETSc AMG preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCPC_AMG_OPTION_H_
+#define PETSCPC_AMG_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+/*!
+ PETSc AMG preconditioner options
+
+*/
+struct PETScPC_AMG_Option
+{
+ PETScPC_AMG_Option(const boost::property_tree::ptree &option);
+
+ /// Set AMG option
+ void setOption(PC &pc)
+ {
+ PCGAMGSetType(pc, type);
+ if(is_agg)
+ {
+ PCGAMGSetNSmooths(pc, nsmooths);
+ }
+ }
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_AMG_Option& opt)
+ {
+ nsmooths = opt.nsmooths;
+ type = opt.type;
+ }
+
+ /// Number of smoothing steps
+ PetscInt nsmooths;
+
+ /*!
+ AMG type: PCGAMGAGG or PCGAMGGEO
+ */
+ PCGAMGType type;
+
+ /// Flag of agg type
+ bool is_agg;
+};
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.cpp b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.cpp
new file mode 100644
index 00000000000..9f837fd2c19
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.cpp
@@ -0,0 +1,46 @@
+/*!
+ \file PETScPC_ASM_Option.cpp
+ \brief Define the configuration data for PETSc ASM preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_ASM_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_ASM_Option::
+PETScPC_ASM_Option(const boost::property_tree::ptree &option)
+ : type(PC_ASM_BASIC)
+{
+ auto type_name = option.get_optional<std::string>("type");
+ if(type_name->find("base") != std::string::npos)
+ {
+ type = PC_ASM_BASIC;
+ }
+ if(type_name->find("interpolate") != std::string::npos)
+ {
+ type = PC_ASM_INTERPOLATE;
+ }
+ if(type_name->find("restrict") != std::string::npos)
+ {
+ type = PC_ASM_RESTRICT;
+ }
+ if(type_name->find("PC_ASM_NONE") != std::string::npos)
+ {
+ type = PC_ASM_NONE;
+ }
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.h b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.h
new file mode 100644
index 00000000000..dc7001328c0
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ASM_Option.h
@@ -0,0 +1,54 @@
+/*!
+ \file PETScPC_ASM_Option.h
+ \brief Define the configuration data for PETSc ASM preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCPC_ASM_OPTION_H_
+#define PETSCPC_ASM_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+/*!
+ PETSc ASM preconditioner options
+
+*/
+struct PETScPC_ASM_Option
+{
+ PETScPC_ASM_Option(const boost::property_tree::ptree &option);
+
+ /// Set ASM option
+ void setOption(PC &pc)
+ {
+ PCASMSetType(pc, type);
+ }
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_ASM_Option& opt)
+ {
+ type = opt.type;
+ }
+
+ /*!
+ ASM type:
+ PC_ASM_BASIC,
+ PC_ASM_INTERPOLATE,
+ PC_ASM_RESTRICT,
+ PC_ASM_NONE
+ */
+ PCASMType type;
+};
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.cpp b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.cpp
new file mode 100644
index 00000000000..e73cab80108
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.cpp
@@ -0,0 +1,69 @@
+/*!
+ \file PETScPC_ILU_Option.cpp
+ \brief Define the configuration data for PETSc ILU preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_ILU_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_ILU_Option::
+PETScPC_ILU_Option(const boost::property_tree::ptree &option)
+ : levels(PETSC_DECIDE), reuse_ordering(false),
+ reuse_fill(false), use_in_place(false), allow_diagonal_fill(false)
+{
+ auto val = option.get_optional<double>("levels");
+ levels = *val;
+
+ auto reuse_order = option.get_optional<bool>("reuse_ordering");
+ reuse_ordering = *reuse_order;
+
+ auto rfill = option.get_optional<bool>("reuse_fill");
+ reuse_fill = *rfill;
+
+ auto inplane = option.get_optional<bool>("use_in_place");
+ use_in_place = *inplane;
+
+ auto low_diag_fill = option.get_optional<bool>("allow_diagonal_fill");
+ allow_diagonal_fill = *low_diag_fill;
+}
+
+void PETScPC_ILU_Option::setOption(PC &pc)
+{
+ PCFactorSetLevels(pc, levels);
+
+ if(reuse_ordering)
+ {
+ PCFactorSetReuseOrdering(pc, PETSC_TRUE);
+ }
+
+ if(reuse_fill)
+ {
+ PCFactorSetReuseFill(pc, PETSC_TRUE);
+ }
+
+ if(use_in_place)
+ {
+ PCFactorSetUseInPlace(pc);
+ }
+
+ if(allow_diagonal_fill)
+ {
+ PCFactorSetAllowDiagonalFill(pc);
+ }
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.h b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.h
new file mode 100644
index 00000000000..e5ebae6b398
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_ILU_Option.h
@@ -0,0 +1,65 @@
+/*!
+ \file PETScPC_ILU_Option.h
+ \brief Define the configuration data for PETSc ILU preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCPC_ILU_OPTION_H_
+#define PETSCPC_ILU_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+
+/*!
+ PETSc ILU preconditioner options
+
+ */
+struct PETScPC_ILU_Option
+{
+ PETScPC_ILU_Option(const boost::property_tree::ptree &option);
+
+ /// Set ILU option
+ void setOption(PC &pc);
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_ILU_Option& opt)
+ {
+ levels = opt.levels;
+ reuse_ordering = opt.reuse_ordering;
+ reuse_fill = opt.reuse_fill;
+ use_in_place = opt.use_in_place;
+ allow_diagonal_fill = opt.allow_diagonal_fill;
+ }
+
+ /// Number of levels of fill for ILU(k)
+ int levels;
+
+ /// Reuse ordering of factorized matrix from previous factorization
+ bool reuse_ordering;
+
+ /// Use the fill ratio computed in the initial factorization.
+ bool reuse_fill;
+
+ /*! for ILU(0) with natural ordering, reuses the space of the matrix
+ for its factorization (overwrites original matrix)
+ */
+ bool use_in_place;
+
+ /// fill in a zero diagonal even if levels of fill indicate it wouldn't be fill
+ bool allow_diagonal_fill;
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.cpp b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.cpp
new file mode 100644
index 00000000000..149d353671f
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.cpp
@@ -0,0 +1,54 @@
+/*!
+ \file PETScPC_LU_Option.cpp
+ \brief Define the configuration data for LU decompoistion preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_LU_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_LU_Option::
+PETScPC_LU_Option(const boost::property_tree::ptree &option)
+ : mat_type(MATORDERINGNATURAL)
+{
+ auto type_name = option.get_optional<std::string>("type");
+ if(type_name->find("natural") != std::string::npos)
+ {
+ mat_type = MATORDERINGNATURAL;
+ }
+ if(type_name->find("nd") != std::string::npos)
+ {
+ mat_type = MATORDERINGND;
+ }
+ if(type_name->find("1wd") != std::string::npos)
+ {
+ mat_type = MATORDERING1WD;
+ }
+ if(type_name->find("rcm") != std::string::npos)
+ {
+ mat_type = MATORDERINGRCM;
+ }
+ if(type_name->find("qmd") != std::string::npos)
+ {
+ mat_type = MATORDERINGQMD;
+ }
+ if(type_name->find("rowlength") != std::string::npos)
+ {
+ mat_type = MATORDERINGROWLENGTH;
+ }
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.h b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.h
new file mode 100644
index 00000000000..07c14471ca8
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_LU_Option.h
@@ -0,0 +1,50 @@
+/*!
+ \file PETScPC_LU_Option.h
+ \brief Define the configuration data for LU decompoistion preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSSPC_LU_OPTION_H_
+#define PETSSPC_LU_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+
+/*!
+ PETSc LU preconditioner options
+
+*/
+struct PETScPC_LU_Option
+{
+ PETScPC_LU_Option(const boost::property_tree::ptree &option);
+
+ /// Set LU option
+ void setOption(PC &pc)
+ {
+ PCFactorSetMatOrderingType(pc, mat_type);
+ }
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_LU_Option& opt)
+ {
+ mat_type = opt.mat_type;
+ }
+
+ /// Ordering routine (to reduce fill) to be used in the LU factorization.
+ MatOrderingType mat_type;
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.cpp b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.cpp
new file mode 100644
index 00000000000..21f9ea19d3e
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.cpp
@@ -0,0 +1,67 @@
+/*!
+ \file PETScPC_SOR_Option.cpp
+ \brief Define the configuration data for PETSc SOR/SSOR preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScPC_SOR_Option.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScPC_SOR_Option::
+PETScPC_SOR_Option(const boost::property_tree::ptree &option)
+ : omega(1.), its(PETSC_DEFAULT),
+ lits(PETSC_DEFAULT), type(SOR_FORWARD_SWEEP)
+{
+ auto val = option.get_optional<double>("omega");
+ omega = *val;
+
+ auto n_its = option.get_optional<int>("local_iterations");
+ lits = *n_its;
+
+ n_its = option.get_optional<int>("parallel_iterations");
+ its = *n_its;
+
+ auto type_name = option.get_optional<std::string>("type");
+ if(type_name->find("pc_sor_symmetric") != std::string::npos)
+ {
+ type = SOR_FORWARD_SWEEP;
+ }
+ if(type_name->find("pc_sor_backward") != std::string::npos)
+ {
+ type = SOR_BACKWARD_SWEEP;
+ }
+ if(type_name->find("pc_sor_local_forward") != std::string::npos)
+ {
+ type = SOR_LOCAL_FORWARD_SWEEP;
+ }
+ if(type_name->find("pc_sor_local_symmetric") != std::string::npos)
+ {
+ type = SOR_LOCAL_SYMMETRIC_SWEEP;
+ }
+ if(type_name->find("pc_sor_local_backward") != std::string::npos)
+ {
+ type = SOR_LOCAL_BACKWARD_SWEEP;
+ }
+}
+
+void PETScPC_SOR_Option::setOption(PC &pc)
+{
+ PCSORSetOmega(pc, omega);
+ PCSORSetIterations(pc, its, lits);
+ PCSORSetSymmetric(pc, type);
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.h b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.h
new file mode 100644
index 00000000000..39a6e86f686
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PC_Option/PETScPC_SOR_Option.h
@@ -0,0 +1,66 @@
+/*!
+ \file PETScPC_SOR_Option.h
+ \brief Define the configuration data for PETSc SOR/SSOR preconditioner.
+
+ \author Wenqing Wang
+ \date 04-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCPC_SOR_OPTION_H_
+#define PETSCPC_SOR_OPTION_H_
+
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+/*!
+ PETSc SOR//SSOR preconditioner options
+
+*/
+struct PETScPC_SOR_Option
+{
+ PETScPC_SOR_Option(const boost::property_tree::ptree &option);
+
+ /// Set SOR/SSOR option
+ void setOption(PC &pc);
+
+ /// Overloaded assign operator
+ void operator = (const PETScPC_SOR_Option& opt)
+ {
+ omega = opt.omega;
+ its = opt.its;
+ lits = opt.lits;
+ type = opt.type;
+ }
+
+ /// Relaxation number
+ PetscReal omega;
+
+ /// Number of parelllel iterations, each parallel iteration
+ /// has 'lits' local iterations
+ PetscInt its;
+
+ /// Number of local iterations
+ PetscInt lits;
+
+ /*!
+ SOR type:
+ SOR_FORWARD_SWEEP
+ SOR_BACKWARD_SWEEP
+ SOR_SYMMETRIC_SWEEP
+ SOR_LOCAL_FORWARD_SWEEP
+ SOR_LOCAL_BACKWARD_SWEEP
+ SOR_LOCAL_SYMMETRIC_SWEEP
+ */
+ MatSORType type;
+};
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
new file mode 100644
index 00000000000..253def0737a
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.cpp
@@ -0,0 +1,183 @@
+/*!
+ \file PETScLinearSolver.cpp
+ \brief Definition of class PETScLinearSolver, which defines a solver object
+ based on PETSc routines.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Mar 2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScLinearSolver.h"
+
+namespace MathLib
+{
+using boost::property_tree::ptree;
+
+PETScLinearSolver::PETScLinearSolver(PETScMatrix &A,
+ const boost::property_tree::ptree &option)
+{
+ KSPCreate(PETSC_COMM_WORLD, &_solver);
+ KSPSetOperators(_solver, A.getRawMatrix(), A.getRawMatrix(), DIFFERENT_NONZERO_PATTERN);
+
+ boost::optional<const ptree&> pt_solver = option.get_child_optional("linear_solver");
+ if(!pt_solver)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"\n*** PETSc linear solver is not specified, bcgs + bjacobi are used.\n");
+ // Base configuration
+ PETScLinearSolverOption opt;
+ opt.setOption(_solver, _pc);
+ KSPSetFromOptions(_solver); // set running time option
+ return;
+ }
+
+ // Preconditioners:
+ boost::optional<const ptree&> pt_pc = option.get_child_optional("preconditioner");
+ if(!pt_pc)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"\n*** PETSc preconditioner is not specified, bjacobi is used.");
+ // Base configuration
+ PETScLinearSolverOption opt(*pt_solver);
+ opt.setOption(_solver, _pc);
+ KSPSetFromOptions(_solver); // set running time option
+ return;
+ }
+
+ // Base configuration
+ PETScLinearSolverOption opt(*pt_solver, *pt_pc);
+
+ opt.setOption(_solver, _pc);
+
+ //----------------------------------------------------------------------
+ // Specific configuration, solver
+ boost::optional<const ptree&> pt_solver_spec = pt_solver->get_child_optional("Richards");
+
+
+ if(pt_solver_spec)
+ {
+ PETScPC_KSP_Richards_Option ksp_opt(*pt_solver_spec);
+ setKSP_Option(ksp_opt);
+ }
+
+ pt_solver_spec = pt_solver->get_child_optional("Chebyshev");
+ if(pt_solver_spec)
+ {
+ PETScPC_KSP_Chebyshev_Option ksp_opt(*pt_solver_spec);
+ setKSP_Option(ksp_opt);
+ }
+
+ pt_solver_spec = pt_solver->get_child_optional("gmres");
+ if(pt_solver_spec)
+ {
+ PETScPC_KSP_GMRES_Option ksp_opt(*pt_solver_spec);
+ setKSP_Option(ksp_opt);
+ }
+
+ //----------------------------------------------------------------------
+ // Specific configuration, preconditioner
+ boost::optional<const ptree&> pt_pc_spec = pt_pc->get_child_optional("sor");
+ if(pt_pc_spec)
+ {
+ PETScPC_SOR_Option pc_opt(*pt_pc_spec);
+ setPC_Option(pc_opt);
+ }
+
+ pt_pc_spec = pt_pc->get_child_optional("asm");
+ if(pt_pc_spec)
+ {
+ PETScPC_ASM_Option pc_opt(*pt_pc_spec);
+ setPC_Option(pc_opt);
+ }
+
+ pt_pc_spec = pt_pc->get_child_optional("amg");
+ if(pt_pc_spec)
+ {
+ PETScPC_AMG_Option pc_opt(*pt_pc_spec);
+ setPC_Option(pc_opt);
+ }
+
+ //----------------------------------------------------------------------------
+ // Only for sub-block preconditioner or sequential computation if it is needed.
+ // ILU or ICC
+ pt_pc_spec = pt_pc->get_child_optional("ilu");
+ if(pt_pc_spec)
+ {
+ PETScPC_ILU_Option pc_opt(*pt_pc_spec);
+ setPC_Option(pc_opt);
+ }
+ else
+ {
+ pt_pc_spec = pt_pc->get_child_optional("icc");
+ if(pt_pc_spec)
+ {
+ PETScPC_ILU_Option pc_opt(*pt_pc_spec);
+ setPC_Option(pc_opt);
+ }
+ }
+
+ pt_pc_spec = pt_pc->get_child_optional("lu");
+ if(pt_pc_spec)
+ {
+ PETScPC_LU_Option pc_opt(*pt_pc_spec);
+ setPC_Option(pc_opt);
+ }
+ //----------------------------------------------------------------------------
+
+ //
+ KSPSetFromOptions(_solver); // set running time option
+}
+
+void PETScLinearSolver::solve(const PETScVector &b, PETScVector &x)
+{
+// define TEST_MEM_PETSC
+#ifdef TEST_MEM_PETSC
+ PetscLogDouble mem1, mem2;
+ PetscMemoryGetCurrentUsage(&mem1);
+#endif
+
+ KSPSolve(_solver, b.getData(), x.getData());
+
+ KSPConvergedReason reason;
+ KSPGetConvergedReason(_solver, &reason);
+
+ if(reason == KSP_DIVERGED_INDEFINITE_PC)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"\nDivergence because of indefinite preconditioner;\n");
+ PetscPrintf(PETSC_COMM_WORLD,"Run the executable again but with -pc_factor_shift_positive_definite option.\n");
+ }
+ else if(reason < 0)
+ {
+ PetscPrintf(PETSC_COMM_WORLD,"\nOther kind of divergence: this should not happen.\n");
+ }
+ else
+ {
+ const char *slv_type;
+ const char *prc_type;
+ KSPGetType(_solver, &slv_type);
+ PCGetType(_pc, &prc_type);
+
+ PetscPrintf(PETSC_COMM_WORLD,"\n================================================");
+ PetscPrintf(PETSC_COMM_WORLD, "\nLinear solver %s with %s preconditioner",
+ slv_type, prc_type);
+
+ int its;
+ KSPGetIterationNumber(_solver, &its);
+ PetscPrintf(PETSC_COMM_WORLD,"\nConvergence in %d iterations.\n", its);
+ PetscPrintf(PETSC_COMM_WORLD,"\n================================================");
+ }
+ PetscPrintf(PETSC_COMM_WORLD,"\n");
+
+#ifdef TEST_MEM_PETSC
+ PetscMemoryGetCurrentUsage(&mem2);
+ PetscPrintf(PETSC_COMM_WORLD, "###Memory usage by solver. Before :%f After:%f Increase:%d\n", mem1, mem2, (int)(mem2 - mem1));
+#endif
+}
+
+} //end of namespace
+
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolver.h b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
new file mode 100644
index 00000000000..d0082a9e9bc
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolver.h
@@ -0,0 +1,83 @@
+/*!
+ \file PETScLinearSolver.h
+ \brief Declaration of class PETScLinearSolver, which defines a solver object
+ based on PETSc routines.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Sep 2013
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#ifndef PETSCLINEARSOLVER_H_
+#define PETSCLINEARSOLVER_H_
+
+#include "PETScMatrix.h"
+#include "PETScVector.h"
+
+#include "PETScLinearSolverOption.h"
+
+#include "KSP_Option/PETScPC_KSP_Chebyshev_Option.h"
+#include "KSP_Option/PETScPC_KSP_Richards_Option.h"
+#include "KSP_Option/PETScPC_KSP_GMRES_Option.h"
+
+#include "PC_Option/PETScPC_ILU_Option.h"
+#include "PC_Option/PETScPC_SOR_Option.h"
+#include "PC_Option/PETScPC_LU_Option.h"
+#include "PC_Option/PETScPC_ASM_Option.h"
+#include "PC_Option/PETScPC_AMG_Option.h"
+
+namespace MathLib
+{
+
+/*!
+ A class of linear solver based on PETSc rountines.
+
+*/
+class PETScLinearSolver
+{
+ public:
+
+ /*!
+ Constructor.
+ \param A Matrix, cannot be constant.
+ \param solver_name Solver name.
+ \param pc_name Preconditioner name.
+ */
+ PETScLinearSolver(PETScMatrix &A, const boost::property_tree::ptree &option);
+
+ template <typename T_KSP_OPTION> void setKSP_Option(T_KSP_OPTION &ksp_opt)
+ {
+ ksp_opt.setOption(_solver);
+ }
+
+ template <typename T_PC_OPTION> void setPC_Option(T_PC_OPTION &pc_opt)
+ {
+ pc_opt.setOption(_pc);
+ }
+
+ ~PETScLinearSolver()
+ {
+ KSPDestroy(&_solver);
+ }
+
+ /*!
+ Solve a system of equations.
+ \param b The right hand of the equations.
+ \param x The solutions to be solved.
+ */
+ void solve(const PETScVector &b, PETScVector &x);
+
+ private:
+ KSP _solver; ///< Slover type.
+ PC _pc; ///< Preconditioner type.
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolverOption.cpp b/MathLib/LinAlg/PETSc/PETScLinearSolverOption.cpp
new file mode 100644
index 00000000000..645f3b91d4d
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolverOption.cpp
@@ -0,0 +1,114 @@
+/*!
+ \file PETScLinearSolverOption.cpp
+ \brief Define members of PETScLinearSolverOption
+
+ \author Wenqing Wang
+ \date 02-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScLinearSolverOption.h"
+
+namespace MathLib
+{
+
+using boost::property_tree::ptree;
+
+PETScLinearSolverOption::
+PETScLinearSolverOption(const boost::property_tree::ptree &ksp_option)
+ : _solver_name("bcgs"), _pc_name("bjacobi"), _preco_side(PC_LEFT),
+ _max_it(2000), _rtol(1.e-5), _atol(PETSC_DEFAULT), _dtol(PETSC_DEFAULT)
+{
+ auto solver_type = ksp_option.get_optional<std::string>("solver_type");
+ if(solver_type)
+ {
+ _solver_name = *solver_type;
+ }
+
+ auto max_iteration_step = ksp_option.get_optional<int>("max_it");
+ if(max_iteration_step)
+ {
+ _max_it = *max_iteration_step;
+ }
+
+ auto error_tolerance = ksp_option.get_optional<double>("rtol");
+ if(error_tolerance)
+ {
+ _rtol = *error_tolerance;
+ }
+
+ error_tolerance = ksp_option.get_optional<double>("atol");
+ if(error_tolerance)
+ {
+ _atol = *error_tolerance;
+ }
+
+ error_tolerance = ksp_option.get_optional<double>("dtol");
+ if(error_tolerance)
+ {
+ _dtol = *error_tolerance;
+ }
+}
+
+PETScLinearSolverOption::
+PETScLinearSolverOption(const boost::property_tree::ptree &ksp_option,
+ const boost::property_tree::ptree &pc_option)
+ : _solver_name("bcgs"), _pc_name("bjacobi"), _preco_side(PC_LEFT),
+ _max_it(2000), _rtol(1.e-5), _atol(PETSC_DEFAULT), _dtol(PETSC_DEFAULT)
+{
+ auto solver_type = ksp_option.get_optional<std::string>("solver_type");
+ if(solver_type)
+ {
+ _solver_name = *solver_type;
+ }
+
+ auto max_iteration_step = ksp_option.get_optional<int>("max_it");
+ if(max_iteration_step)
+ {
+ _max_it = *max_iteration_step;
+ }
+
+ auto error_tolerance = ksp_option.get_optional<double>("rtol");
+ if(error_tolerance)
+ {
+ _rtol = *error_tolerance;
+ }
+
+ error_tolerance = ksp_option.get_optional<double>("atol");
+ if(error_tolerance)
+ {
+ _atol = *error_tolerance;
+ }
+
+ error_tolerance = ksp_option.get_optional<double>("dtol");
+ if(error_tolerance)
+ {
+ _dtol = *error_tolerance;
+ }
+
+ // Preconditioners:
+ auto pc_type = pc_option.get_optional<std::string>("pc_type");
+ if(pc_type)
+ {
+ _pc_name = *pc_type;
+ }
+
+ auto pc_side = pc_option.get_optional<std::string>("pc_side");
+ if(pc_side)
+ {
+ if(pc_side->find("left") != std::string::npos)
+ _preco_side = PC_LEFT;
+ if(pc_side->find("right") != std::string::npos)
+ _preco_side = PC_RIGHT;
+ if(pc_side->find("symmetric") != std::string::npos)
+ _preco_side = PC_SYMMETRIC;
+ }
+}
+
+} // end namespace
+
diff --git a/MathLib/LinAlg/PETSc/PETScLinearSolverOption.h b/MathLib/LinAlg/PETSc/PETScLinearSolverOption.h
new file mode 100644
index 00000000000..2a9b78568e9
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScLinearSolverOption.h
@@ -0,0 +1,161 @@
+/*!
+ \file PETScLinearSolverOption.h
+ \brief Define the configuration data for the PETSc linear solver.
+
+ \author Wenqing Wang
+ \date 02-2014
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+#ifndef PETSCLINEARSOLVEROPTION_H_
+#define PETSCLINEARSOLVEROPTION_H_
+
+#include <string>
+#include <boost/property_tree/ptree.hpp>
+
+#include <petscksp.h>
+
+namespace MathLib
+{
+/*!
+ \brief Class to manage the basic configuration data for
+ a Krylov subspace (KSP) interation method and a procondtioner (PC)
+ of PETSc routines.
+*/
+class PETScLinearSolverOption
+{
+ public:
+
+ /*!
+ Default constructor
+ */
+ PETScLinearSolverOption() : _solver_name("bcgs"), _pc_name("bjacobi"),
+ _preco_side(PC_LEFT), _max_it(2000),
+ _rtol(1.e-5), _atol(PETSC_DEFAULT), _dtol(PETSC_DEFAULT)
+ {}
+
+ /*!
+ \param ksp_option ptree of basic configuration data for solver
+ */
+ PETScLinearSolverOption(const boost::property_tree::ptree &ksp_option);
+
+ /*!
+ \param ksp_option ptree of basic configuration data for solver
+ \param pc_option ptree of basic configuration data for preconditioner
+ */
+ PETScLinearSolverOption(const boost::property_tree::ptree &ksp_option,
+ const boost::property_tree::ptree &pc_option);
+
+ /// Set basic options for a KSP and a PC
+ void setOption(KSP &ksp, PC &pc)
+ {
+ KSPSetType(ksp, _solver_name.c_str());
+ KSPGetPC(ksp, &pc);
+ PCSetType(pc, _pc_name.c_str());
+ KSPSetPCSide(ksp, _preco_side);
+
+ KSPSetTolerances(ksp, _rtol, _atol, _dtol, _max_it);
+ }
+
+ private:
+ /*!
+ The name of solver, and it could be one of the following names
+ "richardson"
+ "chebychev"
+ "cg"
+ "cgne"
+ "nash"
+ "stcg"
+ "gltr"
+ "gmres"
+ "fgmres"
+ "lgmres"
+ "dgmres"
+ "tcqmr"
+ "bcgs"
+ "ibcgs"
+ "bcgsl"
+ "cgs"
+ "tfqmr"
+ "cr"
+ "lsqr"
+ "preonly"
+ "qcg"
+ "bicg"
+ "minres"
+ "symmlq"
+ "lcd"
+ "python"
+ "broyden"
+ "gcr"
+ "ngmres"
+ "specest"
+ */
+ std::string _solver_name;
+
+ /*!
+ The name of preconditioner, and it could be one of the following names
+ none
+ jacobi
+ sor
+ lu
+ shell
+ bjacobi
+ mg
+ eisenstat
+ ilu
+ icc
+ asm
+ gasm
+ ksp
+ composite
+ redundant
+ spai
+ nn
+ cholesky
+ pbjacobi
+ mat
+ hypre
+ parms
+ fieldsplit
+ tfs
+ ml
+ prometheus
+ galerkin
+ exotic
+ hmpi
+ supportgraph
+ asa
+ cp
+ bfbt
+ lsc
+ python
+ pfmg
+ syspfmg
+ redistribute
+ sacusp
+ sacusppoly
+ bicgstabcusp
+ svd
+ ainvcusp
+ gamg
+ */
+ std::string _pc_name;
+
+ /// Flag for which side preconditioning.
+ PCSide _preco_side;
+
+ PetscInt _max_it; ///< Maximum iteration.
+
+ PetscReal _rtol; ///< Tolerance for the relative error, \f$e=|r|/|b|\f$.
+ PetscReal _atol; ///< Tolerance for the absolute error, \f$e=|r|\f$.
+ PetscReal _dtol; ///< Relative increase in the residual.
+};
+
+} // end namespace
+#endif
+
diff --git a/MathLib/LinAlg/PETSc/PETScMatrix.cpp b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
index 2cafefa0445..8e706128cb0 100644
--- a/MathLib/LinAlg/PETSc/PETScMatrix.cpp
+++ b/MathLib/LinAlg/PETSc/PETScMatrix.cpp
@@ -75,7 +75,7 @@ void PETScMatrix::viewer(const std::string &file_name, const PetscViewerFormat v
MatView(_A,viewer);
// This preprocessor is only for debugging, e.g. dump the matrix and exit the program.
-//#define EXIT_TEST
+//#define EXIT_TEST
#ifdef EXIT_TEST
MatDestroy(&_A);
PetscFinalize();
@@ -91,7 +91,7 @@ void PETScMatrix::create(const PetscInt d_nz, const PetscInt o_nz)
MatSetFromOptions(_A);
- // for a dense matrix: MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
+ MatSeqAIJSetPreallocation(_A, d_nz, PETSC_NULL);
MatMPIAIJSetPreallocation(_A, d_nz, PETSC_NULL, o_nz, PETSC_NULL);
MatGetOwnershipRange(_A, &_start_rank, &_end_rank);
diff --git a/MathLib/LinAlg/PETSc/PETScTools.cpp b/MathLib/LinAlg/PETSc/PETScTools.cpp
new file mode 100644
index 00000000000..3d9e0e46213
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScTools.cpp
@@ -0,0 +1,41 @@
+/*!
+ \file PETScTools.cpp
+ \brief Definition of a function related to PETSc solver interface to assign
+ the Dirichlet boundary conditions.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Sep 2013
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#include "PETScTools.h"
+
+namespace MathLib
+{
+
+void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
+ const std::vector<PetscInt> &vec_knownX_id,
+ const std::vector<PetscScalar> &vec_knownX_x)
+{
+ A.setRowsColumnsZero(vec_knownX_id);
+ A.finalizeAssembly();
+
+ if(vec_knownX_id.size() > 0)
+ {
+ x.set(vec_knownX_id, vec_knownX_x);
+ b.set(vec_knownX_id, vec_knownX_x);
+ }
+
+ x.finalizeAssembly();
+ b.finalizeAssembly();
+}
+
+} // end of namespace MathLib
+
+
diff --git a/MathLib/LinAlg/PETSc/PETScTools.h b/MathLib/LinAlg/PETSc/PETScTools.h
new file mode 100644
index 00000000000..bad353b9102
--- /dev/null
+++ b/MathLib/LinAlg/PETSc/PETScTools.h
@@ -0,0 +1,45 @@
+/*!
+ \file PETScTools.h
+ \brief Declaration of a function related to PETSc solver interface to assign
+ the Dirichlet boundary conditions.
+
+ \author Wenqing Wang
+ \version
+ \date Nov 2011 - Sep 2013
+
+ \copyright
+ Copyright (c) 2013, OpenGeoSys Community (http://www.opengeosys.org)
+ Distributed under a Modified BSD License.
+ See accompanying file LICENSE.txt or
+ http://www.opengeosys.org/project/license
+*/
+
+#ifndef PETSCTOOLS_H_
+#define PETSCTOOLS_H_
+
+#include <vector>
+
+#include "PETScMatrix.h"
+#include "PETScVector.h"
+
+namespace MathLib
+{
+/*!
+ \brief apply known solutions to a system of linear equations
+
+ This function introduces the constants into the system by the penalty method.
+
+ \param A Coefficient matrix
+ \param b RHS vector
+ \param vec_knownX_id a vector of known solution entry IDs
+ \param vec_knownX_x a vector of known solutions
+ \param penalty_scaling value for scaling some matrix and right hand side
+ entries to enforce some conditions
+*/
+void applyKnownSolution(PETScMatrix &A, PETScVector &b, PETScVector &x,
+ const std::vector<PetscInt> &vec_knownX_id,
+ const std::vector<PetscScalar> &vec_knownX_x);
+} // end of namespace MathLib
+
+#endif //end of PETSCTOOLS_H_
+
diff --git a/Tests/MathLib/TestLinearSolver.cpp b/Tests/MathLib/TestLinearSolver.cpp
index 6c46e62fe38..cd1d31fdd99 100644
--- a/Tests/MathLib/TestLinearSolver.cpp
+++ b/Tests/MathLib/TestLinearSolver.cpp
@@ -1,7 +1,8 @@
/**
* \file
* \author Norihiro Watanabe
- * \date 2013-04-16
+ * \author Wenqing Wang
+ * \date 2013-04-16, 2014-04
* \brief Implementation tests.
*
* \copyright
@@ -21,11 +22,19 @@
#include "MathLib/LinAlg/Dense/DenseTools.h"
#include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
#include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
+
#ifdef USE_LIS
#include "MathLib/LinAlg/Lis/LisLinearSolver.h"
#include "MathLib/LinAlg/Lis/LisTools.h"
#endif
+#ifdef USE_PETSC
+#include "MathLib/LinAlg/PETSc/PETScMatrix.h"
+#include "MathLib/LinAlg/PETSc/PETScVector.h"
+#include "MathLib/LinAlg/PETSc/PETScLinearSolver.h"
+#include "MathLib/LinAlg/PETSc/PETScTools.h"
+#endif
+
#include "../TestTools.h"
namespace
@@ -34,7 +43,8 @@ namespace
template<class T_Mat>
void setMatrix9x9(T_Mat &mat)
{
- double d_mat[] = {
+ double d_mat[] =
+ {
6.66667e-012, -1.66667e-012, 0, -1.66667e-012, -3.33333e-012, 0, 0, 0, 0,
-1.66667e-012, 1.33333e-011, -1.66667e-012, -3.33333e-012, -3.33333e-012, -3.33333e-012, 0, 0, 0,
0, -1.66667e-012, 6.66667e-012, 0, -3.33333e-012, -1.66667e-012, 0, 0, 0,
@@ -45,9 +55,9 @@ void setMatrix9x9(T_Mat &mat)
0, 0, 0, -3.33333e-012, -3.33333e-012, -3.33333e-012, -1.66667e-012, 1.33333e-011, -1.66667e-012,
0, 0, 0, 0, -3.33333e-012, -1.66667e-012, 0, -1.66667e-012, 6.66667e-012
};
- for (unsigned i = 0; i < 9; i++)
- for (unsigned j = 0; j < 9; j++)
- mat.setValue(i, j, d_mat[i*9+j]);
+ for (unsigned i = 0; i < 9; i++)
+ for (unsigned j = 0; j < 9; j++)
+ mat.setValue(i, j, d_mat[i*9+j]);
}
struct Example1
@@ -59,7 +69,7 @@ struct Example1
double* exH;
Example1()
- : mat(dim_eqs, dim_eqs), exH(new double[dim_eqs])
+ : mat(dim_eqs, dim_eqs), exH(new double[dim_eqs])
{
setMatrix9x9(mat);
std::size_t int_dirichlet_bc_id[] = {2,5,8,0,3,6};
@@ -67,7 +77,8 @@ struct Example1
vec_dirichlet_bc_value.resize(6);
std::fill(vec_dirichlet_bc_value.begin(), vec_dirichlet_bc_value.begin()+3, .0);
std::fill(vec_dirichlet_bc_value.begin()+3, vec_dirichlet_bc_value.end(), 1.0);
- for (std::size_t i=0; i<9; i++) {
+ for (std::size_t i=0; i<9; i++)
+ {
if (i%3==0) exH[i] = 1.0;
if (i%3==1) exH[i] = 0.5;
if (i%3==2) exH[i] = 0.;
@@ -87,8 +98,10 @@ void checkLinearSolverInterface(T_MATRIX &A, boost::property_tree::ptree &ls_opt
// set a coefficient matrix
A.setZero();
- for (size_t i=0; i<ex1.dim_eqs; i++) {
- for (size_t j=0; j<ex1.dim_eqs; j++) {
+ for (size_t i=0; i<ex1.dim_eqs; i++)
+ {
+ for (size_t j=0; j<ex1.dim_eqs; j++)
+ {
double v = ex1.mat(i, j);
if (v!=.0)
A.add(i, j, v);
@@ -112,6 +125,69 @@ void checkLinearSolverInterface(T_MATRIX &A, boost::property_tree::ptree &ls_opt
}
+#ifdef USE_PETSC
+template <class T_MATRIX, class T_VECTOR, class T_LINEAR_SOVLER>
+void checkLinearSolverInterface(T_MATRIX &A, T_VECTOR &b, boost::property_tree::ptree &ls_option)
+{
+ int mrank;
+ MPI_Comm_rank(PETSC_COMM_WORLD, &mrank);
+ // Add entries
+ MathLib::DenseMatrix<double> loc_m(2, 2);
+ loc_m(0, 0) = 1. + mrank;
+ loc_m(0, 1) = 2. + mrank;
+ loc_m(1, 0) = 3. + mrank;
+ loc_m(1, 1) = 4. + mrank;
+
+ std::vector<int> row_pos(2);
+ std::vector<int> col_pos(2);
+ row_pos[0] = 2 * mrank;
+ row_pos[1] = 2 * mrank + 1;
+ col_pos[0] = row_pos[0];
+ col_pos[1] = row_pos[1];
+
+ A.add(row_pos, col_pos, loc_m);
+
+ MathLib::finalizeMatrixAssembly(A);
+
+ const bool deep_copy = false;
+ T_VECTOR x(b, deep_copy);
+
+ std::vector<double> local_vec(2);
+ local_vec[0] = mrank+1;
+ local_vec[1] = 2. * (mrank+1);
+ x.set(col_pos, local_vec);
+
+ double x0[6];
+ double x1[6];
+ x.getGlobalVector(x0);
+
+ A.multiply(x, b);
+
+ // apply BC
+ std::vector<int> bc_id; /// Type must be int to match Petsc_Int
+ std::vector<double> bc_value;
+
+ if(mrank == 1)
+ {
+ bc_id.resize(1);
+ bc_value.resize(1);
+ bc_id[0] = 2 * mrank;
+ bc_value[0] = mrank+1;
+ }
+
+ MathLib::applyKnownSolution(A, b, x, bc_id, bc_value);
+
+ MathLib::finalizeMatrixAssembly(A);
+
+ // solve
+ T_LINEAR_SOVLER ls(A, ls_option);
+ ls.solve(b, x);
+ x.getGlobalVector(x1);
+
+ ASSERT_ARRAY_NEAR(x0, x1, 6, 1e-5);
+}
+#endif
+
} // end namespace
TEST(MathLib, CheckInterface_GaussAlgorithm)
@@ -142,3 +218,201 @@ TEST(Math, CheckInterface_Lis)
}
#endif
+#ifdef USE_PETSC
+TEST(Math, CheckInterface_PETSc_Linear_Solver_default)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ // set solver options using Boost property tree
+ boost::property_tree::ptree t_root;
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, t_root);
+}
+
+TEST(Math, CheckInterface_PETSc_Linear_Solver_basic)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ // set solver options using Boost property tree
+ boost::property_tree::ptree t_root;
+
+ boost::property_tree::ptree t_solver;
+ t_root.put("solver_package", "PETSc");
+ t_solver.put("solver_type", "bcgs");
+ t_solver.put("rtol", 1e-8);
+ t_solver.put("atol", 1e-50);
+ t_solver.put("max_it", 1000);
+ t_solver.put("pc_side", "left");
+
+ t_root.put_child("linear_solver", t_solver);
+
+ boost::property_tree::ptree t_pc;
+ t_pc.put("pc_type", "bjacobi");
+ t_root.put_child("preconditioner", t_pc);
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, t_root);
+}
+
+TEST(Math, CheckInterface_PETSc_Linear_Solver_chebyshev_sor)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ // set solver options using Boost property tree
+ boost::property_tree::ptree t_root;
+
+ boost::property_tree::ptree t_solver;
+ t_root.put("solver_package", "PETSc");
+ t_solver.put("solver_type", "chebyshev");
+ t_solver.put("rtol", 1e-8);
+ t_solver.put("atol", 1e-50);
+ t_solver.put("max_it", 1000);
+
+ // If ommit the following specific configurations,
+ // the default setting will be taken.
+ boost::property_tree::ptree t_solver_spec;
+ t_solver_spec.put("smallest_eignvalue", 0.1);
+ t_solver_spec.put("maximum_eignvalue", 98.);
+ t_solver.put_child("chebyshev", t_solver_spec);
+
+ t_root.put_child("linear_solver", t_solver);
+
+ // pc
+ boost::property_tree::ptree t_pc;
+ t_pc.put("pc_type", "sor");
+
+//#define SEQ_TEST
+#ifdef SEQ_TEST //only for sequential SOR
+ // If ommit the following specific configurations,
+ // the default setting will be taken.
+ boost::property_tree::ptree t_pc_spec;
+ t_pc_spec.put("omega", 0.8);
+ t_pc_spec.put("local_iterations", 5);
+ t_pc_spec.put("parallel_iterations", 2);
+ t_pc_spec.put("type", "pc_sor_backward");
+ t_pc.put_child("sor", t_pc_spec);
+#endif
+
+ t_root.put_child("preconditioner", t_pc);
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, t_root);
+}
+
+TEST(Math, CheckInterface_PETSc_Linear_Solver_gmres_amg)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ // set solver options using Boost property tree
+ boost::property_tree::ptree t_root;
+
+ boost::property_tree::ptree t_solver;
+ t_root.put("solver_package", "PETSc");
+ t_solver.put("solver_type", "gmres");
+ t_solver.put("rtol", 1e-8);
+ t_solver.put("atol", 1e-50);
+ t_solver.put("max_it", 1000);
+
+ // If ommit the following specific configurations,
+ // the default setting will be taken.
+ boost::property_tree::ptree t_solver_spec;
+ t_solver_spec.put("restart_number", 20);
+ t_solver_spec.put("is_modified_ram_schmidt_orthog", false);
+ t_solver_spec.put("refine_type", 1);
+ t_solver.put_child("gmres", t_solver_spec);
+
+ t_root.put_child("linear_solver", t_solver);
+
+ // pc
+ boost::property_tree::ptree t_pc;
+ t_pc.put("pc_type", "gamg");
+
+ // If ommit the following specific configurations,
+ // the default setting will be taken.
+ boost::property_tree::ptree t_pc_spec;
+ t_pc_spec.put("agg_nsmooths", 2.0);
+ t_pc_spec.put("type", "agg");
+ t_pc.put_child("gamg", t_pc_spec);
+
+ t_root.put_child("preconditioner", t_pc);
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, t_root);
+}
+
+TEST(Math, CheckInterface_PETSc_Linear_Solver_cg_asm)
+{
+ MathLib::PETScMatrixOption opt;
+ opt.d_nz = 2;
+ opt.o_nz = 0;
+ opt.is_global_size = false;
+ opt.n_local_cols = 2;
+ MathLib::PETScMatrix A(2, opt);
+
+ const bool is_gloabal_size = false;
+ MathLib::PETScVector b(2, is_gloabal_size);
+
+ // set solver options using Boost property tree
+ boost::property_tree::ptree t_root;
+
+ boost::property_tree::ptree t_solver;
+ t_root.put("solver_package", "PETSc");
+ t_solver.put("solver_type", "cg");
+ t_solver.put("rtol", 1e-8);
+ t_solver.put("atol", 1e-50);
+ t_solver.put("max_it", 1000);
+
+ t_root.put_child("linear_solver", t_solver);
+
+ // pc
+ boost::property_tree::ptree t_pc;
+ t_pc.put("pc_type", "gasm");
+
+ // If ommit the following specific configurations,
+ // the default setting will be taken.
+ boost::property_tree::ptree t_pc_spec;
+ t_pc_spec.put("type", "interpolate");
+ t_pc.put_child("asm", t_pc_spec);
+
+ t_root.put_child("preconditioner", t_pc);
+
+ checkLinearSolverInterface<MathLib::PETScMatrix, MathLib::PETScVector,
+ MathLib::PETScLinearSolver>(A, b, t_root);
+}
+
+#endif
+
+
--
GitLab