PetscSolverPfc.cc 6.32 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
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<PfcData*>(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's avatar
Praetorius, Simon committed
55
    
Praetorius, Simon's avatar
Praetorius, Simon committed
56
    MatMult(data->matM, y2, tmp);			// tmp := M*y2
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
66 67
    
    KSPSolve(data->kspM, tmp, x3); 			// M*x3 = tmp    
Praetorius, Simon's avatar
Praetorius, Simon committed
68 69 70 71 72

    VecDestroy(&y1);
    VecDestroy(&y2);
    VecDestroy(&tmp);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
73
    return 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
104 105 106
    double delta = sqrt(*tau);
    double eps = sqrt(0.5);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
121
    // create MpK-matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
122 123
    DOFMatrix matrixMpK(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
Praetorius, Simon's avatar
Praetorius, Simon committed
124
    Simple_SOT sot2(delta);
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
160
    data.delta = delta;
Praetorius, Simon's avatar
Praetorius, Simon committed
161 162
    data.tau = (*tau);

Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
166 167 168 169 170 171 172
  }


  void PetscSolverPfc::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverPfc::exit()");
    
Praetorius, Simon's avatar
Praetorius, Simon committed
173 174
//     MatDestroy(&matMpK);
//     MatDestroy(&matMpK2);
Praetorius, Simon's avatar
Praetorius, Simon committed
175 176 177 178
    
    solverM->destroyMatrixData();
    solverK->destroyMatrixData();
    solverMpK->destroyMatrixData();
Praetorius, Simon's avatar
Praetorius, Simon committed
179
    solverMpK2->destroyMatrixData();
Praetorius, Simon's avatar
Praetorius, Simon committed
180 181 182 183
    
    solverM->destroyVectorData();
    solverK->destroyVectorData();
    solverMpK->destroyVectorData();
Praetorius, Simon's avatar
Praetorius, Simon committed
184
    solverMpK2->destroyVectorData();
Praetorius, Simon's avatar
Praetorius, Simon committed
185 186 187 188 189 190 191
    
    delete solverM;
    solverM = NULL;
    delete solverK;
    solverK = NULL;
    delete solverMpK;
    solverMpK = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
192 193
    delete solverMpK2;
    solverMpK2 = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
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<const FiniteElemSpace*> 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;
  }
} }