PetscSolverCahnHilliard.cc 8.57 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
/******************************************************************************
 *
 * 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 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.
 * 
 ******************************************************************************/
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
33
34
  /// solve Cahn-Hilliard Preconditioner
  PetscErrorCode pcChShell(PC pc, Vec b, Vec x) // solve Px=b
  {FUNCNAME("PCApply()");
  
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
  void *ctx;
  PCShellGetContext(pc, &ctx);
  CahnHilliardData* data = static_cast<CahnHilliardData*>(ctx);
  
  Vec b1, b2, x1, x2;
  VecNestGetSubVec(b, 0, &b1);
  VecNestGetSubVec(b, 1, &b2);
  
  VecNestGetSubVec(x, 0, &x1);
  VecNestGetSubVec(x, 1, &x2);
  
  Vec y1, y2;
  VecDuplicate(b1, &y1);
  VecDuplicate(b2, &y2);
  
  //     MatGetDiagonal(data->matM, y2);
  //     VecReciprocal(y2);
  //     VecPointwiseMult(y1, y2, b1);
  KSPSolve(data->kspMass, b1, y1); 			// M*y1 = b1    
  MatMultAdd(data->matMinusDeltaK, y1, b2, x1); 	// -> x1 := b2-delta*K*y1
  
  KSPSolve(data->kspLaplace, x1, y2); 			// (M+eps*sqrt(delta))*y2 = x1
  MatMult(data->matM, y2, x1); 			// x1 := M*y2
  
  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
  
  VecDestroy(&y1);
  VecDestroy(&y2);
  
67
  PetscFunctionReturn(0);
68
69
  }
  
70
  
71
  PetscSolverCahnHilliard::PetscSolverCahnHilliard(string name, double *epsPtr, double *deltaPtr)
72
  : PetscSolverGlobalBlockMatrix(name),
73
74
75
76
77
78
79
    massMatrixSolver(NULL),
    laplaceMatrixSolver(NULL),
    deltaKMatrixSolver(NULL),
    useOldInitialGuess(false),
    eps(epsPtr),
    delta(deltaPtr) 
  { 	
80
81
82
83
    Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess);
  }
  
  
84
85
86
87
  void PetscSolverCahnHilliard::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverCahnHilliard::solvePetscMatrix()");
88
89
90
91
    
  /*  if (useOldInitialGuess) {
      if (getVecSolInterior())
      {VecSet(getVecSolInterior(), 0.0);
92
93
94
95
      
      for (int i = 0; i < solution->getSize(); i++)
      {
	Vec tmp;
96
	VecNestGetSubVec(getVecSolInterior(), i, &tmp);
97
98
99
100
101
102
	setDofVector(tmp, solution->getDOFVector(i));
      }
      
      vecSolAssembly();
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
    }
103
104
    KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
    }*/
105
106
    PetscSolverGlobalBlockMatrix::solvePetscMatrix(vec, adaptInfo);
  }
107
  
108
109
110
111
112
  void PetscSolverCahnHilliard::initSolver(KSP &ksp)
  {
    // Create FGMRES based outer solver
    KSPCreate(meshDistributor->getMpiComm(0), &ksp);
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
113
114
115
116
117
    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, "ch_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
118
    KSPSetFromOptions(ksp);
119
  }
120
121
  
  
122
123
124
125
  void PetscSolverCahnHilliard::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverCahnHilliard::initPreconditioner()");
    MSG("PetscSolverCahnHilliard::initPreconditioner()\n");
126
    
127
128
    TEST_EXIT(eps && delta)("eps and/or delta pointers not set!\n");
    
129
130
    //     KSPSetUp(kspInterior);
    
131
132
133
    PCSetType(pc, PCSHELL);
    PCShellSetApply(pc, pcChShell);
    PCShellSetContext(pc, &matShellContext);
134
    
135
136
137
138
139
140
    const FiniteElemSpace *feSpace = componentSpaces[0];
    
    // === matrix M ===	
    DOFMatrix laplaceMatrix(feSpace, feSpace);
    Operator laplaceOp(feSpace, feSpace);
    
141
142
143
    DOFMatrix massMatrix(feSpace, feSpace);
    Operator massOp(feSpace, feSpace);
    
144
145
    DOFMatrix deltaKMatrix(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
146
    
147
    if (phase) {  
148
149
150
      VecAtQP_ZOT zot(phase, NULL);
      massOp.addTerm(&zot);
      laplaceOp.addTerm(&zot); // M
151
      VecAtQP_SOT sot2(phase, NULL, (*eps)*sqrt(*delta));
152
      laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
153
      VecAtQP_SOT sot(phase, NULL, -(*delta));
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
      laplaceOp2.addTerm(&sot); // -delta*K
      massMatrix.assembleOperator(massOp);
      massMatrixSolver = createSubSolver(0, "mass_");
      massMatrixSolver->fillPetscMatrix(&massMatrix);
      
      // === matrix (M + eps*sqrt(delta)*K) ===
      laplaceMatrix.assembleOperator(laplaceOp);
      laplaceMatrixSolver = createSubSolver(0, "laplace_");
      laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
      
      // === matrix (-delta*K) ===
      deltaKMatrix.assembleOperator(laplaceOp2);
      deltaKMatrixSolver = createSubSolver(0, "laplace2_");
      deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
      
      
170
    } else { 
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
      Simple_ZOT zot;
      massOp.addTerm(&zot);
      laplaceOp.addTerm(&zot); // M
      Simple_SOT sot2((*eps)*sqrt(*delta));
      laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
      Simple_SOT sot(-(*delta));    
      laplaceOp2.addTerm(&sot); // -delta*K
      
      massMatrix.assembleOperator(massOp);
      massMatrixSolver = createSubSolver(0, "mass_");
      massMatrixSolver->fillPetscMatrix(&massMatrix);
      
      // === matrix (M + eps*sqrt(delta)*K) ===
      laplaceMatrix.assembleOperator(laplaceOp);
      laplaceMatrixSolver = createSubSolver(0, "laplace_");
      laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
      
      // === matrix (-delta*K) ===
      deltaKMatrix.assembleOperator(laplaceOp2);
      deltaKMatrixSolver = createSubSolver(0, "laplace2_");
      deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
    }
    
    
    
196
197
198
199
200
201
202
203
204
205
    
    // === Setup solver ===
    matShellContext.kspMass = massMatrixSolver->getSolver();
    matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
    matShellContext.matM = massMatrixSolver->getMatInterior();    
    matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();   
    matShellContext.eps = eps;
    matShellContext.delta = delta;
    
    matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0));
206
    
207
    petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
208
209
210
211
212
213
214
    //     petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
    //     {
      //             PC pc;
    //       KSPGetPC(matShellContext.kspMass, &pc);
    //       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
    //     }
    
215
    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 2);
216
217
218
219
220
221
222
    //     petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1);
    //     {
      //             PC pc;
    //       KSPGetPC(matShellContext.kspLaplace, &pc);
    //       PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS);
    //     }
    PCSetFromOptions(pc);
223
  }
224
  
225
226
227
228
229
  
  PetscSolver* PetscSolverCahnHilliard::createSubSolver(int component,
							string kspPrefix)
  {
    FUNCNAME("PetscSolverCahnHilliard::createSubSolver()");
230
    
231
232
    vector<const FiniteElemSpace*> fe;
    fe.push_back(componentSpaces[component]);
233
    
234
235
236
237
    PetscSolver* subSolver = new PetscSolverGlobalMatrix("");
    subSolver->setKspPrefix(kspPrefix.c_str());
    subSolver->setMeshDistributor(meshDistributor, 0);
    subSolver->init(fe, fe);
238
    
239
240
241
    ParallelDofMapping &subDofMap = subSolver->getDofMap();
    subDofMap[0] = dofMap[component];
    subDofMap.update();
242
    
243
244
245
246
247
248
249
250
251
252
    return subSolver;
  }
  
  void PetscSolverCahnHilliard::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()");
    
    massMatrixSolver->destroyMatrixData();
    laplaceMatrixSolver->destroyMatrixData();
    deltaKMatrixSolver->destroyMatrixData();
253
    
254
255
256
    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    deltaKMatrixSolver->destroyVectorData();
257
    
258
259
260
    
    delete massMatrixSolver;
    massMatrixSolver = NULL;
261
    
262
263
    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;
264
    
265
266
267
    delete deltaKMatrixSolver;
    deltaKMatrixSolver = NULL;
  }
268
269
  
  
270
} }