PetscSolverPfc.cc 7.4 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
  }
  
  void PetscSolverPfc::initSolver(KSP &ksp)
  {
    // Create FGMRES based outer solver
    KSPCreate(meshDistributor->getMpiComm(0), &ksp);
80 81 82
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(ksp, getMatInterior(), getMatInterior());
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
83
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
84
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
108 109
    double delta = sqrt(*tau);
    
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
124
    // create MpK-matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
125 126
    DOFMatrix matrixMpK(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
Praetorius, Simon's avatar
Praetorius, Simon committed
127
    Simple_SOT sot2(delta);
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
137 138 139 140
    
    // create MpK2-matrix
    DOFMatrix matrixMpK2(feSpace, feSpace);
    Operator laplaceOp3(feSpace, feSpace);
141
    Simple_SOT sot3(sqrt(delta));
Praetorius, Simon's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
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's avatar
Praetorius, Simon committed
163
    data.delta = delta;
Praetorius, Simon's avatar
Praetorius, Simon committed
164
    data.tau = (*tau);
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's avatar
Praetorius, Simon committed
181

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's avatar
Praetorius, Simon committed
189 190 191 192 193 194 195
  }


  void PetscSolverPfc::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverPfc::exit()");
    
Praetorius, Simon's avatar
Praetorius, Simon committed
196 197
//     MatDestroy(&matMpK);
//     MatDestroy(&matMpK2);
Praetorius, Simon's avatar
Praetorius, Simon committed
198 199 200 201
    
    solverM->destroyMatrixData();
    solverK->destroyMatrixData();
    solverMpK->destroyMatrixData();
Praetorius, Simon's avatar
Praetorius, Simon committed
202
    solverMpK2->destroyMatrixData();
Praetorius, Simon's avatar
Praetorius, Simon committed
203 204 205 206
    
    solverM->destroyVectorData();
    solverK->destroyVectorData();
    solverMpK->destroyVectorData();
Praetorius, Simon's avatar
Praetorius, Simon committed
207
    solverMpK2->destroyVectorData();
Praetorius, Simon's avatar
Praetorius, Simon committed
208 209 210 211 212 213 214
    
    delete solverM;
    solverM = NULL;
    delete solverK;
    solverK = NULL;
    delete solverMpK;
    solverMpK = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
215 216
    delete solverMpK2;
    solverMpK2 = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
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<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;
  }
} }