PetscSolverPfc.cc 6.32 KB
 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 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 `````` } void PetscSolverPfc::initSolver(KSP &ksp) { // Create FGMRES based outer solver KSPCreate(meshDistributor->getMpiComm(0), &ksp); KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); 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 104 105 106 `````` double delta = sqrt(*tau); double eps = sqrt(0.5); `````` Praetorius, Simon committed Oct 23, 2013 107 108 109 110 111 112 113 114 115 116 117 118 119 120 `````` 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 121 `````` // create MpK-matrix `````` Praetorius, Simon committed Oct 23, 2013 122 123 `````` DOFMatrix matrixMpK(feSpace, feSpace); Operator laplaceOp2(feSpace, feSpace); `````` Praetorius, Simon committed Nov 27, 2013 124 `````` Simple_SOT sot2(delta); `````` Praetorius, Simon committed Oct 23, 2013 125 126 127 128 129 130 131 132 133 `````` 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 134 135 136 137 138 139 140 141 142 143 144 145 146 147 `````` // create MpK2-matrix DOFMatrix matrixMpK2(feSpace, feSpace); Operator laplaceOp3(feSpace, feSpace); Simple_SOT sot3(eps*sqrt(delta)); 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 148 149 150 151 152 153 154 155 156 157 158 159 `````` // 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 160 `````` data.delta = delta; `````` Praetorius, Simon committed Oct 23, 2013 161 162 `````` data.tau = (*tau); `````` Praetorius, Simon committed Nov 27, 2013 163 164 165 `````` petsc_helper::setSolver(data.kspM, "M_", KSPCG, PCBJACOBI, 0.0, 1e-8, 5); petsc_helper::setSolver(data.kspMpK, "MpK_", KSPCG, PCBJACOBI, 0.0, 1e-6, 20); petsc_helper::setSolver(data.kspMpK2, "MpK2_", KSPCG, PCBJACOBI, 0.0, 1e-6, 10); `````` Praetorius, Simon committed Oct 23, 2013 166 167 168 169 170 171 172 `````` } void PetscSolverPfc::exitPreconditioner(PC pc) { FUNCNAME("PetscSolverPfc::exit()"); `````` Praetorius, Simon committed Nov 27, 2013 173 174 ``````// MatDestroy(&matMpK); // MatDestroy(&matMpK2); `````` Praetorius, Simon committed Oct 23, 2013 175 176 177 178 `````` solverM->destroyMatrixData(); solverK->destroyMatrixData(); solverMpK->destroyMatrixData(); `````` Praetorius, Simon committed Nov 27, 2013 179 `````` solverMpK2->destroyMatrixData(); `````` Praetorius, Simon committed Oct 23, 2013 180 181 182 183 `````` solverM->destroyVectorData(); solverK->destroyVectorData(); solverMpK->destroyVectorData(); `````` Praetorius, Simon committed Nov 27, 2013 184 `````` solverMpK2->destroyVectorData(); `````` Praetorius, Simon committed Oct 23, 2013 185 186 187 188 189 190 191 `````` delete solverM; solverM = NULL; delete solverK; solverK = NULL; delete solverMpK; solverMpK = NULL; `````` Praetorius, Simon committed Nov 27, 2013 192 193 `````` delete solverMpK2; solverMpK2 = NULL; `````` Praetorius, Simon committed Oct 23, 2013 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 `````` } 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; } } }``````