PetscSolverPfc_diag.cc 7.18 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14
/******************************************************************************
 *
 * 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.
15
 *
16 17 18 19 20 21 22 23 24
 ******************************************************************************/

#include "PetscSolverPfc_diag.h"
#include "parallel/PetscHelper.h"
#include "parallel/PetscSolverGlobalMatrix.h"

namespace AMDiS { namespace Parallel {

  using namespace std;
25

26 27 28
  /// solve Pfc Preconditioner
  PetscErrorCode pcPfcDiagShell(PC pc, Vec b, Vec x) // solve Px=b
  { FUNCNAME("pcPfcDiagShell()");
29

30 31 32
    void *ctx;
    PCShellGetContext(pc, &ctx);
    PfcDiagData* data = static_cast<PfcDiagData*>(ctx);
33

34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
    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);

50
    KSPSolve(data->kspM, b1, y1); 			// M*y1 = b1
51 52
    MatMult(data->matK, y1, tmp); 			// tmp := K*y1
    VecAYPX(tmp, -(data->tau), b2);			// tmp := b2 - tau*tmp
53

54
    KSPSolve(data->kspMpK, tmp, y2); 			// (M+sqrt(tau)K)*y2 = tmp
55 56


57 58 59 60 61 62 63 64
    MatMult(data->matM, y2, tmp);			// tmp := M*y2
    KSPSolve(data->kspS, tmp, x2);

    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
65 66

    KSPSolve(data->kspM, tmp, x3); 			// M*x3 = tmp
67 68 69 70

    VecDestroy(&y1);
    VecDestroy(&y2);
    VecDestroy(&tmp);
71

72 73
    return 0;
  }
74

75 76 77 78
  void PetscSolverPfcDiag::initSolver(KSP &ksp)
  {
    // Create FGMRES based outer solver
    KSPCreate(meshDistributor->getMpiComm(0), &ksp);
79 80 81
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(ksp, getMatInterior(), getMatInterior());
#else
82
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
83
#endif
84 85 86 87 88 89
    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);
90 91


92 93 94
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }
95

96 97 98 99 100 101

  void PetscSolverPfcDiag::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverPfcDiag::initPreconditioner()");

    TEST_EXIT(tau)("tau pointer not set!\n");
102

103 104 105
    PCSetType(getPc(), PCSHELL);
    PCShellSetApply(getPc(), pcPfcDiagShell);
    PCShellSetContext(getPc(), &data);
106

107
    double delta = sqrt(*tau);
108

109
    const FiniteElemSpace *feSpace = componentSpaces[0];
110

111 112 113 114 115 116
    // create mass-matrix
    DOFMatrix matrixM(feSpace, feSpace);
    Operator massOp(feSpace, feSpace);
    Simple_ZOT zot;
    massOp.addTerm(&zot);
    matrixM.assembleOperator(massOp);
117

118 119 120 121
    solverM = createSubSolver(0, "M_");
    solverM->fillPetscMatrix(&matrixM);
    data.matM = solverM->getMatInterior();
    data.kspM = solverM->getSolver();
122

123 124 125 126 127 128 129
    // create MpK-matrix
    DOFMatrix matrixMpK(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
    Simple_SOT sot2(delta);
    laplaceOp2.addTerm(&zot);
    laplaceOp2.addTerm(&sot2);
    matrixMpK.assembleOperator(laplaceOp2);
130

131 132 133 134
    solverMpK = createSubSolver(0, "MpK_");
    solverMpK->fillPetscMatrix(&matrixMpK);
    matMpK = solverMpK->getMatInterior();
    data.kspMpK = solverMpK->getSolver();
135 136


137
    // create MpK2-matrix
138

139 140 141 142 143 144
    // create laplace-matrix
    DOFMatrix matrixK(feSpace, feSpace);
    Operator laplaceOp(feSpace, feSpace);
    Simple_SOT sot;
    laplaceOp.addTerm(&sot);
    matrixK.assembleOperator(laplaceOp);
145

146 147 148
    solverK = createSubSolver(0, "K_");
    solverK->fillPetscMatrix(&matrixK);
    data.matK = solverK->getMatInterior();
149

150 151 152 153 154
    // create Matrix S
    solverS = createSubSolver(0, "MpK2_");
    solverS->fillPetscMatrix(&matrixM); // just to create a dummy matrix
    matS = solverS->getMatInterior();
    data.kspS = solverS->getSolver();
155

156 157 158 159
    Vec x; // need to initialize vector
    solverS->createVec(*solverS->getDofMapping(), x);
    Mat DK;
    MatDuplicate(data.matK, MAT_COPY_VALUES, &DK);
160

161 162 163 164 165 166 167 168
    MatGetDiagonal(data.matM, x);
    VecReciprocal(x);
    MatDiagonalScale(DK, x, PETSC_NULL); 				// DK := M_D^(-1)*K
    MatMatMult(data.matK, DK, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &matS); 	// S := K*DK
    MatAYPX(matS, delta, data.matM, DIFFERENT_NONZERO_PATTERN); 			// S = delta*S + M
    MatAXPY(matS, (-2.0*delta), data.matK, DIFFERENT_NONZERO_PATTERN);	// S = S - 2*sqrt(delta)*K
    VecDestroy(&x);
    MatDestroy(&DK);
169

170 171 172
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(data.kspS, matS, matS);
#else
173
    KSPSetOperators(data.kspS, matS, matS, SAME_NONZERO_PATTERN);
174
#endif
175 176

    // === Setup preconditioner data ===
177 178
    data.delta = delta;
    data.tau = (*tau);
179

180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197
    int nIterM=5, nIterMpL=20, nIterMpL2=10;
    double tolM=1.e-6, tolMpL=1.e-6, tolMpL2=1.e-6;
    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);

    petsc_helper::setSolver(data.kspM, "M_", KSPCG, PCBJACOBI, 0.0, tolM, nIterM);
    petsc_helper::setSolver(data.kspMpK, "MpK_", KSPCG, PCBJACOBI, 0.0, tolMpL, nIterMpL);
    petsc_helper::setSolver(data.kspS, "MpK2_", KSPCG, PCBJACOBI, 0.0, tolMpL2, nIterMpL2);
  }


  void PetscSolverPfcDiag::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverPfcDiag::exit()");
198

199 200
//     MatDestroy(&matMpK);
//     MatDestroy(&matMpK2);
201

202 203 204 205
    solverM->destroyMatrixData();
    solverK->destroyMatrixData();
    solverMpK->destroyMatrixData();
    solverS->destroyMatrixData();
206

207 208 209 210
    solverM->destroyVectorData();
    solverK->destroyVectorData();
    solverMpK->destroyVectorData();
    solverS->destroyVectorData();
211

212 213 214 215 216 217 218 219 220
    delete solverM;
    solverM = NULL;
    delete solverK;
    solverK = NULL;
    delete solverMpK;
    solverMpK = NULL;
    delete solverS;
    solverS = NULL;
  }
221 222


223 224 225 226
  PetscSolver* PetscSolverPfcDiag::createSubSolver(int component,
					      string kspPrefix)
  {
    FUNCNAME("PetscSolverPfcDiag::createSubSolver()");
227

228 229
    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);
230

231 232 233 234
    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
    subSolver->setKspPrefix(kspPrefix.c_str());
    subSolver->setMeshDistributor(meshDistributor, 0);
    subSolver->init(fe, fe);
235

236 237 238
    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();
239

240 241 242
    return subSolver;
  }
} }