PetscSolverCahnHilliard.cc 7.91 KB
Newer Older
1
2
3
4
5
6
7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9
10
11
12
13
14
15
16
17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, 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.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19
 ******************************************************************************/
20
21
22
23
24
25
26



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

27
namespace AMDiS { namespace Parallel {
28

29
  using namespace std;
30

31
32
  /// solve Cahn-Hilliard Preconditioner
  PetscErrorCode pcChShell(PC pc, Vec b, Vec x) // solve Px=b
33
  {
34
35
36
    void *ctx;
    PCShellGetContext(pc, &ctx);
    CahnHilliardData* data = static_cast<CahnHilliardData*>(ctx);
37

38
39
40
    Vec b1, b2, x1, x2;
    VecNestGetSubVec(b, 0, &b1);
    VecNestGetSubVec(b, 1, &b2);
41

42
43
    VecNestGetSubVec(x, 0, &x1);
    VecNestGetSubVec(x, 1, &x2);
44

45
46
47
    Vec y1, y2;
    VecDuplicate(b1, &y1);
    VecDuplicate(b2, &y2);
48

49
50
51
    //     MatGetDiagonal(data->matM, y2);
    //     VecReciprocal(y2);
    //     VecPointwiseMult(y1, y2, b1);
52
    KSPSolve(data->kspMass, b1, y1); 			// M*y1 = b1
53
    MatMultAdd(data->matMinusDeltaK, y1, b2, x1); 	// -> x1 := b2-delta*K*y1
54

55
56
    KSPSolve(data->kspLaplace, x1, y2); 			// (M+eps*sqrt(delta))*y2 = x1
    MatMult(data->matM, y2, x1); 			// x1 := M*y2
57

58
59
60
61
    KSPSolve(data->kspLaplace, x1, x2);			// (M+eps*sqrt(delta))*x2 = x1
    double factor = (*data->eps)/sqrt(*data->delta);
    VecCopy(x2, x1); 					// x1 := x2
    VecAXPBYPCZ(x1, 1.0, factor, -factor, y1, y2);	// x1 = 1*y1 + factor*y2 - factor*x1
62

63
64
    VecDestroy(&y1);
    VecDestroy(&y2);
65

66
    PetscFunctionReturn(0);
67
  }
68
69


70
  PetscSolverCahnHilliard::PetscSolverCahnHilliard(string name, double *epsPtr, double *deltaPtr)
71
  : PetscSolverGlobalBlockMatrix(name),
72
73
74
    massMatrixSolver(NULL),
    laplaceMatrixSolver(NULL),
    deltaKMatrixSolver(NULL),
75
    useOldInitialGuess(false),
76
    phase(NULL),
77
    eps(epsPtr),
78
    delta(deltaPtr),
79
    tau(NULL)
80
  {
81
82
    Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess);
  }
83

84
85
86
87
  void PetscSolverCahnHilliard::initSolver(KSP &ksp)
  {
    // Create FGMRES based outer solver
    KSPCreate(meshDistributor->getMpiComm(0), &ksp);
88
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
89
    if (getInfo() >= 10)
90
      petsc::ksp_monitor_set(ksp, KSPMonitorDefault);
91
    else if (getInfo() >= 20)
92
      petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
93
    petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
94
    KSPSetFromOptions(ksp);
95

96
97
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
98
  }
99
100


101
102
103
104
  void PetscSolverCahnHilliard::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverCahnHilliard::initPreconditioner()");
    MSG("PetscSolverCahnHilliard::initPreconditioner()\n");
105

106
107
108
109
    if (tau) {
      delta = new double;
      *delta = (*gamma) * (*tau);
    }
110

111
    TEST_EXIT(eps && delta)("eps and/or delta pointers not set!\n");
112
113


114
    //     KSPSetUp(kspInterior);
115

116
117
118
    PCSetType(pc, PCSHELL);
    PCShellSetApply(pc, pcChShell);
    PCShellSetContext(pc, &matShellContext);
119

120
    const FiniteElemSpace *feSpace = componentSpaces[0];
121
122

    // === matrix M ===
123
124
    DOFMatrix laplaceMatrix(feSpace, feSpace);
    Operator laplaceOp(feSpace, feSpace);
125

126
127
    DOFMatrix massMatrix(feSpace, feSpace);
    Operator massOp(feSpace, feSpace);
128

129
130
    DOFMatrix deltaKMatrix(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
131
132

    if (phase) {
133
      VecAtQP_ZOT zot(phase, NULL);
134
135
      massOp.addTerm(&zot);
      laplaceOp.addTerm(&zot); // M
136
      VecAtQP_SOT sot2(phase, NULL, (*eps)*sqrt(*delta));
137
      laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
138
      VecAtQP_SOT sot(phase, NULL, -(*delta));
139
140
141
142
      laplaceOp2.addTerm(&sot); // -delta*K
      massMatrix.assembleOperator(massOp);
      massMatrixSolver = createSubSolver(0, "mass_");
      massMatrixSolver->fillPetscMatrix(&massMatrix);
143

144
145
146
147
      // === matrix (M + eps*sqrt(delta)*K) ===
      laplaceMatrix.assembleOperator(laplaceOp);
      laplaceMatrixSolver = createSubSolver(0, "laplace_");
      laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
148

149
150
151
152
      // === matrix (-delta*K) ===
      deltaKMatrix.assembleOperator(laplaceOp2);
      deltaKMatrixSolver = createSubSolver(0, "laplace2_");
      deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
153
154
155


    } else {
156
157
158
159
160
      Simple_ZOT zot;
      massOp.addTerm(&zot);
      laplaceOp.addTerm(&zot); // M
      Simple_SOT sot2((*eps)*sqrt(*delta));
      laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
161
      Simple_SOT sot(-(*delta));
162
      laplaceOp2.addTerm(&sot); // -delta*K
163

164
165
166
      massMatrix.assembleOperator(massOp);
      massMatrixSolver = createSubSolver(0, "mass_");
      massMatrixSolver->fillPetscMatrix(&massMatrix);
167

168
169
170
171
      // === matrix (M + eps*sqrt(delta)*K) ===
      laplaceMatrix.assembleOperator(laplaceOp);
      laplaceMatrixSolver = createSubSolver(0, "laplace_");
      laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
172

173
174
175
176
177
      // === matrix (-delta*K) ===
      deltaKMatrix.assembleOperator(laplaceOp2);
      deltaKMatrixSolver = createSubSolver(0, "laplace2_");
      deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
    }
178
179
180
181




182
183
184
    // === Setup solver ===
    matShellContext.kspMass = massMatrixSolver->getSolver();
    matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
185
186
    matShellContext.matM = massMatrixSolver->getMatInterior();
    matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();
187
188
    matShellContext.eps = eps;
    matShellContext.delta = delta;
189

190
    matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0));
191

192
    petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCBJACOBI, 0.0, 1e-14, 2);
193
194
195
196
197
198
    //     petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
    //     {
      //             PC pc;
    //       KSPGetPC(matShellContext.kspMass, &pc);
    //       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
    //     }
199

200
    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 2);
201
202
203
204
205
206
207
    //     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
    //     {
      //             PC pc;
    //       KSPGetPC(matShellContext.kspLaplace, &pc);
    //       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
    //     }
    PCSetFromOptions(pc);
208
  }
209
210


211
212
213
214
  PetscSolver* PetscSolverCahnHilliard::createSubSolver(int component,
							string kspPrefix)
  {
    FUNCNAME("PetscSolverCahnHilliard::createSubSolver()");
215

216
217
    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);
218

219
220
221
222
    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
    subSolver->setKspPrefix(kspPrefix.c_str());
    subSolver->setMeshDistributor(meshDistributor, 0);
    subSolver->init(fe, fe);
223

224
225
226
    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();
227

228
229
    return subSolver;
  }
230

231
232
233
  void PetscSolverCahnHilliard::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");
234

235
236
237
    massMatrixSolver->destroyMatrixData();
    laplaceMatrixSolver->destroyMatrixData();
    deltaKMatrixSolver->destroyMatrixData();
238

239
240
241
    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    deltaKMatrixSolver->destroyVectorData();
242
243


244
    delete massMatrixSolver;
245
    massMatrixSolver = NULL;
246

247
    delete laplaceMatrixSolver;
248
    laplaceMatrixSolver = NULL;
249

250
    delete deltaKMatrixSolver;
251
    deltaKMatrixSolver = NULL;
252

253
254
255
    if (tau) {
      delete delta;
    }
256
  }
257
258


259
} }