/****************************************************************************** * * Extension of AMDiS - Adaptive multidimensional simulations * * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved. * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis * * Authors: Simon Praetorius et al. * * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. * * * See also license.opensource.txt in the distribution. * ******************************************************************************/ #include "PetscSolverPfc.h" #include "parallel/PetscHelper.h" #include "parallel/PetscSolverGlobalMatrix.h" #include "Expressions.h" namespace AMDiS { namespace Parallel { using namespace std; /// solve Pfc Preconditioner PetscErrorCode pcPfcShell(PC pc, Vec b, Vec x) // solve Px=b { FUNCNAME("pcPfcShell()"); void *ctx; PCShellGetContext(pc, &ctx); PfcData* data = static_cast(ctx); Vec b1, b2, b3, x1, x2, x3; VecNestGetSubVec(b, 0, &b1); VecNestGetSubVec(b, 1, &b2); VecNestGetSubVec(b, 2, &b3); VecNestGetSubVec(x, 0, &x1); VecNestGetSubVec(x, 1, &x2); VecNestGetSubVec(x, 2, &x3); // Hilfsvariablen Vec y1, y2, tmp; VecDuplicate(b1, &y1); VecDuplicate(b1, &y2); VecDuplicate(b1, &tmp); KSPSolve(data->kspM, b1, y1); // M*y1 = b1 MatMult(data->matK, y1, tmp); // tmp := K*y1 VecAYPX(tmp, -sqr(data->delta), b2); // tmp := b2 - M0*tau*tmp KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp MatMult(data->matM, y2, tmp); // tmp := M*y2 KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp MatMult(data->matM, x2, tmp); // tmp := M*x2 KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp VecCopy(x2, x1); // x1 := x2 VecAXPBYPCZ(x1, 1.0, 1.0/(data->delta), -1.0/(data->delta), y1, y2); // x1 = 1*y1 + factor*y2 - factor*x1 MatMult(data->matK, x2, tmp); // tmp := K*x2 VecAYPX(tmp, -1.0, b3); // tmp := b3 - tmp KSPSolve(data->kspM, tmp, x3); // M*x3 = tmp VecDestroy(&y1); VecDestroy(&y2); VecDestroy(&tmp); return 0; } void PetscSolverPfc::initSolver(KSP &ksp) { // Create FGMRES based outer solver KSPCreate(meshDistributor->getMpiComm(0), &ksp); petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior()); if (getInfo() >= 10) petsc::ksp_monitor_set(ksp, KSPMonitorDefault); else if (getInfo() >= 20) petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm); petsc_helper::setSolver(ksp, "pfc_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations()); KSPSetFromOptions(ksp); if (useOldInitialGuess) KSPSetInitialGuessNonzero(ksp, PETSC_TRUE); } void PetscSolverPfc::initPreconditioner(PC pc) { FUNCNAME("PetscSolverPfc::initPreconditioner()"); TEST_EXIT(tau)("tau pointer not set!\n"); PCSetType(getPc(), PCSHELL); PCShellSetApply(getPc(), pcPfcShell); PCShellSetContext(getPc(), &data); double delta = std::sqrt((*tau) * M0); const FiniteElemSpace *feSpace = componentSpaces[0]; // create mass-matrix DOFMatrix matrixM(feSpace, feSpace); Operator opM(feSpace, feSpace); addZOT(opM, 1.0); matrixM.assembleOperator(opM); solverM = createSubSolver(0, "M_"); solverM->fillPetscMatrix(&matrixM); data.matM = solverM->getMatInterior(); data.kspM = solverM->getSolver(); // create MpK-matrix DOFMatrix matrixMpK(feSpace, feSpace); Operator opMpK(feSpace, feSpace); addZOT(opMpK, 1.0); addSOT(opMpK, delta); matrixMpK.assembleOperator(opMpK); solverMpK = createSubSolver(0, "MpK_"); solverMpK->fillPetscMatrix(&matrixMpK); matMpK = solverMpK->getMatInterior(); data.kspMpK = solverMpK->getSolver(); // create MpK2-matrix DOFMatrix matrixMpK2(feSpace, feSpace); Operator opMpK2(feSpace, feSpace); addZOT(opMpK2, 1.0); addSOT(opMpK2, std::sqrt(delta)); matrixMpK2.assembleOperator(opMpK2); solverMpK2 = createSubSolver(0, "MpK2_"); solverMpK2->fillPetscMatrix(&matrixMpK2); matMpK2 = solverMpK2->getMatInterior(); data.kspMpK2 = solverMpK2->getSolver(); // create laplace-matrix DOFMatrix matrixK(feSpace, feSpace); Operator opK(feSpace, feSpace); addSOT(opK, 1.0); matrixK.assembleOperator(opK); solverK = createSubSolver(0, "K_"); solverK->fillPetscMatrix(&matrixK); data.matK = solverK->getMatInterior(); // === Setup preconditioner data === data.delta = delta; data.tau = (*tau); int nIterM=5, nIterMpL=20, nIterMpL2=10; double tolM=PETSC_DEFAULT, tolMpL=PETSC_DEFAULT, tolMpL2=PETSC_DEFAULT; double rtolM=PETSC_DEFAULT, rtolMpL=PETSC_DEFAULT, rtolMpL2=PETSC_DEFAULT; bool directM = false, directMpL = false, directMpL2 = false; Parameters::get("precon_pfc_M->max iteration", nIterM); Parameters::get("precon_pfc_MpL->max iteration", nIterMpL); Parameters::get("precon_pfc_MpL2->max iteration", nIterMpL2); Parameters::get("precon_pfc_M->tolerance", tolM); Parameters::get("precon_pfc_MpL->tolerance", tolMpL); Parameters::get("precon_pfc_MpL2->tolerance", tolMpL2); Parameters::get("precon_pfc_M->relative tolerance", rtolM); Parameters::get("precon_pfc_MpL->relative tolerance", rtolMpL); Parameters::get("precon_pfc_MpL2->relative tolerance", rtolMpL2); Parameters::get("precon_pfc_M->use direct solver", directM); Parameters::get("precon_pfc_MpL->use direct solver", directMpL); Parameters::get("precon_pfc_MpL2->use direct solver", directMpL2); bool useAMG = false; Parameters::get("precon_pfc_MpL2->use AMG", useAMG); if (directM) petsc_helper::setSolverWithLu(data.kspM, "M_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1); else petsc_helper::setSolver(data.kspM, "M_", KSPCG, PCBJACOBI, rtolM, tolM, nIterM); if (directMpL) petsc_helper::setSolverWithLu(data.kspMpK, "MpL_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1); else petsc_helper::setSolver(data.kspMpK, "MpL_", KSPCG, PCBJACOBI, rtolMpL, tolMpL, nIterMpL); if (directMpL2) petsc_helper::setSolverWithLu(data.kspMpK2, "MpL2_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1); else if (useAMG) petsc_helper::setSolver(data.kspMpK2, "MpL2_", KSPRICHARDSON, PCHYPRE, rtolMpL2, tolMpL2, nIterMpL2); else petsc_helper::setSolver(data.kspMpK2, "MpL2_", KSPCG, PCBJACOBI, rtolMpL2, tolMpL2, nIterMpL2); } void PetscSolverPfc::exitPreconditioner(PC pc) { FUNCNAME("PetscSolverPfc::exit()"); // MatDestroy(&matMpK); // MatDestroy(&matMpK2); solverM->destroyMatrixData(); solverK->destroyMatrixData(); solverMpK->destroyMatrixData(); solverMpK2->destroyMatrixData(); solverM->destroyVectorData(); solverK->destroyVectorData(); solverMpK->destroyVectorData(); solverMpK2->destroyVectorData(); delete solverM; solverM = NULL; delete solverK; solverK = NULL; delete solverMpK; solverMpK = NULL; delete solverMpK2; solverMpK2 = NULL; } PetscSolver* PetscSolverPfc::createSubSolver(int component, string kspPrefix) { FUNCNAME("PetscSolverPfc::createSubSolver()"); vector fe; fe.push_back(componentSpaces[component]); PetscSolver* subSolver = new PetscSolverGlobalMatrix(""); subSolver->setKspPrefix(kspPrefix.c_str()); subSolver->setMeshDistributor(meshDistributor, 0); subSolver->init(fe, fe); ParallelDofMapping &subDofMap = subSolver->getDofMap(); subDofMap[0] = dofMap[component]; subDofMap.update(); return subSolver; } } }