OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

changed: refactor LinSolParams class

Browse Source 
- change from a struct-of-arrays to an array-of-structs storage scheme
- group ML/GAMG/Hypre settings in their own class
- generalize block syntax. any number of blocks, on any basis,
  is now supported (just on the configuration side for now, still
  hardcoded to 2 blocks on first basis as earlier)
- add unit tests for parsing
This commit is contained in:
Arne Morten Kvarving
2015-10-29 16:53:50 +01:00
committed by Knut Morten Okstad
parent bc0e326880
commit c32274a3c2
14 changed files with 781 additions and 792 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
cmake/Modules/IFEMTesting.cmake
View File

@@ -7,7 +7,7 @@ macro(IFEM_add_test_app path workdir name)
endmacro()
macro(IFEM_add_unittests IFEM_PATH)
IFEM_add_test_app("${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C"
IFEM_add_test_app("${IFEM_PATH}/src/Utility/Test/*.C;${IFEM_PATH}/src/ASM/Test/*.C;${IFEM_PATH}/src/LinAlg/Test/*.C"
${IFEM_PATH}
IFEM
${IFEM_LIBRARIES} ${IFEM_DEPLIBS})

812
src/LinAlg/LinSolParams.C
View File

@@ -21,440 +21,204 @@
#include <iterator>
void LinSolParams::setDefault ()
bool LinSolParams::read (std::istream& is, int nparam)
{
#ifdef HAS_PETSC
nResetSolver = 0;
// Use GMRES with ILU preconditioner as default
method = KSPGMRES;
prec = PCILU;
subprec.resize(1,PCILU);
hypretype.resize(1,"boomeramg");
package.resize(1,MATSOLVERPETSC);
levels.resize(1,0);
mglevels.resize(1,1);
presmoother.resize(1,PCILU);
noPreSmooth.resize(1,1);
postsmoother.resize(1,PCILU);
noPostSmooth.resize(1,1);
noFineSmooth.resize(1,1);
maxCoarseSize.resize(1,1000);
mgKsp.resize(1,"defrichardson");
GAMGtype.resize(1,"agg");
GAMGrepartition.resize(1,PETSC_FALSE);
GAMGuseAggGasm.resize(1,PETSC_FALSE);
overlap.resize(1,1);
nx.resize(1,0);
ny.resize(1,0);
nz.resize(1,0);
nullspc.resize(1,NONE);
nblock = 1;
schur = false;
schurPrec = SIMPLE;
ncomps.resize(2,1);
ncomps.front() = 2;
atol = 1.0e-6;
rtol = 1.0e-6;
dtol = 1.0e3;
maxIts = 1000;
#endif
return false;
}
bool LinSolParams::read (std::istream& is, int nparam)
bool LinSolParams::BlockParams::read(const TiXmlElement* elem)
{
char* cline = 0;
for (int i = 0; i < nparam && (cline = utl::readLine(is)); i++)
utl::getAttribute(elem, "basis", basis);
utl::getAttribute(elem, "components", comps);
const char* value;
const TiXmlElement* child = elem->FirstChildElement();
for (; child; child = child->NextSiblingElement())
if ((value = utl::getValue(child,"pc")))
prec = value;
else if ((value = utl::getValue(child,"package")))
package = value;
else if ((value = utl::getValue(child,"ilu_fill_level")))
ilu_fill_level = atoi(value);
else if ((value = utl::getValue(child,"mglevels")))
mglevel = atoi(value);
else if ((value = utl::getValue(child,"noPreSmooth")))
noPreSmooth = atoi(value);
else if ((value = utl::getValue(child,"noPostSmooth")))
noPostSmooth = atoi(value);
else if ((value = utl::getValue(child,"noFineSmooth")))
noFineSmooth = atoi(value);
else if ((value = utl::getValue(child,"presmoother")))
presmoother = value;
else if ((value = utl::getValue(child,"postsmoother")))
postsmoother = value;
else if ((value = utl::getValue(child,"finesmoother")))
finesmoother = value;
else if ((value = utl::getValue(child,"mgksp")))
mgKSP = value;
else if (!strcasecmp(child->Value(),"dirsmoother")) {
size_t order;
std::string type;
if (!utl::getAttribute(child,"type",type))
return false;
if (!utl::getAttribute(child,"order",order))
return false;
dirSmoother.push_back(DirSmoother(order, type));
}
else if ((value = utl::getValue(child,"overlap")))
overlap = atoi(value);
else if ((value = utl::getValue(child,"nx")))
subdomains[0] = atoi(value);
else if ((value = utl::getValue(child,"ny")))
subdomains[1] = atoi(value);
else if ((value = utl::getValue(child,"nz")))
subdomains[2] = atoi(value);
else if (!strcasecmp(child->Value(),"gamg"))
gamg.read(child);
else if (!strcasecmp(child->Value(),"ml"))
ml.read(child);
else if (!strcasecmp(child->Value(),"hypre"))
hypre.read(child);
else if ((value = utl::getValue(child,"maxCoarseSize")))
maxCoarseSize = atoi(value);
#ifdef HAS_PETSC
if (!strncasecmp(cline,"type",4)) {
char* c = strchr(cline,'=');
for (++c; *c == ' '; c++);
int j = 0;
while (c[j] != EOF && c[j] != '\n' && c[j] != ' ' && c[j] != '\0') j++;
method.assign(c,j);
}
else if (!strncasecmp(cline,"pc",2)) {
char* c = strchr(cline,'=');
for (++c; *c == ' '; c++);
int j = 0;
while (c[j] != EOF && c[j] != '\n' && c[j] != ' ' && c[j] != '\0') j++;
prec.assign(c,j);
}
else if (!strncasecmp(cline,"hypretype",9)) {
char* c = strchr(cline,'=');
for (++c; *c == ' '; c++);
int j = 0;
while (c[j] != EOF && c[j] != '\n' && c[j] != ' ' && c[j] != '\0') j++;
hypretype[0].assign(c,j);
}
else if (!strncasecmp(cline,"package",7)) {
char* c = strchr(cline,'=');
for (++c; *c == ' '; c++);
int j = 0;
while (c[j] != EOF && c[j] != '\n' && c[j] != ' ' && c[j] != '\0') j++;
package[0].assign(c,j);
}
else if (!strncasecmp(cline,"levels",6)) {
char* c = strchr(cline,'=');
levels[0] = atoi(++c);
}
else if (!strncasecmp(cline,"overlap",7)) {
char* c = strchr(cline,'=');
overlap[0] = atoi(++c);
}
else if (!strncasecmp(cline,"atol",4)) {
char* c = strchr(cline,'=');
atol = atof(++c);
}
else if (!strncasecmp(cline,"rtol",4)) {
char* c = strchr(cline,'=');
rtol = atof(++c);
}
else if (!strncasecmp(cline,"dtol",4)) {
char* c = strchr(cline,'=');
dtol = atof(++c);
}
else if (!strncasecmp(cline,"maxits",6)) {
char* c = strchr(cline,'=');
maxIts = atoi(++c);
}
else if (!strncasecmp(cline,"nx",2)) {
char* c = strchr(cline,'=');
nx[0] = atoi(++c);
}
else if (!strncasecmp(cline,"ny",2)) {
char* c = strchr(cline,'=');
ny[0] = atoi(++c);
}
else if (!strncasecmp(cline,"nz",2)) {
char* c = strchr(cline,'=');
nz[0] = atoi(++c);
}
else if (!strncasecmp(cline,"ncomponents",11)) {
char* c = strchr(cline,'=');
std::istringstream this_line(c);
std::istream_iterator<int> begin(this_line), end;
ncomps.assign(begin, end);
nblock = ncomps.size();
}
else if (!strncasecmp(cline,"nullspace",6)) {
char* c = strchr(cline,'=');
if (!strncasecmp(c,"CONSTANT",8))
nullspc[0] = CONSTANT;
else if (!strncasecmp(c,"RIGID_BODY",10))
nullspc[0] = RIGID_BODY;
else if ((value = utl::getValue(child,"nullspace"))) {
if (!strcasecmp(value,"constant"))
nullspace = CONSTANT;
else if (!strcasecmp(value,"rigid_body"))
nullspace = RIGID_BODY;
else
nullspc[0] = NONE;
nullspace = NONE;
}
else {
std::cerr <<" *** LinSolParams::read: Unknown keyword: "
<< cline << std::endl;
return false;
}
#else
;
#endif
else
return false;
return true;
}
bool LinSolParams::read (const TiXmlElement* child)
void LinSolParams::BlockParams::GAMGSettings::read(const TiXmlElement* elem)
{
const char* value;
const TiXmlElement* child = elem->FirstChildElement();
for (; child; child = child->NextSiblingElement())
if ((value = utl::getValue(child,"type")))
type = value;
else if ((value = utl::getValue(child,"repartition")))
repartition = atoi(value);
else if ((value = utl::getValue(child,"use_agg_smoother")))
useAggGasm = atoi(value);
else if ((value = utl::getValue(child,"proc_eq_limit")))
procEqLimit = atoi(value);
else if ((value = utl::getValue(child,"reuse_interpolation")))
reuseInterp = atoi(value);
else if ((value = utl::getValue(child,"threshold")))
threshold = atof(value);
}
void LinSolParams::BlockParams::MLSettings::read(const TiXmlElement* elem)
{
const char* value;
const TiXmlElement* child = elem->FirstChildElement();
for (; child; child = child->NextSiblingElement())
if ((value = utl::getValue(child,"coarse_package")))
coarsePackage = value;
else if ((value = utl::getValue(child,"coarse_solver")))
coarseSolver = value;
else if ((value = utl::getValue(child,"coarsen_scheme")))
coarsenScheme = value;
else if ((value = utl::getValue(child,"symmetrize")))
symmetrize = atoi(value);
else if ((value = utl::getValue(child,"repartition")))
repartition = atoi(value);
else if ((value = utl::getValue(child,"block_scaling")))
blockScaling = atoi(value);
else if ((value = utl::getValue(child,"put_on_single_proc")))
putOnSingleProc = atoi(value);
else if ((value = utl::getValue(child,"reuse_interpolation")))
reuseInterp = atoi(value);
else if ((value = utl::getValue(child,"reusable")))
reusable = atoi(value);
else if ((value = utl::getValue(child,"keep_agg_info")))
keepAggInfo = atoi(value);
else if ((value = utl::getValue(child,"threshold")))
threshold = atof(value);
else if ((value = utl::getValue(child,"damping_factor")))
dampingFactor = atof(value);
else if ((value = utl::getValue(child,"repartition_ratio")))
repartitionRatio = atof(value);
else if ((value = utl::getValue(child,"aux")))
aux = atoi(value);
else if ((value = utl::getValue(child,"aux_threshold")))
auxThreshold = atof(value);
}
void LinSolParams::BlockParams::HypreSettings::read(const TiXmlElement* elem)
{
const char* value;
const TiXmlElement* child = elem->FirstChildElement();
for (; child; child = child->NextSiblingElement())
if ((value = utl::getValue(child,"type")))
type = value;
else if ((value = utl::getValue(child,"no_agg_coarse")))
noAggCoarse = atoi(value);
else if ((value = utl::getValue(child,"no_path_agg_coarse")))
noPathAggCoarse = atof(value);
else if ((value = utl::getValue(child,"truncation")))
truncation = atof(value);
else if ((value = utl::getValue(child,"threshold")))
threshold = atof(value);
else if ((value = utl::getValue(child,"coarsen_scheme")))
coarsenScheme = value;
}
bool LinSolParams::read (const TiXmlElement* elem)
{
#ifdef HAS_PETSC
const char* value = 0;
if ((value = utl::getValue(child,"type")))
method = value;
else if ((value = utl::getValue(child,"noStepBeforeReset")))
nResetSolver = atoi(value);
else if ((value = utl::getValue(child,"pc")))
prec = value;
else if ((value = utl::getValue(child,"schurpc"))) {
if (!strncasecmp(value,"SIMPLE",6))
schurPrec = SIMPLE;
else if (!strncasecmp(value,"MSIMPLER",8))
schurPrec = MSIMPLER;
else if (!strncasecmp(value,"PCD",3))
schurPrec = PCD;
}
else if ((value = utl::getValue(child,"subpc"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
subprec.assign(begin,end);
}
else if ((value = utl::getValue(child,"hypretype"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
hypretype.assign(begin, end);
}
else if ((value = utl::getValue(child,"package"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
package.assign(begin, end);
}
else if ((value = utl::getValue(child,"levels"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
levels.assign(begin, end);
}
else if ((value = utl::getValue(child,"mglevels"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
mglevels.assign(begin, end);
}
else if ((value = utl::getValue(child,"noPreSmooth"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
noPreSmooth.assign(begin, end);
}
else if ((value = utl::getValue(child,"noPostSmooth"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
noPostSmooth.assign(begin, end);
}
else if ((value = utl::getValue(child,"noFineSmooth"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
noFineSmooth.assign(begin, end);
}
else if ((value = utl::getValue(child,"presmoother"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
presmoother.assign(begin, end);
}
else if ((value = utl::getValue(child,"postsmoother"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
postsmoother.assign(begin, end);
}
else if ((value = utl::getValue(child,"finesmoother"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
finesmoother.assign(begin, end);
}
else if ((value = utl::getValue(child,"mgksp"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
mgKsp.assign(begin, end);
}
else if (!strcasecmp(child->Value(),"dirsmoother")) {
size_t block, order;
std::string type;
if (!utl::getAttribute(child,"block",block))
return false;
if (!utl::getAttribute(child,"type",type))
return false;
if (!utl::getAttribute(child,"order",order))
return false;
if (block < 1)
return false;
if (dirsmoother.empty() || (dirsmoother.size() < block))
dirsmoother.resize(std::max(subprec.size(),block));
dirsmoother[block-1].push_back(type);
if (dirOrder.empty() || (dirOrder.size() < block))
dirOrder.resize(std::max(subprec.size(),block));
dirOrder[block-1].push_back(order);
}
else if ((value = utl::getValue(child,"overlap"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
overlap.assign(begin, end);
}
else if ((value = utl::getValue(child,"nx"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
nx.assign(begin, end);
}
else if ((value = utl::getValue(child,"ny"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
ny.assign(begin, end);
}
else if ((value = utl::getValue(child,"nz"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
nz.assign(begin, end);
}
else if ((value = utl::getValue(child,"atol")))
atol = atof(value);
else if ((value = utl::getValue(child,"rtol")))
rtol = atof(value);
else if ((value = utl::getValue(child,"dtol")))
dtol = atof(value);
else if ((value = utl::getValue(child,"maxits")))
maxIts = atoi(value);
else if ((value = utl::getValue(child,"maxCoarseSize"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
maxCoarseSize.assign(begin, end);
}
else if ((value = utl::getValue(child,"GAMGtype"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
GAMGtype.assign(begin, end);
}
else if ((value = utl::getValue(child,"GAMGrepartition"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
GAMGrepartition.assign(begin, end);
}
else if ((value = utl::getValue(child,"GAMGuseAggGasm"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
GAMGuseAggGasm.assign(begin, end);
}
else if ((value = utl::getValue(child,"GAMGprocEqLimit"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
GAMGprocEqLimit.assign(begin, end);
}
else if ((value = utl::getValue(child,"GAMGreuseInterpolation"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
GAMGreuseInterp.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLCoarsePackage"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
MLCoarsePackage.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLCoarseSolver"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
MLCoarseSolver.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLCoarsenScheme"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
MLCoarsenScheme.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLSymmetrize"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLSymmetrize.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLRepartition"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLRepartition.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLBlockScaling"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLBlockScaling.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLPutOnSingleProc"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLPutOnSingleProc.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLReuseInterpolation"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLReuseInterp.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLReuseable"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLReusable.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLKeepAggInfo"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLKeepAggInfo.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLThreshold"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscReal> begin(this_line), end;
MLThreshold.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLDampingFactor"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscReal> begin(this_line), end;
MLDampingFactor.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLRepartitionRatio"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscReal> begin(this_line), end;
MLRepartitionRatio.assign(begin, end);
}
else if ((value = utl::getValue(child,"MLAux"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
MLAux.assign(begin, end);
}
else if ((value = utl::getValue(child,"HypreNoAggCoarse"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
HypreNoAggCoarse.assign(begin, end);
}
else if ((value = utl::getValue(child,"HypreNoPathAggCoarse"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscInt> begin(this_line), end;
HypreNoPathAggCoarse.assign(begin, end);
}
else if ((value = utl::getValue(child,"HypreTruncation"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscReal> begin(this_line), end;
HypreTruncation.assign(begin, end);
}
else if ((value = utl::getValue(child,"HypreThreshold"))) {
std::istringstream this_line(value);
std::istream_iterator<PetscReal> begin(this_line), end;
HypreThreshold.assign(begin, end);
}
else if ((value = utl::getValue(child,"HypreCoarsenScheme"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
HypreCoarsenScheme.assign(begin, end);
}
else if ((value = utl::getValue(child,"ncomponents"))) {
std::istringstream this_line(value);
std::istream_iterator<int> begin(this_line), end;
ncomps.assign(begin, end);
nblock = ncomps.size();
}
else if ((value = utl::getValue(child,"nullspace"))) {
std::istringstream this_line(value);
std::istream_iterator<std::string> begin(this_line), end;
std::istream_iterator<std::string> it;
nullspc.clear();
for (it = begin;it != end;it++) {
if (!strcasecmp(it->c_str(),"constant"))
nullspc.push_back(CONSTANT);
else if (!strcasecmp(value,"rigid_body"))
nullspc.push_back(RIGID_BODY);
else
nullspc.push_back(NONE);
const TiXmlElement* child = elem->FirstChildElement();
for (; child; child = child->NextSiblingElement())
if ((value = utl::getValue(child,"type")))
method = value;
else if ((value = utl::getValue(child,"noStepBeforeReset")))
nResetSolver = atoi(value);
else if ((value = utl::getValue(child,"pc")))
prec = value;
#ifdef HAS_PETSC
else if ((value = utl::getValue(child,"schurpc"))) {
if (!strncasecmp(value,"SIMPLE",6))
schurPrec = SIMPLE;
else if (!strncasecmp(value,"MSIMPLER",8))
schurPrec = MSIMPLER;
else if (!strncasecmp(value,"PCD",3))
schurPrec = PCD;
}
}
else
{
std::cerr <<" *** LinSolParams::read: Unknown keyword: "
<< child->Value() << std::endl;
return false;
}
#endif
else if (!strcasecmp(child->Value(),"block")) {
blocks.resize(blocks.size()+1);
blocks.back().read(child);
}
else if ((value = utl::getValue(child,"atol")))
atol = atof(value);
else if ((value = utl::getValue(child,"rtol")))
rtol = atof(value);
else if ((value = utl::getValue(child,"dtol")))
dtol = atof(value);
else if ((value = utl::getValue(child,"maxits")))
maxIts = atoi(value);
if (blocks.size() == 0) {
blocks.resize(1);
blocks.back().read(elem);
}
return true;
}
@@ -467,25 +231,6 @@ bool LinSolParams::read (const char* filename)
#ifdef HAS_PETSC
const char* LinSolParams::getPreconditioner (int i) const
{
const char* precon = i < 0 ? prec.c_str() : subprec[i].c_str();
return strncasecmp(precon,"asmlu",5) == 0 ? "asm" : precon;
}
int LinSolParams::getLocalPartitioning(size_t dir, size_t i) const
{
switch (dir) {
case 0: return nx[i];
case 1: return ny[i];
case 2: return nz[i];
default: return 0;
}
}
void LinSolParams::setParams (KSP& ksp, PetscIntMat& locSubdDofs,
PetscIntMat& subdDofs, PetscRealVec& coords,
ISMat& dirIndexSet) const
@@ -525,7 +270,7 @@ void LinSolParams::setParams (KSP& ksp, PetscIntMat& locSubdDofs,
}
if (!strncasecmp(prec.c_str(),"asm",3) || !strncasecmp(prec.c_str(),"gasm",4))
setupAdditiveSchwarz(pc, overlap[0], !strncasecmp(prec.c_str(),"asmlu",5),
setupAdditiveSchwarz(pc, blocks[0].overlap, !strncasecmp(prec.c_str(),"asmlu",5),
locSubdDofs, subdDofs, false);
else if (!strncasecmp(prec.c_str(),"ml",2) || !strncasecmp(prec.c_str(),"gamg",4)) {
if (!strncasecmp(prec.c_str(),"ml",2))
@@ -537,7 +282,7 @@ void LinSolParams::setParams (KSP& ksp, PetscIntMat& locSubdDofs,
PCSetUp(pc);
// Settings for coarse solver
if (!MLCoarseSolver.empty())
if (!blocks[0].ml.coarseSolver.empty())
setupCoarseSolver(pc, "", 0);
setupSmoothers(pc, 0, dirIndexSet, locSubdDofs, subdDofs);
@@ -561,9 +306,9 @@ bool LinSolParams::addDirSmoother(PC pc, Mat P, int block,
{
PCSetType(pc,"composite");
PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE);
for (size_t k = 0;k < dirsmoother[block].size();k++)
for (size_t k = 0;k < blocks[block].dirSmoother.size();k++)
PCCompositeAddPC(pc,"shell");
for (size_t k = 0;k < dirsmoother[block].size();k++) {
for (size_t k = 0;k < blocks[block].dirSmoother.size();k++) {
PC dirpc;
PCCompositeGetPC(pc,k,&dirpc);
PCPerm *pcperm;
@@ -572,7 +317,7 @@ bool LinSolParams::addDirSmoother(PC pc, Mat P, int block,
PCShellSetContext(dirpc,pcperm);
PCShellSetDestroy(dirpc,PCPermDestroy);
PCShellSetName(dirpc,"dir");
PCPermSetUp(dirpc,&dirIndexSet[block][k],P,dirsmoother[block][k].c_str());
PCPermSetUp(dirpc,&dirIndexSet[block][k],P,blocks[block].dirSmoother[k].type.c_str());
}
return true;
@@ -602,111 +347,112 @@ static std::string ToString(T data)
void LinSolParams::setMLOptions(const std::string& prefix, int block) const
{
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_maxNLevels").c_str(),
ToString(mglevels[block]).c_str());
ToString(blocks[block].mglevel).c_str());
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_maxCoarseSize").c_str(),
ToString(maxCoarseSize[block]).c_str());
if (!MLCoarsenScheme.empty())
ToString(blocks[block].maxCoarseSize).c_str());
if (!blocks[block].ml.coarsenScheme.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_CoarsenScheme").c_str(),
MLCoarsenScheme[block].c_str());
if (!MLThreshold.empty())
blocks[block].ml.coarsenScheme.c_str());
if (blocks[block].ml.threshold > -1.0)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_Threshold").c_str(),
ToString(MLThreshold[block]).c_str());
if (!MLDampingFactor.empty())
ToString(blocks[block].ml.threshold).c_str());
if (blocks[block].ml.dampingFactor > -1.0)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_DampingFactor").c_str(),
ToString(MLDampingFactor[block]).c_str());
if (!MLRepartitionRatio.empty())
ToString(blocks[block].ml.dampingFactor).c_str());
if (blocks[block].ml.repartitionRatio > -1.0)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_repartitionMaxMinRatio").c_str(),
ToString(MLRepartitionRatio[block]).c_str());
if (!MLSymmetrize.empty())
ToString(blocks[block].ml.repartitionRatio).c_str());
if (blocks[block].ml.symmetrize > -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_Symmetrize").c_str(),
ToString(MLSymmetrize[block]).c_str());
if (!MLRepartition.empty())
ToString(blocks[block].ml.symmetrize).c_str());
if (blocks[block].ml.repartition > -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_repartition").c_str(),
ToString(MLRepartition[block]).c_str());
if (!MLBlockScaling.empty())
ToString(blocks[block].ml.repartition).c_str());
if (blocks[block].ml.blockScaling > -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_BlockScaling").c_str(),
ToString(MLBlockScaling[block]).c_str());
if (!MLPutOnSingleProc.empty())
ToString(blocks[block].ml.blockScaling).c_str());
if (blocks[block].ml.putOnSingleProc > -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_repartitionPutOnSingleProc").c_str(),
ToString(MLPutOnSingleProc[block]).c_str());
if (!MLReuseInterp.empty())
ToString(blocks[block].ml.putOnSingleProc).c_str());
if (blocks[block].ml.reuseInterp > -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_reuse_interpolation").c_str(),
ToString(MLReuseInterp[block]).c_str());
if (!MLKeepAggInfo.empty())
ToString(blocks[block].ml.reuseInterp).c_str());
if (blocks[block].ml.keepAggInfo> -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_KeepAggInfo").c_str(),
ToString(MLKeepAggInfo[block]).c_str());
if (!MLReusable.empty())
ToString(blocks[block].ml.keepAggInfo).c_str());
if (blocks[block].ml.reusable > -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_Reusable").c_str(),
ToString(MLReusable[block]).c_str());
if (!MLAux.empty()) {
ToString(blocks[block].ml.reusable).c_str());
if (blocks[block].ml.aux > -1) {
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_Aux").c_str(),
ToString(MLAux[block]).c_str());
if (!MLThreshold.empty())
ToString(blocks[block].ml.aux).c_str());
if (blocks[block].ml.auxThreshold > -1)
PetscOptionsSetValue(AddPrefix(prefix,"pc_ml_AuxThreshold").c_str(),
ToString(MLThreshold[block]).c_str());
ToString(blocks[block].ml.auxThreshold).c_str());
}
}
void LinSolParams::setGAMGOptions(const std::string& prefix, int block) const
template<class T>
static void condSetup(const std::string& prefix, const std::string& option, const T& var)
{
if (!maxCoarseSize.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_gamg_coarse_eq_limit").c_str(),
ToString(maxCoarseSize[block]).c_str());
if (!GAMGprocEqLimit.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_gamg_process_eq_limit").c_str(),
ToString(GAMGprocEqLimit[block]).c_str());
if (!GAMGtype.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_gamg_type").c_str(),
GAMGtype[block].c_str());
if (!GAMGrepartition.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_gamg_repartition").c_str(),
ToString(GAMGrepartition[block]).c_str());
if (!GAMGuseAggGasm.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_gamg_use_agg_gasm").c_str(),
ToString(GAMGuseAggGasm[block]).c_str());
if (!GAMGreuseInterp.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_gamg_reuse_interpolation").c_str(),
ToString(GAMGreuseInterp[block]).c_str());
if (!MLThreshold.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_gamg_threshold").c_str(),
ToString(MLThreshold[block]).c_str());
if (!mglevels.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_mg_levels").c_str(),
ToString(mglevels[block]).c_str());
if (var > -1)
PetscOptionsSetValue(AddPrefix(prefix,option).c_str(), ToString(var).c_str());
};
template<>
void condSetup(const std::string& prefix, const std::string& option, const std::string& var)
{
if (!var.empty())
PetscOptionsSetValue(AddPrefix(prefix,option).c_str(), var.c_str());
};
void LinSolParams::setGAMGOptions(const std::string& prefix, size_t iblock) const
{
if (iblock >= blocks.size())
return;
const BlockParams& block = blocks[iblock];
condSetup(prefix, "pc_gamg_type", block.gamg.type);
condSetup(prefix, "pc_gamg_coarse_eq_limit", block.maxCoarseSize);
condSetup(prefix, "pc_gamg_process_eq_limit", block.gamg.procEqLimit);
condSetup(prefix, "pc_gamg_repartition", block.gamg.repartition);
condSetup(prefix, "pc_gamg_use_agg_gasm", block.gamg.useAggGasm);
condSetup(prefix, "pc_gamg_reuse_interpolation", block.gamg.reuseInterp);
condSetup(prefix, "pc_gamg_threshold", block.gamg.threshold);
condSetup(prefix, "pc_mg_levels", block.mglevel);
}
void LinSolParams::setHypreOptions(const std::string& prefix, int block) const
void LinSolParams::setHypreOptions(const std::string& prefix, size_t iblock) const
{
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_type").c_str(),
hypretype[block].c_str());
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_max_levels").c_str(),
ToString(mglevels[block]).c_str());
if (iblock >= blocks.size())
return;
const BlockParams& block = blocks[iblock];
condSetup(prefix,"pc_hypre_type", block.hypre.type);
condSetup(prefix, "pc_hypre_boomeramg_max_levels", block.mglevel);
// TODO: Investigate why these are hardcoded.
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_max_iter").c_str(), "1");
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_tol").c_str(), "0.0");;
if (!HypreThreshold.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_strong_threshold").c_str(),
ToString(HypreThreshold[block]).c_str());
if (!HypreCoarsenScheme.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_coarsen_type").c_str(),
ToString(HypreCoarsenScheme[block]).c_str());
if (!HypreNoAggCoarse.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_agg_nl").c_str(),
ToString(HypreNoAggCoarse[block]).c_str());
if (!HypreNoPathAggCoarse.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_agg_num_paths").c_str(),
ToString(HypreNoPathAggCoarse[block]).c_str());
if (!HypreTruncation.empty())
PetscOptionsSetValue(AddPrefix(prefix,"pc_hypre_boomeramg_truncfactor").c_str(),
ToString(HypreTruncation[block]).c_str());
condSetup(prefix, "pc_hypre_boomeramg_strong_threshold", block.hypre.threshold);
condSetup(prefix, "pc_hypre_boomeramg_coarsen_type", block.hypre.coarsenScheme);
condSetup(prefix, "pc_hypre_boomeramg_agg_nl", block.hypre.noAggCoarse);
condSetup(prefix, "pc_hypre_boomeramg_agg_num_paths", block.hypre.noPathAggCoarse);
condSetup(prefix, "pc_hypre_boomeramg_truncfactor", block.hypre.truncation);
}
void LinSolParams::setupCoarseSolver(PC& pc, const std::string& prefix, int block) const
void LinSolParams::setupCoarseSolver(PC& pc, const std::string& prefix, size_t iblock) const
{
if (MLCoarseSolver[block] == "OneLevelSchwarz" ||
MLCoarseSolver[block] == "TwoLevelSchwarz") {
if (iblock >= blocks.size())
return;
const BlockParams& block = blocks[iblock];
if (block.ml.coarseSolver == "OneLevelSchwarz" ||
block.ml.coarseSolver == "TwoLevelSchwarz") {
KSP cksp;
PC cpc;
PCMGGetCoarseSolve(pc,&cksp);
@@ -720,7 +466,7 @@ void LinSolParams::setupCoarseSolver(PC& pc, const std::string& prefix, int bloc
PCRedistributeGetKSP(cpc,&sksp);
KSPSetTolerances(sksp,1.0e-2,PETSC_DEFAULT,PETSC_DEFAULT,10);
KSPGetPC(sksp,&spc);
if (MLCoarseSolver[block] == "OneLevelSchwarz") {
if (block.ml.coarseSolver == "OneLevelSchwarz") {
PCSetType(spc,PCGASM);
PCSetUp(spc);
} else {
@@ -749,22 +495,26 @@ void LinSolParams::setupCoarseSolver(PC& pc, const std::string& prefix, int bloc
}
}
if (!MLCoarsePackage.empty()) {
if (!block.ml.coarsePackage.empty()) {
KSP cksp;
PC cpc;
PCMGGetCoarseSolve(pc,&cksp);
KSPGetPC(cksp,&cpc);
PCSetType(cpc,PCLU);
PCFactorSetMatSolverPackage(cpc,MLCoarsePackage[block].c_str());
PCFactorSetMatSolverPackage(cpc,block.ml.coarsePackage.c_str());
PCSetUp(cpc);
}
}
void LinSolParams::setupSmoothers(PC& pc, int block, ISMat& dirIndexSet,
void LinSolParams::setupSmoothers(PC& pc, size_t iblock, ISMat& dirIndexSet,
const PetscIntMat& locSubdDofs,
const PetscIntMat& subdDofs) const
{
if (iblock >= blocks.size())
return;
const BlockParams& block = blocks[iblock];
PetscInt n;
PCMGGetLevels(pc,&n);
@@ -779,34 +529,32 @@ void LinSolParams::setupSmoothers(PC& pc, int block, ISMat& dirIndexSet,
PCMGGetSmoother(pc,i,&preksp);
std::string compare = mgKsp[(size_t)block<mgKsp.size()-1?block:0];
// warn that richardson might break symmetry if the KSP is CG
if (compare == "defrichardson" && method == KSPCG)
if (block.mgKSP == "defrichardson" && method == KSPCG)
std::cerr << "WARNING: Using multigrid with Richardson on sublevels.\n"
<< "If you get divergence with KSP_DIVERGED_INDEFINITE_PC, try\n"
<< "adding <mgksp>chebyshev</mgksp. Add <mgksp>richardson</mgksp>\n"
<< "to quell this warning." << std::endl;
if (compare == "richardson" || compare == "defrichardson")
if (block.mgKSP == "richardson" || block.mgKSP == "defrichardson")
KSPSetType(preksp,KSPRICHARDSON);
else if (compare == "chebyshev")
else if (block.mgKSP == "chebyshev")
KSPSetType(preksp,KSPCHEBYSHEV);
if ((i == n-1) && (!finesmoother.empty())) {
smoother = finesmoother[block];
noSmooth = noFineSmooth[block];
if ((i == n-1) && (!block.finesmoother.empty())) {
smoother = block.finesmoother;
noSmooth = block.noFineSmooth;
}
else {
smoother = presmoother[block];
noSmooth = noPreSmooth[block];
smoother = block.presmoother;
noSmooth = block.noPreSmooth;
}
KSPSetTolerances(preksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,noSmooth);
KSPGetPC(preksp,&prepc);
if (smoother == "asm" || smoother == "asmlu")
setupAdditiveSchwarz(prepc, overlap[0], smoother == "asmlu",
setupAdditiveSchwarz(prepc, block.overlap, smoother == "asmlu",
locSubdDofs, subdDofs, true);
else if (smoother == "compositedir" && (i==n-1)) {
Mat mat;
@@ -818,12 +566,12 @@ void LinSolParams::setupSmoothers(PC& pc, int block, ISMat& dirIndexSet,
PCGetOperators(prepc,&mat,&Pmat);
#endif
addDirSmoother(prepc,Pmat,block,dirIndexSet);
addDirSmoother(prepc,Pmat,iblock,dirIndexSet);
}
else
PCSetType(prepc,smoother.c_str());
PCFactorSetLevels(prepc,levels[block]);
PCFactorSetLevels(prepc,block.ilu_fill_level);
KSPSetUp(preksp);
}
}

288
src/LinAlg/LinSolParams.h
View File

@@ -16,10 +16,12 @@
#ifndef _LINSOLPARAMS_H
#define _LINSOLPARAMS_H
#include "PETScSupport.h"
#include <array>
#include <iostream>
#include <string>
#include <vector>
#include "PETScSupport.h"
#ifdef HAS_PETSC
@@ -41,6 +43,13 @@ enum SchurPrec { SIMPLE, MSIMPLER, PCD };
class TiXmlElement;
#ifdef HAS_PETSC
#define BLANK_IF_NO_PETSC(a) a
#else
#define BLANK_IF_NO_PETSC(a) ""
#endif
/*!
\brief Class for linear solver parameters.
\details It contains information about solver method, preconditioner
@@ -51,10 +60,18 @@ class LinSolParams
{
public:
//! \brief Default constructor.
LinSolParams(PetscInt n, int m = 1) : nsd(n), msgLev(m) { this->setDefault(); }
//! \brief Set default values.
void setDefault();
LinSolParams(PetscInt n, int m = 1) :
atol(1e-6),
rtol(1e-6),
dtol(1e3),
maxIts(1000),
nResetSolver(0),
#ifdef HAS_PETSC
schurPrec(SIMPLE),
#endif
nsd(n),
msgLev(m)
{}
//! \brief Read linear solver parameters from stream.
bool read(std::istream& is, int nparams = 10);
@@ -65,49 +82,151 @@ public:
//! \brief Read linear solver parameters from file.
bool read(const char* filename);
//! \brief Linear solver settings for a block of the linear system
struct BlockParams {
BlockParams() :
package(BLANK_IF_NO_PETSC(MATSOLVERPETSC)),
basis(1),
comps(0),
overlap(1),
ilu_fill_level(0),
mglevel(-1),
noPreSmooth(-1),
noFineSmooth(-1),
noPostSmooth(-1),
presmoother(BLANK_IF_NO_PETSC(PCILU)),
postsmoother(BLANK_IF_NO_PETSC(PCILU)),
mgKSP("defrichardson"),
maxCoarseSize(-1)
#ifdef HAS_PETSC
//! \brief Set linear solver method
void setMethod(KSPType type) { method.assign(type); }
, nullspace(NONE)
#endif
{}
//! \brief Set preconditioner
void setPreconditioner(PCType type) { prec.assign(type); }
//! \brief Settings for a directional smoother
struct DirSmoother {
int order; //!< Ordering of DOFs
std::string type; //!< Directional smoother types
//! \brief Set preconditioner for sub-block
void setSubPreconditioner(PCType type, size_t i = 0) { subprec[i].assign(type); }
DirSmoother(int o, const std::string& t) : order(o), type(t) {}
};
//! \brief Set linear solver package
void setPackage(MatSolverPackage stype, size_t i = 0) { package[i].assign(stype); }
//! \brief Structure holding settings for a GAMG preconditioner
struct GAMGSettings {
//! \brief Default constructor
GAMGSettings() :
procEqLimit(-1),
repartition(-1),
useAggGasm(-1),
reuseInterp(-1),
threshold(-1)
{}
//! \brief Set absolute convergence tolerance
void setAbsTolerance(Real eps) { atol = eps; }
std::string type; //!< Type for GAMG multigrid method ("geo" or "agg")
PetscInt procEqLimit; //!< Limit number of equations on each process on coarse grid
PetscInt repartition; //!< Repartition coarse grid for GAMG preconditioner
PetscInt useAggGasm; //!< Use aggregation aggregates for GASM smoother
PetscInt reuseInterp; //!< Reuse interpolation
PetscReal threshold; //!< Smoother drop toleranse
//! \brief Set relative convergence tolerance
void setRelTolerance(Real eps) { rtol = eps; }
//! \brief Read linear solver parameters from XML document.
void read(const TiXmlElement* elem);
};
//! \brief Set divergence tolerance
void setDivTolerance(Real eps) { dtol = eps; }
//! \brief Structure holding settings for a ML preconditioner
struct MLSettings {
//! \brief Default constructor
MLSettings() :
symmetrize(-1),
reuseInterp(-1),
blockScaling(-1),
putOnSingleProc(-1),
threshold(-1),
dampingFactor(-1),
repartitionRatio(-1),
reusable(-1),
repartition(-1),
keepAggInfo(-1),
aux(-1),
auxThreshold(-1)
{}
//! \brief Set maximum number of iterations
void setMaxIterations(int its) { maxIts = its; }
std::string coarsePackage; //!< Coarse matrix solver package
std::string coarseSolver; //!< DD type coarse solver for ML
std::string coarsenScheme; //!< Coarsening scheme for ML
PetscInt symmetrize; //!< Symmetrize aggregation
PetscInt reuseInterp; //!< Reuse interpolation operators between solves
PetscInt blockScaling; //!< Scale all dofs at each node together
PetscInt putOnSingleProc; //!< If below threshold, assign to a single processor
PetscReal threshold; //!< Smoother drop toleranse for ML
PetscReal dampingFactor; //!< Damping factor
PetscInt repartitionRatio; //!< Max-min ratio for repartitioning
PetscInt reusable; //!< Store intermediate data structures such that the MG hierarchy is reusable
PetscInt repartition; //!< Repartitioning
PetscInt keepAggInfo; //!< Allow ML preconditioner to be reused
PetscInt aux; //!< Use auxillary coordinate based laplacian for aggregation
PetscReal auxThreshold; //!< Tolerance for auxillary laplacian aggregation
//! \brief Set number of overlap
void setOverlap(int olap, size_t i = 0) { overlap[i] = olap; }
//! \brief Read linear solver parameters from XML document.
void read(const TiXmlElement* elem);
};
//! \brief Set number of local subdomains for each patch
void setLocalPartitioning(size_t NX = 0, size_t NY = 0, size_t NZ = 0, size_t i = 0)
{ nx[i] = NX; ny[i] = NY; nz[i] = NZ; }
//! \brief Structure holding settings for a Hypre preconditioner
struct HypreSettings {
//! \brief Default constructor
HypreSettings() :
type("boomeramg"),
noAggCoarse(-1),
noPathAggCoarse(-1),
truncation(-1),
threshold(-1)
{}
//! \brief Set null-space of matrix
void setNullSpace(NullSpace nspc, size_t i = 0) { nullspc[i] = nspc; }
std::string type; //!< Type of hypre preconditioner
PetscInt noAggCoarse; //!< Number of levels of agressive coarsening
PetscInt noPathAggCoarse; //!< Number of paths for aggressive coarsening
PetscReal truncation; //!< Truncation factor for interpolation
PetscReal threshold; //!< Drop tolerance for Hypre
std::string coarsenScheme; //!< Coarsening scheme for Hypre
//! \brief Read linear solver parameters from XML document.
void read(const TiXmlElement* elem);
};
//! \brief Read settings from XML block
//! \param[in] child XML block
bool read(const TiXmlElement* child);
std::string prec; //!< Preconditioner for block
std::string package; //!< Package providing solver
size_t basis; //!< Basis for block
size_t comps; //!< Components from basis (1, 2, 3, 12, 13, 23, 123, ..., 0 = all)
PetscInt overlap; //!< Overlap in ASM
PetscInt ilu_fill_level; //!< Fill level for ILU
PetscInt mglevel; //!< Levels for multigrid
PetscInt noPreSmooth; //!< Number of presmoothings for AMG
PetscInt noFineSmooth; //!< Number of fine smoothings for AMG
PetscInt noPostSmooth; //!< Number of postsmoothings for AMG
std::string presmoother; //!< Presmoother for AMG
std::string postsmoother; //!< Postsmoother for AMG
std::string finesmoother; //!< Smoother on finest grid
std::string mgKSP; //!< KSP type on MG levels
int maxCoarseSize; //!< Max number of DOFs for coarse AMG system
GAMGSettings gamg; //!< Settings for GAMG preconditioner
MLSettings ml; //!< Settings for ML preconditioner
HypreSettings hypre; //!< Settings for Hypre preconditioner
std::array<size_t,3> subdomains; //!< Number of local subdomains in each parameter direction
std::vector<DirSmoother> dirSmoother; //!< Direction smoothers
#ifdef HAS_PETSC
NullSpace nullspace; //!< Nullspace for matrix
#endif
};
//! \brief Get linear solver method
const char* getMethod() const { return method.c_str(); }
//! \brief Get preconditioner
const char* getPreconditioner(int i = -1) const;
//! \brief Get linear solver package
const char* getPackage(size_t i = 0) const { return package[i].c_str(); }
const std::string& getMethod() const { return method; }
//! \brief Get absolute convergence tolerance
Real getAbsTolerance() const { return atol; }
@@ -118,24 +237,36 @@ public:
//! \brief Get divergence tolerance
Real getDivTolerance() const { return dtol; }
//! \brief Get number of iterations before reset (GMRES)
PetscInt getResetSolver() const { return nResetSolver; }
//! \brief Get maximum number of iterations
int getMaxIterations() const { return maxIts; }
//! \brief Get number of overlaps
int getOverlap(int i = 0) const { return overlap[i]; }
//! \brief Get local partitioning
int getLocalPartitioning(size_t dir = 0, size_t i = 0) const;
//! \brief Number of blocks in matrix system
int getNoBlocks() const { return nblock; }
size_t getNoBlocks() const { return blocks.size(); }
//! \brief Number of components in a matrix block
const IntVec& getComponents() const { return ncomps; }
//! \brief Obtain settings for a given block
const BlockParams& getBlock(size_t i) const { return blocks[i]; }
//! \brief Get number of overlaps
NullSpace getNullSpace(size_t i = 0) const { return i < nullspc.size()?nullspc[i]:NONE; }
//! \brief Get main preconditioner
const char* getPreconditioner() const
{
return (prec == "asmlu" ? "asm" : prec.c_str());
}
#ifdef HAS_PETSC
//! \brief Get schur complement type
SchurPrec getSchurType() const { return schurPrec; }
#endif
//! \brief Get dimensionality of node coordinates
size_t getDimension() const { return nsd; }
//! \brief Get the message level
int getMessageLevel() const { return msgLev; }
#ifdef HAS_PETSC
//! \brief Set linear solver parameters for KSP object
void setParams(KSP& ksp, PetscIntMat& locSubdDofs,
PetscIntMat& subdDofs, PetscRealVec& coords,
@@ -156,18 +287,18 @@ public:
//! \brief Set GAMG options
//! \param[in] prefix The prefix of the block to set parameters for
//! \param[in] block The index of the block to set parameters for
void setGAMGOptions(const std::string& prefix, int block) const;
void setGAMGOptions(const std::string& prefix, size_t block) const;
//! \brief Set Hypre options
//! \param[in] prefix The prefix of the block to set parameters for
//! \param[in] block The index of the block to set parameters for
void setHypreOptions(const std::string& prefix, int block) const;
void setHypreOptions(const std::string& prefix, size_t block) const;
//! \brief Setup the coarse solver in a multigrid
//! \param[in] PC The preconditioner to set coarse solver for
//! \param[in] prefix The prefix of the block to set parameters for
//! \param[in] block The index of the block to set parameters for
void setupCoarseSolver(PC& pc, const std::string& prefix, int block) const;
void setupCoarseSolver(PC& pc, const std::string& prefix, size_t block) const;
//! \brief Setup the smoothers in a multigrid
//! \param[in] PC The preconditioner to set coarse solver for
@@ -175,7 +306,7 @@ public:
//! \param[in] dirIndexSet The index set for direction smoothers
//! \param[in] locSubdDofs Local subdomain DOFs for ASM preconditioners
//! \param[in] subdDofs Subdomain DOFs for ASM preconditioners
void setupSmoothers(PC& pc, int block, ISMat& dirIndexSet,
void setupSmoothers(PC& pc, size_t block, ISMat& dirIndexSet,
const PetscIntMat& locSubdDofs,
const PetscIntMat& subdDofs) const;
@@ -189,69 +320,22 @@ public:
void setupAdditiveSchwarz(PC& pc, int overlap, bool asmlu,
const PetscIntMat& locSubdDofs,
const PetscIntMat& subdDofs, bool smoother) const;
#endif
private:
PetscReal atol; // Absolute tolerance
PetscReal rtol; // Relative tolerance
PetscReal dtol; // Divergence tolerance
PetscInt maxIts; // Maximum number of iterations
PetscInt nResetSolver; // Number of linear solves before reset of KSP/PC
int nblock; // Number of block
bool schur; // Schur complement solver
#ifdef HAS_PETSC
SchurPrec schurPrec; // Preconditioner for Schur system
#endif
std::string method; // Linear solver method
std::string prec; // Preconditioner
StringVec subprec; // Preconditioners for block-system
StringVec hypretype; // Type of hypre preconditioner
StringVec package; // Linear software package (petsc, superlu_dist, ...)
IntVec ncomps; // Components for each fields in block-vector
PetscIntVec overlap; // Number of overlap in ASM
PetscIntVec levels; // Number of levels of fill to use
PetscIntVec mglevels; // Number of levels for MG
PetscIntVec noPreSmooth; // Number of presmoothings for AMG
PetscIntVec noPostSmooth; // Number of postsmoothings for AMG
PetscIntVec noFineSmooth; // Number of fine grid smoothings for AMG
StringVec presmoother; // Presmoother for AMG
StringVec postsmoother; // Postsmoother for AMG
StringVec finesmoother; // Smoother on finest grid
StringVec mgKsp; // KSP type on MG levels ("richardson")
IntVec maxCoarseSize; // Max number of DOFS for coarse AMG system
StringVec MLCoarsePackage; // Coarse matrix solver package
StringVec MLCoarseSolver; // DD type coarse solver for ML
StringVec MLCoarsenScheme; // Coarsening scheme for ML
PetscIntVec MLSymmetrize; // Symmetrize aggregation
PetscIntVec MLReuseInterp; // Reuse interpolation operators between solves
PetscIntVec MLBlockScaling; // Scale all dofs at each node together
PetscIntVec MLPutOnSingleProc; // If below assign to a single processor
PetscRealVec MLThreshold; // Smoother drop toleranse for ML
PetscRealVec MLDampingFactor; // Damping factor
PetscRealVec MLRepartitionRatio;// Max-min ratio for repartitioning
PetscIntVec MLRepartition; // Repartitioning
PetscIntVec MLReusable; // Store intermediate datastructures such that the MG hierarchy is reusable
PetscIntVec MLKeepAggInfo; // Allow ML preconditioner to be reused
PetscIntVec MLAux; // Use auxillary coordinate based laplacian for aggregation
StringVec GAMGtype; // Type for GAMG multigrid method ("geo" or "agg") ("geo")
PetscIntVec GAMGprocEqLimit; // Limit number of equation on each process on coarse grid
PetscIntVec GAMGrepartition; // Repartition coarse grid for GAMG preconditioner (0 = false, 1 = true)
PetscIntVec GAMGuseAggGasm; // Use aggregation aggregates for GASM smoother (0 = false, 1 = true)
PetscIntVec GAMGreuseInterp; // Reuse interpolation
PetscIntVec HypreNoAggCoarse; // Number of levels of aggressive coarsening
PetscIntVec HypreNoPathAggCoarse; // Number of paths for aggressive coarsening
PetscRealVec HypreTruncation; // Truncation factor for interpolation
PetscRealVec HypreThreshold; // Drop tolerance for Hypre
StringVec HypreCoarsenScheme;// Coarsening scheme for Hypre
IntVec nx; // Number of local subdomains in first parameter direction
IntVec ny; // Number of local subdomains in second parameter direction
IntVec nz; // Number of local subdomains in third parameter direction
StringMat dirsmoother; // Directional smoother types
IntMat dirOrder; // Direction smoother orders/numberings
std::vector<NullSpace> nullspc; // Null-space for matrix
#endif // HAS_PETSC
std::vector<BlockParams> blocks; //!< Parameters for each block
PetscInt nsd; //!< Number of space dimensions
int msgLev; //!< Flag for extra output
friend class PETScMatrix;
friend class PETScBlockMatrix;
};
#endif

193
src/LinAlg/PETScBlockMatrix.C
View File

@@ -29,16 +29,15 @@
PETScBlockMatrix::PETScBlockMatrix (const ProcessAdm& padm, const LinSolParams& par) :
PETScMatrix(padm,par,LinAlg::GENERAL_MATRIX)
{
ncomps.resize(par.ncomps.size());
nblocks = ncomps.size();
ncomps.resize(par.getNoBlocks());
nblocks = par.getNoBlocks();
for (size_t i = 0;i < nblocks;i++)
ncomps[i] = par.ncomps[i];
ncomps[i] = par.getBlock(i).comps/10 + 1;
isvec = (IS*) malloc(sizeof(IS)*nblocks);
matvec = (Mat*) malloc(sizeof(Mat)*nblocks*nblocks);
for (size_t m = 0;m < nblocks*nblocks;m++) {
MatCreate(*adm.getCommunicator(),&matvec[m]);
MatSetOption(matvec[m],MAT_STRUCTURALLY_SYMMETRIC,PETSC_TRUE);
MatSetFromOptions(matvec[m]);
}
@@ -50,34 +49,8 @@ PETScBlockMatrix::PETScBlockMatrix (const ProcessAdm& padm, const LinSolParams&
}
PETScBlockMatrix::PETScBlockMatrix (const ProcessAdm& padm, const IntVec& nc, const LinSolParams& par)
: PETScMatrix(padm,par,LinAlg::GENERAL_MATRIX)
{
nblocks = nc.size();
ncomps.resize(nblocks);
for (size_t i = 0;i < nblocks;i++)
ncomps[i] = nc[i];
isvec = (IS*) malloc(sizeof(IS)*nblocks);
matvec = (Mat*) malloc(sizeof(Mat)*nblocks*nblocks);
for (size_t m = 0;m < nblocks*nblocks;m++) {
MatCreate(*adm.getCommunicator(),&matvec[m]);
MatSetFromOptions(matvec[m]);
}
LinAlgInit::increfs();
}
PETScBlockMatrix::~PETScBlockMatrix ()
{
// Deallocation of null space
if (solParams.getNullSpace() == CONSTANT)
MatNullSpaceDestroy(PETSCMANGLE(nsp));
// Deallocation of linear solver object.
KSPDestroy(PETSCMANGLE(ksp));
for (size_t m = 0;m < nblocks*nblocks;m++)
MatDestroy(&matvec[m]);
free(matvec);
@@ -94,8 +67,6 @@ PETScBlockMatrix::~PETScBlockMatrix ()
for (PetscInt i = 0; i < ISsize; i++)
ISDestroy(PETSCMANGLE(elmIS[i]));
delete elmIS;
}
@@ -216,23 +187,31 @@ void PETScBlockMatrix::initAssembly (const SAM& sam, bool)
{
const SAMpatchPara* sampch = dynamic_cast<const SAMpatchPara*>(&sam);
if (!strncasecmp(solParams.prec.c_str(),"gamg",4) || !strncasecmp(solParams.prec.c_str(),"ml",2))
if (!strncasecmp(solParams.getPreconditioner(),"gamg",4) ||
!strncasecmp(solParams.getPreconditioner(),"ml",2))
sampch->getLocalNodeCoordinates(coords);
else for (size_t i = 0; i < solParams.subprec.size(); i++)
if (!strncasecmp(solParams.subprec[i].c_str(),"gamg",4) || !strncasecmp(solParams.subprec[i].c_str(),"ml",2)) {
else for (size_t i = 0; i < solParams.getNoBlocks(); i++)
if (!strncasecmp(solParams.getBlock(i).prec.c_str(),"gamg",4) ||
!strncasecmp(solParams.getBlock(i).prec.c_str(),"ml",2)) {
sampch->getLocalNodeCoordinates(coords);
break;
}
if (!solParams.dirOrder.empty()) {
dirIndexSet.resize(solParams.dirOrder.size());
for (size_t i = 0;i < solParams.dirOrder.size();i++)
for (size_t j = 0;j < solParams.dirOrder[i].size();j++) {
bool dirsmoothers=false;
for (size_t i=0; i < solParams.getNoBlocks();++i)
if (!solParams.getBlock(i).dirSmoother.empty())
dirsmoothers = true;
if (dirsmoothers) {
dirIndexSet.resize(solParams.getNoBlocks());
for (size_t i = 0;i < solParams.getNoBlocks();i++)
for (size_t j = 0;j < solParams.getBlock(i).dirSmoother.size();j++) {
IS permIndex;
if (solParams.dirOrder[i][j] != 123) {
if (solParams.getBlock(i).dirSmoother[j].order != 123) {
PetscIntVec perm;
sampch->getDirOrdering(perm,solParams.dirOrder[i][j],solParams.ncomps[i]);
ISCreateGeneral(*adm.getCommunicator(),perm.size(),&(perm[0]),PETSC_COPY_VALUES,&permIndex);
sampch->getDirOrdering(perm,solParams.getBlock(i).dirSmoother[j].order,
solParams.getBlock(i).comps);
ISCreateGeneral(*adm.getCommunicator(),perm.size(),&(perm[0]),
PETSC_COPY_VALUES,&permIndex);
ISSetPermutation(permIndex);
}
else {
@@ -243,10 +222,10 @@ void PETScBlockMatrix::initAssembly (const SAM& sam, bool)
}
}
int nx = solParams.getLocalPartitioning(0);
int ny = solParams.getLocalPartitioning(1);
int nz = solParams.getLocalPartitioning(2);
int overlap = solParams.getOverlap();
int nx = solParams.getBlock(0).subdomains[0];
int ny = solParams.getBlock(0).subdomains[1];
int nz = solParams.getBlock(0).subdomains[2];
int overlap = solParams.getBlock(1).overlap;
if (nx+ny+nz > 0) {
locSubdDofsBlock.resize(nblocks);
@@ -359,7 +338,7 @@ void PETScBlockMatrix::initAssembly (const SAM& sam, bool)
MatSetOption(matvec[i],MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
#endif
if (solParams.getNullSpace(1) == CONSTANT) {
if (solParams.getBlock(1).nullspace == CONSTANT) {
Vec const_mode;
VecCreate(*adm.getCommunicator(), &const_mode);
VecSetSizes(const_mode, ilast-ifirst, PETSC_DECIDE);
@@ -375,11 +354,12 @@ void PETScBlockMatrix::initAssembly (const SAM& sam, bool)
VecDuplicate(const_mode,&x);
this->renumberRHS(const_mode,x,true);
VecNormalize(x, NULL);
MatNullSpaceCreate(*adm.getCommunicator(),PETSC_FALSE,1,&x,&nsp);
nsp = new MatNullSpace;
MatNullSpaceCreate(*adm.getCommunicator(),PETSC_FALSE,1,&x,nsp);
#if PETSC_VERSION_MAJOR > 6
KSPSetNullSpace(ksp,nsp);
KSPSetNullSpace(ksp,*nsp);
#else
MatSetNullSpace(matvec[1],nsp);
MatSetNullSpace(matvec[1],*nsp);
#endif
VecDestroy(&x);
VecDestroy(&const_mode);
@@ -423,9 +403,9 @@ bool PETScBlockMatrix::solve (SystemVector& B, bool newLHS, Real*)
newLHS = nLinSolves==0; // don't reset preconditioner unless requested
// Reset linear solver
if (nLinSolves && solParams.nResetSolver)
if (nLinSolves%solParams.nResetSolver == 0) {
if (solParams.schur)
if (nLinSolves && solParams.getResetSolver())
if (nLinSolves%solParams.getResetSolver()== 0) {
if (solParams.getNoBlocks() > 1)
MatDestroy(&Sp);
newLHS = true;
}
@@ -455,9 +435,9 @@ bool PETScBlockMatrix::solve (const SystemVector& B, SystemVector& X, bool newLH
newLHS = nLinSolves==0; // don't reset preconditioner unless requested
// Reset linear solver
if (nLinSolves && solParams.nResetSolver)
if (nLinSolves%solParams.nResetSolver == 0) {
if (solParams.schur)
if (nLinSolves && solParams.getResetSolver())
if (nLinSolves%solParams.getResetSolver() == 0) {
if (solParams.getNoBlocks() > 1)
MatDestroy(&Sp);
newLHS = true;
}
@@ -541,19 +521,22 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
return false;
// Set linear solver method
KSPSetType(ksp,solParams.method.c_str());
KSPSetTolerances(ksp,solParams.rtol,solParams.atol,solParams.dtol,solParams.maxIts);
KSPSetType(ksp,solParams.getMethod().c_str());
KSPSetTolerances(ksp,solParams.getRelTolerance(),
solParams.getAbsTolerance(),
solParams.getDivTolerance(),
solParams.getMaxIterations());
// Set preconditioner
PC pc;
KSPGetPC(ksp,&pc);
PCSetType(pc,solParams.getPreconditioner());
if (!strncasecmp(solParams.prec.c_str(),"gamg",4) || !strncasecmp(solParams.prec.c_str(),"ml",2)) {
PetscInt nloc = coords.size()/solParams.nsd;
PCSetCoordinates(pc,solParams.nsd,nloc,&coords[0]);
if (!strncasecmp(solParams.getPreconditioner(),"gamg",4) ||
!strncasecmp(solParams.getPreconditioner(),"ml",2)) {
PetscInt nloc = coords.size()/solParams.getDimension();
PCSetCoordinates(pc,solParams.getDimension(),nloc,&coords[0]);
}
//PCFactorSetMatSolverPackage(pc,solParams.package[0].c_str());
if (!strncasecmp(solParams.prec.c_str(),"fieldsplit",10)) {
if (!strncasecmp(solParams.getPreconditioner(),"fieldsplit",10)) {
PetscInt m1, m2, n1, n2, nr, nc, nsplit;
KSP *subksp;
PC subpc[2];
@@ -565,12 +548,10 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
else
VecCreateSeq(PETSC_COMM_SELF,m1,&diagA00);
if (!P || (solParams.schurPrec == SIMPLE)) {
if (!P || (solParams.getSchurType() == SIMPLE))
MatGetDiagonal(matvec[0],diagA00);
}
else {
else
MatGetDiagonal(P->getMatrixBlock(0,0),diagA00);
}
VecReciprocal(diagA00);
VecScale(diagA00,-1.0);
@@ -597,14 +578,14 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
KSPSetReusePreconditioner(subksp[1], PETSC_TRUE);
#endif
if (P && Pb && (solParams.schurPrec == PCD)) {
if (P && Pb && (solParams.getSchurType() == PCD)) {
KSPSetType(subksp[1],"preonly");
KSPGetPC(subksp[1],&subpc[1]);
PCSetType(subpc[1],PCSHELL);
PCCreate(*adm.getCommunicator(),&S);
PCCreate(*adm.getCommunicator(),&Fp);
const char* prec1 = solParams.getPreconditioner(1);
const char* prec1 = solParams.getBlock(1).prec.c_str();
if (strncasecmp(prec1,"compositedir",12))
PCSetType(S,prec1);
else if (!solParams.addDirSmoother(S,Sp,1,dirIndexSet))
@@ -642,14 +623,14 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
else if (!solParams.addDirSmoother(S,Sp,1,dirIndexSet))
return false;
PCFactorSetLevels(S,solParams.levels[1]);
if (!strncasecmp(solParams.subprec[1].c_str(),"asm",3)) {
solParams.setupAdditiveSchwarz(S, solParams.overlap[1],
!strncasecmp(solParams.subprec[1].c_str(),"asmlu",5),
PCFactorSetLevels(S,solParams.getBlock(1).ilu_fill_level);
if (!strncasecmp(solParams.getBlock(1).prec.c_str(),"asm",3)) {
solParams.setupAdditiveSchwarz(S, solParams.getBlock(1).overlap,
!strncasecmp(solParams.getBlock(1).prec.c_str(),"asmlu",5),
locSubdDofsBlock.empty()?PetscIntMat():locSubdDofsBlock[1],
subdDofsBlock.empty()?PetscIntMat():subdDofsBlock[1],false);
}
else if (!strncasecmp(solParams.subprec[1].c_str(),"ml",2)) {
else if (!strncasecmp(solParams.getBlock(1).prec.c_str(),"ml",2)) {
PetscInt n;
solParams.setMLOptions("fieldsplit_p", 1);
@@ -657,12 +638,12 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
PCSetUp(S);
// Set coarse solver package
if (!solParams.MLCoarsePackage.empty()) {
if (!solParams.getBlock(1).ml.coarsePackage.empty()) {
KSP cksp;
PC cpc;
PCMGGetCoarseSolve(S,&cksp);
KSPGetPC(cksp,&cpc);
PCFactorSetMatSolverPackage(cpc,solParams.MLCoarsePackage[1].c_str());
PCFactorSetMatSolverPackage(cpc,solParams.getBlock(1).ml.coarsePackage.c_str());
}
PCMGGetLevels(S,&n);
@@ -678,20 +659,20 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
PCMGGetSmoother(S,i,&preksp);
KSPSetType(preksp,"richardson");
if ((i == n-1) && (!solParams.finesmoother.empty())) {
smoother = solParams.finesmoother[1];
noSmooth = solParams.noFineSmooth[1];
if ((i == n-1) && (!solParams.getBlock(1).finesmoother.empty())) {
smoother = solParams.getBlock(1).finesmoother;
noSmooth = solParams.getBlock(1).noFineSmooth;
}
else {
smoother = solParams.presmoother[1];
noSmooth = solParams.noPreSmooth[1];
smoother = solParams.getBlock(1).presmoother;
noSmooth = solParams.getBlock(1).noPreSmooth;
}
KSPSetTolerances(preksp,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,noSmooth);
KSPGetPC(preksp,&prepc);
if (smoother == "asm" || smoother == "asmlu" ) {
solParams.setupAdditiveSchwarz(prepc, solParams.overlap[1],
solParams.setupAdditiveSchwarz(prepc, solParams.getBlock(1).overlap,
smoother == "asmlu",
locSubdDofs, subdDofs, true);
}
@@ -702,13 +683,14 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
else
PCSetType(prepc,smoother.c_str());
PCFactorSetLevels(prepc,solParams.levels[1]);
PCFactorSetLevels(prepc,solParams.getBlock(1).ilu_fill_level);
KSPSetUp(preksp);
}
}
else if (!strncasecmp(solParams.subprec[1].c_str(),"gamg",4) || !strncasecmp(solParams.subprec[1].c_str(),"ml",2) ) {
PetscInt nloc = coords.size()/solParams.nsd;
PCSetCoordinates(pc,solParams.nsd,nloc,&coords[0]);
else if (!strncasecmp(solParams.getBlock(1).prec.c_str(),"gamg",4) ||
!strncasecmp(solParams.getBlock(1).prec.c_str(),"ml",2) ) {
PetscInt nloc = coords.size()/solParams.getDimension();
PCSetCoordinates(pc,solParams.getDimension(),nloc,&coords[0]);
}
PCProdSetUp(subpc[1],&QpL,&Fp,&S);
@@ -716,16 +698,20 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
}
// Preconditioner for blocks
for (PetscInt m = 0; m < nsplit; m++) {
char pchar='1';
for (PetscInt m = 0; m < nsplit; m++, pchar++) {
std::string prefix;
if (m == 0)
prefix = "fieldsplit_u";
else
prefix = "fieldsplit_p";
if (nsplit == 2) {
if (m == 0)
prefix = "fieldsplit_u";
else
prefix = "fieldsplit_p";
} else
prefix = std::string("fieldsplit_b")+pchar;
KSPSetType(subksp[m],"preonly");
KSPGetPC(subksp[m],&subpc[m]);
if (!strncasecmp(solParams.subprec[m].c_str(),"compositedir",12)) {
if (!strncasecmp(solParams.getBlock(m).prec.c_str(),"compositedir",12)) {
Mat mat;
Mat Pmat;
#if PETSC_VERSION_MINOR < 5
@@ -739,35 +725,36 @@ bool PETScBlockMatrix::setParameters(PETScMatrix *P2, PETScVector *Pb)
return false;
}
else
PCSetType(subpc[m],solParams.getPreconditioner(m));
PCSetType(subpc[m],solParams.getBlock(m).prec.c_str());
if (!strncasecmp(solParams.subprec[m].c_str(),"gamg",4) || !strncasecmp(solParams.subprec[m].c_str(),"ml",2) ) {
PetscInt nloc = coords.size()/solParams.nsd;
PCSetCoordinates(subpc[m],solParams.nsd,nloc,&coords[0]);
if (!strncasecmp(solParams.getBlock(m).prec.c_str(),"gamg",4) ||
!strncasecmp(solParams.getBlock(m).prec.c_str(),"ml",2) ) {
PetscInt nloc = coords.size()/solParams.getDimension();
PCSetCoordinates(subpc[m],solParams.getDimension(),nloc,&coords[0]);
PCGAMGSetType(subpc[m], PCGAMGAGG); // TODO?
}
PCFactorSetLevels(subpc[m],solParams.levels[m]);
if (!strncasecmp(solParams.subprec[m].c_str(),"asm",3)) {
solParams.setupAdditiveSchwarz(subpc[m], solParams.overlap[m],
!strncasecmp(solParams.subprec[m].c_str(),"asmlu",5),
PCFactorSetLevels(subpc[m],solParams.getBlock(m).ilu_fill_level);
if (!strncasecmp(solParams.getBlock(m).prec.c_str(),"asm",3)) {
solParams.setupAdditiveSchwarz(subpc[m], solParams.getBlock(m).overlap,
!strncasecmp(solParams.getBlock(m).prec.c_str(),"asmlu",5),
locSubdDofsBlock.empty()?PetscIntMat():locSubdDofsBlock[m],
subdDofsBlock.empty()?PetscIntMat():subdDofsBlock[m],false);
} else if (!strncasecmp(solParams.subprec[m].c_str(),"ml",2)) {
} else if (!strncasecmp(solParams.getBlock(m).prec.c_str(),"ml",2)) {
solParams.setMLOptions(prefix, m);
PCSetFromOptions(subpc[m]);
PCSetUp(subpc[m]);
// Settings for coarse solver
if (!solParams.MLCoarseSolver.empty())
if (!solParams.getBlock(m).ml.coarseSolver.empty())
solParams.setupCoarseSolver(subpc[m], prefix, m);
solParams.setupSmoothers(subpc[m], m, dirIndexSet,
locSubdDofsBlock.empty()?PetscIntMat():locSubdDofsBlock[m],
subdDofsBlock.empty()?PetscIntMat():subdDofsBlock[m]);
}
else if (!strncasecmp(solParams.subprec[m].c_str(),"hypre",5))
else if (!strncasecmp(solParams.getBlock(m).prec.c_str(),"hypre",5))
solParams.setHypreOptions(prefix,m);
}
}

4
src/LinAlg/PETScBlockMatrix.h
View File

@@ -36,10 +36,6 @@ public:
#ifdef HAS_PETSC
//! \brief Constructor.
PETScBlockMatrix(const ProcessAdm& padm, const LinSolParams& spar);
//! \brief Constructor defining the blocks.
//! \param[in] ncomp Number of components in each block
//! \param[in] spar Linear solver parameters
PETScBlockMatrix(const ProcessAdm& padm, const IntVec& ncomp, const LinSolParams& spar);
//! \brief Destructor
virtual ~PETScBlockMatrix();
#else

62
src/LinAlg/PETScMatrix.C
View File

@@ -161,7 +161,7 @@ Real PETScVector::Linfnorm() const
PETScMatrix::PETScMatrix (const ProcessAdm& padm, const LinSolParams& spar,
LinAlg::LinearSystemType ltype) :
adm(padm), solParams(spar), linsysType(ltype)
nsp(nullptr), adm(padm), solParams(spar), linsysType(ltype)
{
// Create matrix object, by default the matrix type is AIJ
MatCreate(*adm.getCommunicator(),&A);
@@ -170,12 +170,13 @@ PETScMatrix::PETScMatrix (const ProcessAdm& padm, const LinSolParams& spar,
KSPCreate(*adm.getCommunicator(),&ksp);
// Create null space if any
if (solParams.getNullSpace() == CONSTANT) {
MatNullSpaceCreate(*adm.getCommunicator(),PETSC_TRUE,0,0,&nsp);
if (solParams.getBlock(0).nullspace == CONSTANT) {
nsp = new MatNullSpace;
MatNullSpaceCreate(*adm.getCommunicator(),PETSC_TRUE,0,0,nsp);
#if PETSC_VERSION_MINOR < 6
KSPSetNullSpace(ksp,nsp);
KSPSetNullSpace(ksp,*nsp);
#else
MatSetNullSpace(A,nsp);
MatSetNullSpace(A,*nsp);
#endif
}
@@ -188,7 +189,8 @@ PETScMatrix::PETScMatrix (const ProcessAdm& padm, const LinSolParams& spar,
}
PETScMatrix::PETScMatrix (const PETScMatrix& B) : adm(B.adm), solParams(B.solParams)
PETScMatrix::PETScMatrix (const PETScMatrix& B) :
nsp(nullptr), adm(B.adm), solParams(B.solParams)
{
MatDuplicate(B.A,MAT_COPY_VALUES,&A);
@@ -196,12 +198,13 @@ PETScMatrix::PETScMatrix (const PETScMatrix& B) : adm(B.adm), solParams(B.solPar
KSPCreate(*adm.getCommunicator(),&ksp);
// Create null space, if any
if (solParams.getNullSpace() == CONSTANT) {
MatNullSpaceCreate(*adm.getCommunicator(),PETSC_TRUE,0,0,&nsp);
if (solParams.getBlock(0).nullspace == CONSTANT) {
nsp = new MatNullSpace;
MatNullSpaceCreate(*adm.getCommunicator(),PETSC_TRUE,0,0,nsp);
#if PETSC_VERSION_MINOR < 6
KSPSetNullSpace(ksp,nsp);
KSPSetNullSpace(ksp,*nsp);
#else
MatSetNullSpace(A,nsp);
MatSetNullSpace(A,*nsp);
#endif
}
LinAlgInit::increfs();
@@ -215,8 +218,10 @@ PETScMatrix::PETScMatrix (const PETScMatrix& B) : adm(B.adm), solParams(B.solPar
PETScMatrix::~PETScMatrix ()
{
// Deallocation of null space
if (solParams.getNullSpace() == CONSTANT)
MatNullSpaceDestroy(&nsp);
if (nsp)
MatNullSpaceDestroy(nsp);
delete nsp;
nsp = nullptr;
// Deallocation of linear solver object.
KSPDestroy(&ksp);
@@ -307,7 +312,7 @@ void PETScMatrix::initAssembly (const SAM& sam, bool)
!strncasecmp(solParams.getPreconditioner(),"ml",2))
sampch->getLocalNodeCoordinates(coords);
else if (solParams.getNullSpace(1) == RIGID_BODY) {
else if (solParams.getBlock(0).nullspace == RIGID_BODY) {
int nsd = sampch->getLocalNodeCoordinates(coords);
#ifdef PARALLEL_PETSC
std::cerr << "WARNING: Rigid body null space not implemented in parallel, ignoring" << std::endl;
@@ -319,22 +324,23 @@ void PETScMatrix::initAssembly (const SAM& sam, bool)
VecSetFromOptions(coordVec);
for (size_t i=0;i<coords.size();++i)
VecSetValue(coordVec, i, coords[i], INSERT_VALUES);
MatNullSpaceCreateRigidBody(coordVec, &nsp);
nsp = new MatNullSpace;
MatNullSpaceCreateRigidBody(coordVec, nsp);
#if PETSC_VERSION_MINOR < 6
KSPSetNullSpace(ksp,nsp);
KSPSetNullSpace(ksp,*nsp);
#else
MatSetNullSpace(A,nsp);
MatSetNullSpace(A,*nsp);
#endif
#endif
}
if (!solParams.dirOrder.empty()) {
if (!solParams.getBlock(0).dirSmoother.empty()) {
dirIndexSet.resize(1);
for (size_t j = 0;j < solParams.dirOrder[0].size();j++) {
for (size_t j = 0;j < solParams.getBlock(0).dirSmoother.size();j++) {
IS permIndex;
if (solParams.dirOrder[0][j] != 123) {
if (solParams.getBlock(0).dirSmoother[j].order != 123) {
PetscIntVec perm;
sampch->getDirOrdering(perm,solParams.dirOrder[0][j]);
sampch->getDirOrdering(perm,solParams.getBlock(0).dirSmoother[j].order);
ISCreateGeneral(*adm.getCommunicator(),perm.size(),&(perm[0]),PETSC_COPY_VALUES,&permIndex);
ISSetPermutation(permIndex);
}
@@ -346,10 +352,10 @@ void PETScMatrix::initAssembly (const SAM& sam, bool)
}
}
int nx = solParams.getLocalPartitioning(0);
int ny = solParams.getLocalPartitioning(1);
int nz = solParams.getLocalPartitioning(2);
int overlap = solParams.getOverlap();
int nx = solParams.getBlock(0).subdomains[0];
int ny = solParams.getBlock(0).subdomains[1];
int nz = solParams.getBlock(0).subdomains[2];
int overlap = solParams.getBlock(0).overlap;
if (nx+ny+nz > 0) {
sampch->getLocalSubdomains(locSubdDofs,nx,ny,nz);
@@ -534,8 +540,8 @@ bool PETScMatrix::solve (const SystemVector& b, SystemVector& x, bool newLHS)
bool PETScMatrix::solve (const Vec& b, Vec& x, bool newLHS, bool knoll)
{
// Reset linear solver
if (nLinSolves && solParams.nResetSolver)
if (nLinSolves%solParams.nResetSolver == 0) {
if (nLinSolves && solParams.getResetSolver())
if (nLinSolves%solParams.getResetSolver() == 0) {
KSPDestroy(&ksp);
KSPCreate(*adm.getCommunicator(),&ksp);
setParams = true;
@@ -564,10 +570,10 @@ bool PETScMatrix::solve (const Vec& b, Vec& x, bool newLHS, bool knoll)
return false;
}
if (solParams.msgLev > 1) {
if (solParams.getMessageLevel() > 1) {
PetscInt its;
KSPGetIterationNumber(ksp,&its);
PetscPrintf(PETSC_COMM_WORLD,"\n Iterations for %s = %D\n",solParams.getMethod(),its);
PetscPrintf(PETSC_COMM_WORLD,"\n Iterations for %s = %D\n",solParams.getMethod().c_str(),its);
}
nLinSolves++;

28
src/LinAlg/PETScMatrix.h
View File

@@ -245,20 +245,20 @@ protected:
//! \brief Solve a linear system
bool solve(const Vec& b, Vec& x, bool newLHS, bool knoll);
Mat A; //!< The actual PETSc matrix
KSP ksp; //!< Linear equation solver
MatNullSpace nsp; //!< Null-space of linear operator
const ProcessAdm& adm; //!< Process administrator
const LinSolParams& solParams; //!< Linear solver parameters
bool setParams; //!< If linear solver parameters are set
IS* elmIS; //!< Element index sets
PetscInt ISsize; //!< Number of index sets/elements
PetscIntMat locSubdDofs; //!< Degrees of freedom for unique subdomains
PetscIntMat subdDofs; //!< Degrees of freedom for subdomains
PetscRealVec coords; //!< Coordinates of local nodes (x0,y0,z0,x1,y1,...)
ISMat dirIndexSet; //!< Direction ordering
int nLinSolves; //!< Number of linear solves
LinAlg::LinearSystemType linsysType; //!< Linear system type
Mat A; //!< The actual PETSc matrix
KSP ksp; //!< Linear equation solver
MatNullSpace* nsp; //!< Null-space of linear operator
const ProcessAdm& adm; //!< Process administrator
const LinSolParams& solParams; //!< Linear solver parameters
bool setParams; //!< If linear solver parameters are set
IS* elmIS; //!< Element index sets
PetscInt ISsize; //!< Number of index sets/elements
PetscIntMat locSubdDofs; //!< Degrees of freedom for unique subdomains
PetscIntMat subdDofs; //!< Degrees of freedom for subdomains
PetscRealVec coords; //!< Coordinates of local nodes (x0,y0,z0,x1,y1,...)
ISMat dirIndexSet; //!< Direction ordering
int nLinSolves; //!< Number of linear solves
LinAlg::LinearSystemType linsysType; //!< Linear system type
#else // dummy implementation when PETSc is not included
virtual SystemMatrix* copy() const { return 0; }

2
src/LinAlg/SystemMatrix.C
View File

@@ -73,7 +73,7 @@ SystemMatrix* SystemMatrix::create (const ProcessAdm& padm, Type matrixType,
#ifdef HAS_PETSC
if (matrixType == PETSC) {
if (spar.getNoBlocks() > 1)
return new PETScBlockMatrix(padm,spar.getComponents(),spar);
return new PETScBlockMatrix(padm,spar);
else
return new PETScMatrix(padm,spar,ltype);
}

107
src/LinAlg/Test/TestLinSolParams.C Normal file
View File

@@ -0,0 +1,107 @@
//==============================================================================
//!
//! \file TestLinSolParams.C
//!
//! \date Nov 2 2015
//!
//! \author Arne Morten Kvarving / SINTEF
//!
//! \brief Tests for LinSolParams
//!
//==============================================================================
#include "LinSolParams.h"
#include "tinyxml.h"
#include "gtest/gtest.h"
TEST(TestLinSolParams, ParseBlocks)
{
TiXmlDocument doc;
doc.LoadFile("src/LinAlg/Test/refdata/linsolver_blocks.xml");
LinSolParams params(2);
params.read(doc.RootElement());
// global settings
ASSERT_EQ(params.getDimension(), 2u);
ASSERT_EQ(params.getNoBlocks(), 2u);
ASSERT_STREQ(params.getPreconditioner(), "fieldsplit");
ASSERT_EQ(params.getMethod(), "bcgs");
ASSERT_FLOAT_EQ(params.getAbsTolerance(), 1e-12);
ASSERT_FLOAT_EQ(params.getRelTolerance(), 1e-4);
ASSERT_FLOAT_EQ(params.getDivTolerance(), 1e6);
ASSERT_EQ(params.getMaxIterations(), 123);
// first block
ASSERT_EQ(params.getBlock(0).basis, 1u);
ASSERT_EQ(params.getBlock(0).comps, 12u);
ASSERT_EQ(params.getBlock(0).prec, "sor");
// second block
ASSERT_EQ(params.getBlock(1).basis, 2u);
ASSERT_EQ(params.getBlock(1).comps, 3u);
ASSERT_EQ(params.getBlock(1).prec, "ilu");
ASSERT_EQ(params.getBlock(1).ilu_fill_level, 1);
}
TEST(TestLinSolParams, ParseGAMG)
{
TiXmlDocument doc;
doc.LoadFile("src/LinAlg/Test/refdata/linsolver_gamg.xml");
LinSolParams params(2);
params.read(doc.RootElement());
ASSERT_EQ(params.getBlock(0).gamg.type, "agg");
ASSERT_EQ(params.getBlock(0).gamg.procEqLimit, 1000);
ASSERT_EQ(params.getBlock(0).gamg.repartition, 1);
ASSERT_EQ(params.getBlock(0).gamg.useAggGasm, 1);
ASSERT_EQ(params.getBlock(0).gamg.reuseInterp, 1);
ASSERT_FLOAT_EQ(params.getBlock(0).gamg.threshold, 0.1);
ASSERT_EQ(params.getBlock(0).finesmoother, "compositedir");
ASSERT_EQ(params.getBlock(0).dirSmoother.size(), 1u);
ASSERT_EQ(params.getBlock(0).dirSmoother[0].type, "ilu");
ASSERT_EQ(params.getBlock(0).dirSmoother[0].order, 12);
}
TEST(TestLinSolParams, ParseML)
{
TiXmlDocument doc;
doc.LoadFile("src/LinAlg/Test/refdata/linsolver_ml.xml");
LinSolParams params(2);
params.read(doc.RootElement());
ASSERT_EQ(params.getBlock(0).ml.coarsePackage, "petsc");
ASSERT_EQ(params.getBlock(0).ml.coarseSolver, "lu");
ASSERT_EQ(params.getBlock(0).ml.coarsenScheme, "agg");
ASSERT_EQ(params.getBlock(0).ml.symmetrize, 1);
ASSERT_EQ(params.getBlock(0).ml.blockScaling, 1);
ASSERT_EQ(params.getBlock(0).ml.putOnSingleProc, 1000);
ASSERT_EQ(params.getBlock(0).ml.reuseInterp, 1);
ASSERT_EQ(params.getBlock(0).ml.reusable, 1);
ASSERT_EQ(params.getBlock(0).ml.keepAggInfo, 1);
ASSERT_EQ(params.getBlock(0).ml.aux, 1);
ASSERT_FLOAT_EQ(params.getBlock(0).ml.threshold, 0.1);
ASSERT_FLOAT_EQ(params.getBlock(0).ml.auxThreshold, 0.05);
}
TEST(TestLinSolParams, ParseHypre)
{
TiXmlDocument doc;
doc.LoadFile("src/LinAlg/Test/refdata/linsolver_hypre.xml");
LinSolParams params(2);
params.read(doc.RootElement());
ASSERT_EQ(params.getBlock(0).hypre.type, "boom");
ASSERT_EQ(params.getBlock(0).hypre.noAggCoarse, 1);
ASSERT_EQ(params.getBlock(0).hypre.noPathAggCoarse, 2);
ASSERT_FLOAT_EQ(params.getBlock(0).hypre.truncation, 0.005);
ASSERT_FLOAT_EQ(params.getBlock(0).hypre.threshold, 0.01);
ASSERT_EQ(params.getBlock(0).hypre.coarsenScheme, "agg");
}

16
src/LinAlg/Test/refdata/linsolver_blocks.xml Normal file
View File

@@ -0,0 +1,16 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<linearsolver>
<type>bcgs</type>
<pc>fieldsplit</pc>
<block basis="1" components="12">
<pc>sor</pc>
</block>
<block basis="2" components="3">
<pc>ilu</pc>
<ilu_fill_level>1</ilu_fill_level>
</block>
<atol>1.0e-12</atol>
<rtol>1.0e-4</rtol>
<dtol>1.0e6</dtol>
<maxits>123</maxits>
</linearsolver>

13
src/LinAlg/Test/refdata/linsolver_gamg.xml Normal file
View File

@@ -0,0 +1,13 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<linearsolver>
<gamg>
<type>agg</type>
<repartition>1</repartition>
<use_agg_smoother>1</use_agg_smoother>
<proc_eq_limit>1000</proc_eq_limit>
<reuse_interpolation>1</reuse_interpolation>
<threshold>0.1</threshold>
</gamg>
<finesmoother>compositedir</finesmoother>
<dirsmoother type="ilu" order="12"/>
</linearsolver>

11
src/LinAlg/Test/refdata/linsolver_hypre.xml Normal file
View File

@@ -0,0 +1,11 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<linearsolver>
<hypre>
<type>boom</type>
<no_agg_coarse>1</no_agg_coarse>
<no_path_agg_coarse>2</no_path_agg_coarse>
<truncation>0.005</truncation>
<threshold>0.01</threshold>
<coarsen_scheme>agg</coarsen_scheme>
</hypre>
</linearsolver>

19
src/LinAlg/Test/refdata/linsolver_ml.xml Normal file
View File

@@ -0,0 +1,19 @@
<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
<linearsolver>
<ml>
<coarse_package>petsc</coarse_package>
<coarse_solver>lu</coarse_solver>
<coarsen_scheme>agg</coarsen_scheme>
<symmetrize>1</symmetrize>
<repartition>1</repartition>
<block_scaling>1</block_scaling>
<put_on_single_proc>1000</put_on_single_proc>
<reuse_interpolation>1</reuse_interpolation>
<reusable>1</reusable>
<keep_agg_info>1</keep_agg_info>
<damping_factor>0.9</damping_factor>
<aux>1</aux>
<threshold>0.1</threshold>
<aux_threshold>0.05</aux_threshold>
</ml>
</linearsolver>

16
src/SIM/SIMbase.C
View File

@@ -491,15 +491,21 @@ bool SIMbase::parse (const TiXmlElement* elem)
myModel.resize(1,this->createDefaultGeometry(elem));
}
if (!strcasecmp(elem->Value(),"linearsolver")) {
if (!mySolParams)
mySolParams = new LinSolParams(nsd);
result &= mySolParams->read(elem);
}
const TiXmlElement* child = elem->FirstChildElement();
for (; child; child = child->NextSiblingElement())
if (!strcasecmp(elem->Value(),"geometry"))
result &= this->parseGeometryTag(child);
else if (!strcasecmp(elem->Value(),"boundaryconditions"))
result &= this->parseBCTag(child);
else if (!strcasecmp(elem->Value(),"linearsolver"))
else if (!strcasecmp(elem->Value(),"linearsolver")) {
result &= this->parseLinSolTag(child);
else if (!strcasecmp(elem->Value(),"eigensolver"))
} else if (!strcasecmp(elem->Value(),"eigensolver"))
result &= opt.parseEigSolTag(child);
else if (!strcasecmp(elem->Value(),"postprocessing"))
result &= this->parseOutputTag(child);
@@ -764,13 +770,9 @@ bool SIMbase::parseLinSolTag (const TiXmlElement* elem)
if (!strcasecmp(elem->Value(),"class")) {
if (elem->FirstChild())
opt.setLinearSolver(elem->FirstChild()->Value());
return true;
}
if (!mySolParams)
mySolParams = new LinSolParams(nsd);
return mySolParams->read(elem);
return true;
}