PetscSolverPfc.cc 7.65 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
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
 *
Praetorius, Simon's avatar
Praetorius, Simon committed
16 17 18 19 20
 ******************************************************************************/

#include "PetscSolverPfc.h"
#include "parallel/PetscHelper.h"
#include "parallel/PetscSolverGlobalMatrix.h"
21
#include "Expressions.h"
Praetorius, Simon's avatar
Praetorius, Simon committed
22 23 24 25

namespace AMDiS { namespace Parallel {

  using namespace std;
26

Praetorius, Simon's avatar
Praetorius, Simon committed
27 28 29
  /// solve Pfc Preconditioner
  PetscErrorCode pcPfcShell(PC pc, Vec b, Vec x) // solve Px=b
  { FUNCNAME("pcPfcShell()");
30

Praetorius, Simon's avatar
Praetorius, Simon committed
31 32 33
    void *ctx;
    PCShellGetContext(pc, &ctx);
    PfcData* data = static_cast<PfcData*>(ctx);
34

Praetorius, Simon's avatar
Praetorius, Simon committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50
    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);

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

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

57 58 59 60
    MatMult(data->matM, y2, tmp);		// tmp := M*y2
    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
61

62
    VecCopy(x2, x1); 				// x1 := x2
Praetorius, Simon's avatar
Praetorius, Simon committed
63 64
    VecAXPBYPCZ(x1, 1.0, 1.0/(data->delta), -1.0/(data->delta), y1, y2);	// x1 = 1*y1 + factor*y2 - factor*x1

65 66
    MatMult(data->matK, x2, tmp); 		// tmp := K*x2
    VecAYPX(tmp, -1.0, b3);			// tmp := b3 - tmp
67 68

    KSPSolve(data->kspM, tmp, x3); 		// M*x3 = tmp
Praetorius, Simon's avatar
Praetorius, Simon committed
69 70 71 72

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

Praetorius, Simon's avatar
Praetorius, Simon committed
74
    return 0;
Praetorius, Simon's avatar
Praetorius, Simon committed
75
  }
76

Praetorius, Simon's avatar
Praetorius, Simon committed
77 78 79 80
  void PetscSolverPfc::initSolver(KSP &ksp)
  {
    // Create FGMRES based outer solver
    KSPCreate(meshDistributor->getMpiComm(0), &ksp);
81
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
Praetorius, Simon's avatar
Praetorius, Simon committed
82
    if (getInfo() >= 10)
83
      petsc::ksp_monitor_set(ksp, KSPMonitorDefault);
Praetorius, Simon's avatar
Praetorius, Simon committed
84
    else if (getInfo() >= 20)
85
      petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
Praetorius, Simon's avatar
Praetorius, Simon committed
86 87
    petsc_helper::setSolver(ksp, "pfc_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
    KSPSetFromOptions(ksp);
88 89


Praetorius, Simon's avatar
Praetorius, Simon committed
90 91 92
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }
93

Praetorius, Simon's avatar
Praetorius, Simon committed
94 95 96 97 98 99

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
101 102 103
    PCSetType(getPc(), PCSHELL);
    PCShellSetApply(getPc(), pcPfcShell);
    PCShellSetContext(getPc(), &data);
104

105
    double delta = std::sqrt((*tau) * M0);
106

Praetorius, Simon's avatar
Praetorius, Simon committed
107
    const FiniteElemSpace *feSpace = componentSpaces[0];
108

Praetorius, Simon's avatar
Praetorius, Simon committed
109 110
    // create mass-matrix
    DOFMatrix matrixM(feSpace, feSpace);
111 112 113
    Operator opM(feSpace, feSpace);
    addZOT(opM, 1.0);
    matrixM.assembleOperator(opM);
114

Praetorius, Simon's avatar
Praetorius, Simon committed
115 116 117 118
    solverM = createSubSolver(0, "M_");
    solverM->fillPetscMatrix(&matrixM);
    data.matM = solverM->getMatInterior();
    data.kspM = solverM->getSolver();
119

Praetorius, Simon's avatar
Praetorius, Simon committed
120
    // create MpK-matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
121
    DOFMatrix matrixMpK(feSpace, feSpace);
122 123 124 125
    Operator opMpK(feSpace, feSpace);
    addZOT(opMpK, 1.0);
    addSOT(opMpK, delta);
    matrixMpK.assembleOperator(opMpK);
126

Praetorius, Simon's avatar
Praetorius, Simon committed
127 128 129 130
    solverMpK = createSubSolver(0, "MpK_");
    solverMpK->fillPetscMatrix(&matrixMpK);
    matMpK = solverMpK->getMatInterior();
    data.kspMpK = solverMpK->getSolver();
131 132


Praetorius, Simon's avatar
Praetorius, Simon committed
133 134
    // create MpK2-matrix
    DOFMatrix matrixMpK2(feSpace, feSpace);
135 136 137 138
    Operator opMpK2(feSpace, feSpace);
    addZOT(opMpK2, 1.0);
    addSOT(opMpK2, std::sqrt(delta));
    matrixMpK2.assembleOperator(opMpK2);
139

Praetorius, Simon's avatar
Praetorius, Simon committed
140 141 142 143
    solverMpK2 = createSubSolver(0, "MpK2_");
    solverMpK2->fillPetscMatrix(&matrixMpK2);
    matMpK2 = solverMpK2->getMatInterior();
    data.kspMpK2 = solverMpK2->getSolver();
144

Praetorius, Simon's avatar
Praetorius, Simon committed
145 146
    // create laplace-matrix
    DOFMatrix matrixK(feSpace, feSpace);
147 148 149
    Operator opK(feSpace, feSpace);
    addSOT(opK, 1.0);
    matrixK.assembleOperator(opK);
150

Praetorius, Simon's avatar
Praetorius, Simon committed
151 152 153
    solverK = createSubSolver(0, "K_");
    solverK->fillPetscMatrix(&matrixK);
    data.matK = solverK->getMatInterior();
154 155

    // === Setup preconditioner data ===
Praetorius, Simon's avatar
Praetorius, Simon committed
156
    data.delta = delta;
Praetorius, Simon's avatar
Praetorius, Simon committed
157
    data.tau = (*tau);
158

159 160 161
    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;
162
    bool directM = false, directMpL = false, directMpL2 = false;
163 164 165 166 167 168 169 170 171
    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);
172

173 174 175
    Parameters::get("precon_pfc_M->use direct solver", directM);
    Parameters::get("precon_pfc_MpL->use direct solver", directMpL);
    Parameters::get("precon_pfc_MpL2->use direct solver", directMpL2);
176

177 178
    bool useAMG = false;
    Parameters::get("precon_pfc_MpL2->use AMG", useAMG);
Praetorius, Simon's avatar
Praetorius, Simon committed
179

180 181 182 183
    if (directM)
      petsc_helper::setSolverWithLu(data.kspM, "M_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
    else
      petsc_helper::setSolver(data.kspM, "M_", KSPCG, PCBJACOBI, rtolM, tolM, nIterM);
184

185 186 187 188
    if (directMpL)
      petsc_helper::setSolverWithLu(data.kspMpK, "MpL_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
    else
      petsc_helper::setSolver(data.kspMpK, "MpL_", KSPCG, PCBJACOBI, rtolMpL, tolMpL, nIterMpL);
189

190 191 192 193 194 195
    if (directMpL2)
      petsc_helper::setSolverWithLu(data.kspMpK2, "MpL2_", KSPRICHARDSON, PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
    else if (useAMG)
      petsc_helper::setSolver(data.kspMpK2, "MpL2_", KSPRICHARDSON, PCHYPRE, rtolMpL2, tolMpL2, nIterMpL2);
    else
      petsc_helper::setSolver(data.kspMpK2, "MpL2_", KSPCG, PCBJACOBI, rtolMpL2, tolMpL2, nIterMpL2);
Praetorius, Simon's avatar
Praetorius, Simon committed
196 197 198 199 200 201
  }


  void PetscSolverPfc::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverPfc::exit()");
202

Praetorius, Simon's avatar
Praetorius, Simon committed
203 204
//     MatDestroy(&matMpK);
//     MatDestroy(&matMpK2);
205

Praetorius, Simon's avatar
Praetorius, Simon committed
206 207 208
    solverM->destroyMatrixData();
    solverK->destroyMatrixData();
    solverMpK->destroyMatrixData();
Praetorius, Simon's avatar
Praetorius, Simon committed
209
    solverMpK2->destroyMatrixData();
210

Praetorius, Simon's avatar
Praetorius, Simon committed
211 212 213
    solverM->destroyVectorData();
    solverK->destroyVectorData();
    solverMpK->destroyVectorData();
Praetorius, Simon's avatar
Praetorius, Simon committed
214
    solverMpK2->destroyVectorData();
215

Praetorius, Simon's avatar
Praetorius, Simon committed
216 217 218 219 220 221
    delete solverM;
    solverM = NULL;
    delete solverK;
    solverK = NULL;
    delete solverMpK;
    solverMpK = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
222 223
    delete solverMpK2;
    solverMpK2 = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
224
  }
225 226


Praetorius, Simon's avatar
Praetorius, Simon committed
227 228 229 230
  PetscSolver* PetscSolverPfc::createSubSolver(int component,
					      string kspPrefix)
  {
    FUNCNAME("PetscSolverPfc::createSubSolver()");
231

Praetorius, Simon's avatar
Praetorius, Simon committed
232 233
    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);
234

Praetorius, Simon's avatar
Praetorius, Simon committed
235 236 237 238
    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
    subSolver->setKspPrefix(kspPrefix.c_str());
    subSolver->setMeshDistributor(meshDistributor, 0);
    subSolver->init(fe, fe);
239

Praetorius, Simon's avatar
Praetorius, Simon committed
240 241 242
    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();
243

Praetorius, Simon's avatar
Praetorius, Simon committed
244 245 246
    return subSolver;
  }
} }