PetscSolverPfc.cc 7.84 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
82
83
#if (PETSC_VERSION_MINOR >= 5)
    KSPSetOperators(ksp, getMatInterior(), getMatInterior());
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
84
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
85
#endif
Praetorius, Simon's avatar
Praetorius, Simon committed
86
    if (getInfo() >= 10)
87
      KSPMonitorSet(ksp, PETSC_MONITOR_CAST(KSPMonitorDefault), PETSC_NULL, PETSC_NULL);
Praetorius, Simon's avatar
Praetorius, Simon committed
88
    else if (getInfo() >= 20)
Praetorius, Simon's avatar
Praetorius, Simon committed
89
      KSPMonitorSet(ksp, PETSC_MONITOR_CAST(KSPMonitorTrueResidualNorm), PETSC_NULL, PETSC_NULL);
Praetorius, Simon's avatar
Praetorius, Simon committed
90
91
    petsc_helper::setSolver(ksp, "pfc_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
    KSPSetFromOptions(ksp);
92
93


Praetorius, Simon's avatar
Praetorius, Simon committed
94
95
96
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }
97

Praetorius, Simon's avatar
Praetorius, Simon committed
98
99
100
101
102
103

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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
105
106
107
    PCSetType(getPc(), PCSHELL);
    PCShellSetApply(getPc(), pcPfcShell);
    PCShellSetContext(getPc(), &data);
108

109
    double delta = std::sqrt((*tau) * M0);
110

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

Praetorius, Simon's avatar
Praetorius, Simon committed
113
114
    // create mass-matrix
    DOFMatrix matrixM(feSpace, feSpace);
115
116
117
    Operator opM(feSpace, feSpace);
    addZOT(opM, 1.0);
    matrixM.assembleOperator(opM);
118

Praetorius, Simon's avatar
Praetorius, Simon committed
119
120
121
122
    solverM = createSubSolver(0, "M_");
    solverM->fillPetscMatrix(&matrixM);
    data.matM = solverM->getMatInterior();
    data.kspM = solverM->getSolver();
123

Praetorius, Simon's avatar
Praetorius, Simon committed
124
    // create MpK-matrix
Praetorius, Simon's avatar
Praetorius, Simon committed
125
    DOFMatrix matrixMpK(feSpace, feSpace);
126
127
128
129
    Operator opMpK(feSpace, feSpace);
    addZOT(opMpK, 1.0);
    addSOT(opMpK, delta);
    matrixMpK.assembleOperator(opMpK);
130

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


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

Praetorius, Simon's avatar
Praetorius, Simon committed
144
145
146
147
    solverMpK2 = createSubSolver(0, "MpK2_");
    solverMpK2->fillPetscMatrix(&matrixMpK2);
    matMpK2 = solverMpK2->getMatInterior();
    data.kspMpK2 = solverMpK2->getSolver();
148

Praetorius, Simon's avatar
Praetorius, Simon committed
149
150
    // create laplace-matrix
    DOFMatrix matrixK(feSpace, feSpace);
151
152
153
    Operator opK(feSpace, feSpace);
    addSOT(opK, 1.0);
    matrixK.assembleOperator(opK);
154

Praetorius, Simon's avatar
Praetorius, Simon committed
155
156
157
    solverK = createSubSolver(0, "K_");
    solverK->fillPetscMatrix(&matrixK);
    data.matK = solverK->getMatInterior();
158
159

    // === Setup preconditioner data ===
Praetorius, Simon's avatar
Praetorius, Simon committed
160
    data.delta = delta;
Praetorius, Simon's avatar
Praetorius, Simon committed
161
    data.tau = (*tau);
162

163
164
165
    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;
166
    bool directM = false, directMpL = false, directMpL2 = false;
167
168
169
170
171
172
173
174
175
    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);
176

177
178
179
    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);
180

181
182
    bool useAMG = false;
    Parameters::get("precon_pfc_MpL2->use AMG", useAMG);
Praetorius, Simon's avatar
Praetorius, Simon committed
183

184
185
186
187
    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);
188

189
190
191
192
    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);
193

194
195
196
197
198
199
    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
200
201
202
203
204
205
  }


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

Praetorius, Simon's avatar
Praetorius, Simon committed
207
208
//     MatDestroy(&matMpK);
//     MatDestroy(&matMpK2);
209

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

Praetorius, Simon's avatar
Praetorius, Simon committed
215
216
217
    solverM->destroyVectorData();
    solverK->destroyVectorData();
    solverMpK->destroyVectorData();
Praetorius, Simon's avatar
Praetorius, Simon committed
218
    solverMpK2->destroyVectorData();
219

Praetorius, Simon's avatar
Praetorius, Simon committed
220
221
222
223
224
225
    delete solverM;
    solverM = NULL;
    delete solverK;
    solverK = NULL;
    delete solverMpK;
    solverMpK = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
226
227
    delete solverMpK2;
    solverMpK2 = NULL;
Praetorius, Simon's avatar
Praetorius, Simon committed
228
  }
229
230


Praetorius, Simon's avatar
Praetorius, Simon committed
231
232
233
234
  PetscSolver* PetscSolverPfc::createSubSolver(int component,
					      string kspPrefix)
  {
    FUNCNAME("PetscSolverPfc::createSubSolver()");
235

Praetorius, Simon's avatar
Praetorius, Simon committed
236
237
    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);
238

Praetorius, Simon's avatar
Praetorius, Simon committed
239
240
241
242
    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
    subSolver->setKspPrefix(kspPrefix.c_str());
    subSolver->setMeshDistributor(meshDistributor, 0);
    subSolver->init(fe, fe);
243

Praetorius, Simon's avatar
Praetorius, Simon committed
244
245
246
    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();
247

Praetorius, Simon's avatar
Praetorius, Simon committed
248
249
250
    return subSolver;
  }
} }