Skip to content
Snippets Groups Projects
Commit 8c7baf0b authored by wenqing's avatar wenqing
Browse files

Set PETSc solve use command line option only according to the suggestion by Nori

parent 3d34af8b
No related branches found
No related tags found
No related merge requests found
Showing
with 17 additions and 1124 deletions
......@@ -31,10 +31,6 @@ 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)
......
/*!
\file PETScKSP_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 "PETScKSP_Chebyshev_Option.h"
namespace MathLib
{
using boost::property_tree::ptree;
PETScKSP_Chebyshev_Option::
PETScKSP_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
/*!
\file PETScKSP_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 PETSCKSP_CHEBYSHEV_OPTION_H_
#define PETSCKSP_CHEBYSHEV_OPTION_H_
#include <boost/property_tree/ptree.hpp>
#include <petscksp.h>
namespace MathLib
{
/*!
PETSc KSP Chebyshev options
*/
struct PETScKSP_Chebyshev_Option
{
PETScKSP_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 PETScKSP_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
/*!
\file PETScKSP_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 "PETScKSP_GMRES_Option.h"
namespace MathLib
{
using boost::property_tree::ptree;
PETScKSP_GMRES_Option::
PETScKSP_GMRES_Option(const boost::property_tree::ptree &option)
: restart_number(30), is_modified_gram_schmidt(false),
refine_type(KSP_GMRES_CGS_REFINE_NEVER)
{
auto val = option.get_optional<double>("restart_number");
restart_number = *val;
boost::optional<bool> bool_vals = option.get_optional<bool>("modified_gram_schmidt");
is_modified_gram_schmidt = *bool_vals;
auto refine_type = option.get_optional<int>("refine_type");
switch(*refine_type)
{
case 0:
refine_type = KSP_GMRES_CGS_REFINE_NEVER;
break;
case 1:
refine_type = KSP_GMRES_CGS_REFINE_IFNEEDED;
break;
case 2:
refine_type = KSP_GMRES_CGS_REFINE_ALWAYS;
break;
default:
refine_type = KSP_GMRES_CGS_REFINE_NEVER;
break;
}
}
/// Set Chebyshev option
void PETScKSP_GMRES_Option::setOption(KSP &solver)
{
KSPGMRESSetRestart(solver, restart_number);
if(is_modified_gram_schmidt)
{
KSPGMRESSetOrthogonalization(solver, KSPGMRESClassicalGramSchmidtOrthogonalization);
}
KSPGMRESSetCGSRefinementType(solver, refine_type);
}
} // end namespace
/*!
\file PETScKSP_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 PETSCKSP_GMRES_OPTION_H_
#define PETSCKSP_GMRES_OPTION_H_
#include <boost/property_tree/ptree.hpp>
#include <petscksp.h>
namespace MathLib
{
/*!
PETSc KSP GMRES options
*/
struct PETScKSP_GMRES_Option
{
PETScKSP_GMRES_Option(const boost::property_tree::ptree &option);
/// Set GMRES option
void setOption(KSP &solver);
/// Overloaded assign operator
void operator = (const PETScKSP_GMRES_Option& opt)
{
restart_number = opt.restart_number;
is_modified_gram_schmidt = opt.is_modified_gram_schmidt;
refine_type = opt.refine_type;
}
/// Restart number of GMRES.
PetscInt restart_number;
/// Flag for the modified Gram-Schmidt orthogonalization.
bool is_modified_gram_schmidt;
/*!
\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;
};
} // end namespace
#endif
/*!
\file PETScKSP_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 "PETScKSP_Richards_Option.h"
namespace MathLib
{
using boost::property_tree::ptree;
PETScKSP_Richards_Option::
PETScKSP_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
/*!
\file PETScKSP_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 PETSCKSP_RICHARDS_OPTION_H_
#define PETSCKSP_RICHARDS_OPTION_H_
#include <boost/property_tree/ptree.hpp>
#include <petscksp.h>
namespace MathLib
{
/*!
PETSc KSP Richards options
*/
struct PETScKSP_Richards_Option
{
PETScKSP_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 PETScKSP_Richards_Option& opt)
{
damping_factor_richards = opt.damping_factor_richards;
}
/// Damping factor for Richards.
PetscReal damping_factor_richards;
};
} // end namespace
#endif
/*!
\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
/*!
\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
/*!
\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("none") != std::string::npos)
{
type = PC_ASM_NONE;
}
}
} // end namespace
/*!
\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
/*!
\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
/*!
\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
/*!
\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
/*!
\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
/*!
\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
/*!
\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
......@@ -18,118 +18,19 @@
namespace MathLib
{
using boost::property_tree::ptree;
PETScLinearSolver::PETScLinearSolver(PETScMatrix &A,
const boost::property_tree::ptree &option)
PETScLinearSolver::PETScLinearSolver(PETScMatrix &A, const std::string prefix)
{
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)
{
PETScKSP_Richards_Option ksp_opt(*pt_solver_spec);
setKSP_Option(ksp_opt);
}
pt_solver_spec = pt_solver->get_child_optional("Chebyshev");
if(pt_solver_spec)
{
PETScKSP_Chebyshev_Option ksp_opt(*pt_solver_spec);
setKSP_Option(ksp_opt);
}
pt_solver_spec = pt_solver->get_child_optional("gmres");
if(pt_solver_spec)
{
PETScKSP_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);
}
//----------------------------------------------------------------------------
//
KSPSetOptionsPrefix(_solver, prefix.c_str());
KSPSetFromOptions(_solver); // set running time option
KSPGetPC(_solver, &_pc);
PCSetOptionsPrefix(_pc, prefix.c_str());
PCSetFromOptions(_pc); // set running time option
}
void PETScLinearSolver::solve(const PETScVector &b, PETScVector &x)
......
......@@ -17,20 +17,12 @@
#ifndef PETSCLINEARSOLVER_H_
#define PETSCLINEARSOLVER_H_
#include "PETScMatrix.h"
#include "PETScVector.h"
#include "PETScLinearSolverOption.h"
#include<string>
#include "KSP_Option/PETScKSP_Chebyshev_Option.h"
#include "KSP_Option/PETScKSP_Richards_Option.h"
#include "KSP_Option/PETScKSP_GMRES_Option.h"
#include "petscksp.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"
#include "PETScMatrix.h"
#include "PETScVector.h"
namespace MathLib
{
......@@ -38,6 +30,7 @@ namespace MathLib
/*!
A class of linear solver based on PETSc rountines.
All options for KSP and PC must given in the command line.
*/
class PETScLinearSolver
{
......@@ -45,28 +38,12 @@ class PETScLinearSolver
/*!
Constructor
\param A Matrix, cannot be constant.
\param option Configuration data for solver and preconditioner.
\param A Matrix, cannot be constant.
\param prefix Name used to distinguish the options in the command
line for this solver. It can be the name of the PDE
that owns an instance of this class.
*/
PETScLinearSolver(PETScMatrix &A, const boost::property_tree::ptree &option);
/*!
Set Krylov subspace parameters with given option
\param ksp_opt Given parameters
*/
template <typename T_KSP_OPTION> void setKSP_Option(T_KSP_OPTION &ksp_opt)
{
ksp_opt.setOption(_solver);
}
/*!
Set preconditioner parameters with given option
\param pc_opt Given parameters
*/
template <typename T_PC_OPTION> void setPC_Option(T_PC_OPTION &pc_opt)
{
pc_opt.setOption(_pc);
}
PETScLinearSolver(PETScMatrix &A, const std::string prefix);
~PETScLinearSolver()
{
......
/*!
\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
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment