PetscSolverPfc.cc 7.4 KB
Newer Older
 Praetorius, Simon committed Oct 23, 2013 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 ``````/****************************************************************************** * * 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" 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, -(data->tau), b2); // tmp := b2 - tau*tmp KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp `````` Praetorius, Simon committed Nov 27, 2013 55 `````` `````` Praetorius, Simon committed Oct 23, 2013 56 `````` MatMult(data->matM, y2, tmp); // tmp := M*y2 `````` Praetorius, Simon committed Nov 27, 2013 57 58 59 `````` KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp MatMult(data->matM, x2, tmp); // tmp := M*x2 KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp `````` Praetorius, Simon committed Oct 23, 2013 60 61 62 63 64 65 `````` 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 `````` Praetorius, Simon committed Nov 27, 2013 66 67 `````` KSPSolve(data->kspM, tmp, x3); // M*x3 = tmp `````` Praetorius, Simon committed Oct 23, 2013 68 69 70 71 72 `````` VecDestroy(&y1); VecDestroy(&y2); VecDestroy(&tmp); `````` Praetorius, Simon committed Nov 27, 2013 73 `````` return 0; `````` Praetorius, Simon committed Oct 23, 2013 74 75 76 77 78 79 `````` } void PetscSolverPfc::initSolver(KSP &ksp) { // Create FGMRES based outer solver KSPCreate(meshDistributor->getMpiComm(0), &ksp); `````` Praetorius, Simon committed Nov 11, 2014 80 81 82 ``````#if (PETSC_VERSION_MINOR >= 5) KSPSetOperators(ksp, getMatInterior(), getMatInterior()); #else `````` Praetorius, Simon committed Oct 23, 2013 83 `````` KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); `````` Praetorius, Simon committed Nov 11, 2014 84 ``````#endif `````` Praetorius, Simon committed Oct 23, 2013 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 `````` if (getInfo() >= 10) KSPMonitorSet(ksp, KSPMonitorDefault, PETSC_NULL, PETSC_NULL); else if (getInfo() >= 20) KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL); 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); `````` Praetorius, Simon committed Nov 27, 2013 108 109 `````` double delta = sqrt(*tau); `````` Praetorius, Simon committed Oct 23, 2013 110 111 112 113 114 115 116 117 118 119 120 121 122 123 `````` const FiniteElemSpace *feSpace = componentSpaces[0]; // create mass-matrix DOFMatrix matrixM(feSpace, feSpace); Operator massOp(feSpace, feSpace); Simple_ZOT zot; massOp.addTerm(&zot); matrixM.assembleOperator(massOp); solverM = createSubSolver(0, "M_"); solverM->fillPetscMatrix(&matrixM); data.matM = solverM->getMatInterior(); data.kspM = solverM->getSolver(); `````` Praetorius, Simon committed Nov 27, 2013 124 `````` // create MpK-matrix `````` Praetorius, Simon committed Oct 23, 2013 125 126 `````` DOFMatrix matrixMpK(feSpace, feSpace); Operator laplaceOp2(feSpace, feSpace); `````` Praetorius, Simon committed Nov 27, 2013 127 `````` Simple_SOT sot2(delta); `````` Praetorius, Simon committed Oct 23, 2013 128 129 130 131 132 133 134 135 136 `````` laplaceOp2.addTerm(&zot); laplaceOp2.addTerm(&sot2); matrixMpK.assembleOperator(laplaceOp2); solverMpK = createSubSolver(0, "MpK_"); solverMpK->fillPetscMatrix(&matrixMpK); matMpK = solverMpK->getMatInterior(); data.kspMpK = solverMpK->getSolver(); `````` Praetorius, Simon committed Nov 27, 2013 137 138 139 140 `````` // create MpK2-matrix DOFMatrix matrixMpK2(feSpace, feSpace); Operator laplaceOp3(feSpace, feSpace); `````` Praetorius, Simon committed Aug 06, 2014 141 `````` Simple_SOT sot3(sqrt(delta)); `````` Praetorius, Simon committed Nov 27, 2013 142 143 144 145 146 147 148 149 150 `````` laplaceOp3.addTerm(&zot); laplaceOp3.addTerm(&sot3); matrixMpK2.assembleOperator(laplaceOp3); solverMpK2 = createSubSolver(0, "MpK2_"); solverMpK2->fillPetscMatrix(&matrixMpK2); matMpK2 = solverMpK2->getMatInterior(); data.kspMpK2 = solverMpK2->getSolver(); `````` Praetorius, Simon committed Oct 23, 2013 151 152 153 154 155 156 157 158 159 160 161 162 `````` // create laplace-matrix DOFMatrix matrixK(feSpace, feSpace); Operator laplaceOp(feSpace, feSpace); Simple_SOT sot; laplaceOp.addTerm(&sot); matrixK.assembleOperator(laplaceOp); solverK = createSubSolver(0, "K_"); solverK->fillPetscMatrix(&matrixK); data.matK = solverK->getMatInterior(); // === Setup preconditioner data === `````` Praetorius, Simon committed Nov 27, 2013 163 `````` data.delta = delta; `````` Praetorius, Simon committed Oct 23, 2013 164 `````` data.tau = (*tau); `````` Praetorius, Simon committed Aug 06, 2014 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 `````` 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; 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); bool useAMG = false; Parameters::get("precon_pfc_MpL2->use AMG", useAMG); `````` Praetorius, Simon committed Oct 23, 2013 181 `````` `````` Praetorius, Simon committed Aug 06, 2014 182 183 184 185 186 187 188 `````` petsc_helper::setSolver(data.kspM, "M_", KSPCG, PCJACOBI, rtolM, tolM, nIterM); petsc_helper::setSolver(data.kspMpK, "MpK_", KSPCG, PCJACOBI, rtolMpL, tolMpL, nIterMpL); if (!useAMG) { petsc_helper::setSolver(data.kspMpK2, "MpK2_", KSPCG, PCJACOBI, rtolMpL2, tolMpL2, nIterMpL2); } else { petsc_helper::setSolver(data.kspMpK2, "MpK2_", KSPRICHARDSON, PCHYPRE, rtolMpL2, tolMpL2, nIterMpL2); } `````` Praetorius, Simon committed Oct 23, 2013 189 190 191 192 193 194 195 `````` } void PetscSolverPfc::exitPreconditioner(PC pc) { FUNCNAME("PetscSolverPfc::exit()"); `````` Praetorius, Simon committed Nov 27, 2013 196 197 ``````// MatDestroy(&matMpK); // MatDestroy(&matMpK2); `````` Praetorius, Simon committed Oct 23, 2013 198 199 200 201 `````` solverM->destroyMatrixData(); solverK->destroyMatrixData(); solverMpK->destroyMatrixData(); `````` Praetorius, Simon committed Nov 27, 2013 202 `````` solverMpK2->destroyMatrixData(); `````` Praetorius, Simon committed Oct 23, 2013 203 204 205 206 `````` solverM->destroyVectorData(); solverK->destroyVectorData(); solverMpK->destroyVectorData(); `````` Praetorius, Simon committed Nov 27, 2013 207 `````` solverMpK2->destroyVectorData(); `````` Praetorius, Simon committed Oct 23, 2013 208 209 210 211 212 213 214 `````` delete solverM; solverM = NULL; delete solverK; solverK = NULL; delete solverMpK; solverMpK = NULL; `````` Praetorius, Simon committed Nov 27, 2013 215 216 `````` delete solverMpK2; solverMpK2 = NULL; `````` Praetorius, Simon committed Oct 23, 2013 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 `````` } 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; } } }``````