PetscSolverPfc_diag.cc 6.94 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
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
80
    if (getInfo() >= 10)
81
      petsc::ksp_monitor_set(ksp, KSPMonitorDefault);
82
    else if (getInfo() >= 20)
83
      petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
84
85
    petsc_helper::setSolver(ksp, "pfc_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
    KSPSetFromOptions(ksp);
86
87


88
89
90
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }
91

92
93
94
95
96
97

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

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

99
100
101
    PCSetType(getPc(), PCSHELL);
    PCShellSetApply(getPc(), pcPfcDiagShell);
    PCShellSetContext(getPc(), &data);
102

103
    double delta = sqrt(*tau);
104

105
    const FiniteElemSpace *feSpace = componentSpaces[0];
106

107
108
109
110
111
112
    // create mass-matrix
    DOFMatrix matrixM(feSpace, feSpace);
    Operator massOp(feSpace, feSpace);
    Simple_ZOT zot;
    massOp.addTerm(&zot);
    matrixM.assembleOperator(massOp);
113

114
115
116
117
    solverM = createSubSolver(0, "M_");
    solverM->fillPetscMatrix(&matrixM);
    data.matM = solverM->getMatInterior();
    data.kspM = solverM->getSolver();
118

119
120
121
122
123
124
125
    // create MpK-matrix
    DOFMatrix matrixMpK(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
    Simple_SOT sot2(delta);
    laplaceOp2.addTerm(&zot);
    laplaceOp2.addTerm(&sot2);
    matrixMpK.assembleOperator(laplaceOp2);
126

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


133
    // create MpK2-matrix
134

135
136
137
138
139
140
    // create laplace-matrix
    DOFMatrix matrixK(feSpace, feSpace);
    Operator laplaceOp(feSpace, feSpace);
    Simple_SOT sot;
    laplaceOp.addTerm(&sot);
    matrixK.assembleOperator(laplaceOp);
141

142
143
144
    solverK = createSubSolver(0, "K_");
    solverK->fillPetscMatrix(&matrixK);
    data.matK = solverK->getMatInterior();
145

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

152
153
154
155
    Vec x; // need to initialize vector
    solverS->createVec(*solverS->getDofMapping(), x);
    Mat DK;
    MatDuplicate(data.matK, MAT_COPY_VALUES, &DK);
156

157
158
159
160
161
162
163
164
    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);
165

166
    petsc::ksp_set_operators(data.kspS, matS, matS);
167
168

    // === Setup preconditioner data ===
169
170
    data.delta = delta;
    data.tau = (*tau);
171

172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
    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()");
190

191
192
//     MatDestroy(&matMpK);
//     MatDestroy(&matMpK2);
193

194
195
196
197
    solverM->destroyMatrixData();
    solverK->destroyMatrixData();
    solverMpK->destroyMatrixData();
    solverS->destroyMatrixData();
198

199
200
201
202
    solverM->destroyVectorData();
    solverK->destroyVectorData();
    solverMpK->destroyVectorData();
    solverS->destroyVectorData();
203

204
205
206
207
208
209
210
211
212
    delete solverM;
    solverM = NULL;
    delete solverK;
    solverK = NULL;
    delete solverMpK;
    solverMpK = NULL;
    delete solverS;
    solverS = NULL;
  }
213
214


215
216
217
218
  PetscSolver* PetscSolverPfcDiag::createSubSolver(int component,
					      string kspPrefix)
  {
    FUNCNAME("PetscSolverPfcDiag::createSubSolver()");
219

220
221
    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);
222

223
224
225
226
    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
    subSolver->setKspPrefix(kspPrefix.c_str());
    subSolver->setMeshDistributor(meshDistributor, 0);
    subSolver->init(fe, fe);
227

228
229
230
    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();
231

232
233
234
    return subSolver;
  }
} }