PetscSolverNSCH.cc 14.8 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
 ******************************************************************************/
Sebastian Aland's avatar
Sebastian Aland committed
20
21
22
23
24
25
26
27
28
29



#include "parallel/PetscSolverNSCH.h"
#include "parallel/PetscHelper.h"
#include "TransformDOF.h"




30
namespace AMDiS { namespace Parallel {
31

Sebastian Aland's avatar
Sebastian Aland committed
32
  using namespace std;
33
34


Sebastian Aland's avatar
Sebastian Aland committed
35
36
  PetscErrorCode pcShell(PC pc, Vec b, Vec x) // solve Px=b
  {FUNCNAME("PCApply()");
37

Sebastian Aland's avatar
Sebastian Aland committed
38
39
  void *ctx;
  PCShellGetContext(pc, &ctx);
40
  NSCHData* data = static_cast<NSCHData*>(ctx);
41

Sebastian Aland's avatar
Sebastian Aland committed
42
  /// extract vectors
43
44
45
  Vec b1, b2, b34, b5, x1, x2, x34, x5;
  data->globalMatrixSolver->extractVectorComponent(b, data->dim+1, &b1);
  data->globalMatrixSolver->extractVectorComponent(b, data->dim+2, &b2);
46
  data->globalMatrixSolver->extractVectorComponent(b, 0, &b34, data->dim);
47
  data->globalMatrixSolver->extractVectorComponent(b, data->dim, &b5);
48

49
50
  data->globalMatrixSolver->extractVectorComponent(x, data->dim+1, &x1);
  data->globalMatrixSolver->extractVectorComponent(x, data->dim+2, &x2);
51
  data->globalMatrixSolver->extractVectorComponent(x, 0, &x34, data->dim);
52
  data->globalMatrixSolver->extractVectorComponent(x, data->dim, &x5);
53
54
55


  /// solve Cahn-Hilliard Preconditioner
Sebastian Aland's avatar
Sebastian Aland committed
56
57
58
    Vec y1, y2;
     VecDuplicate(b1, &y1);
     VecDuplicate(b2, &y2);
59
60

     KSPSolve(data->kspMassCH, b1, y1); 			// M*y1 = b1
Sebastian Aland's avatar
Sebastian Aland committed
61
     MatMultAdd(data->matMinusDeltaK, y1, b2, x1); 	// -> x1 := b2-delta*K*y1
62

Sebastian Aland's avatar
Sebastian Aland committed
63
64
     KSPSolve(data->kspLaplaceCH, x1, y2); 			// (M+eps*sqrt(delta))*y2 = x1
     MatMult(data->matMassCH, y2, x1); 			// x1 := M*y2
65

Sebastian Aland's avatar
Sebastian Aland committed
66
67
68
69
     KSPSolve(data->kspLaplaceCH, 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
70

Sebastian Aland's avatar
Sebastian Aland committed
71
     VecDestroy(&y1);
72
73
74
75
     VecDestroy(&y2);
    /**/

  /// solve Navier-Stokes Preconditioner
Sebastian Aland's avatar
Sebastian Aland committed
76
77
78
79
80
  Vec tmp34, tmp5, tmp5_2, tmp34_2;
  VecDuplicate(b34, &tmp34);
  VecDuplicate(b34, &tmp34_2);
  VecDuplicate(b5, &tmp5);
  VecDuplicate(b5, &tmp5_2);
81

Sebastian Aland's avatar
Sebastian Aland committed
82
83
  KSPSolve(data->kspVelocity, b34, tmp34);
  VecScale(tmp34, -1.0);
84
85
  MatMultAdd(data->matDiv, tmp34, b5, tmp5);

Sebastian Aland's avatar
Sebastian Aland committed
86
87
88
89
90
91
92
93
  /// approximierte Schur Komplement
   KSPSolve(data->kspLaplace, tmp5, x5);
   { //project out constant Null-space
      int vecSize;
      VecGetSize(x5, &vecSize);
      PetscScalar vecSum;
      VecSum(x5, &vecSum);
      vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
94
      VecShift(x5, vecSum);
Sebastian Aland's avatar
Sebastian Aland committed
95
96
97
98
99
      //VecView(y, PETSC_VIEWER_STDOUT_WORLD);
    }
    MatMult(data->matConDif, x5, tmp5);
    KSPSolve(data->kspMass, tmp5, x5);
    VecScale(x5,-1.0);
100
101

  MatMultAdd(data->matGrad, x5, b34, tmp34);
Sebastian Aland's avatar
Sebastian Aland committed
102
103
  KSPSolve(data->kspVelocity, tmp34, x34);
  VecScale(x5,-1.0);
104
/**/
Sebastian Aland's avatar
Sebastian Aland committed
105
106
107
108
109
110
111
112
113
114
115
116
117

  VecDestroy(&tmp34);
  VecDestroy(&tmp5);
  VecDestroy(&tmp5_2);
  VecDestroy(&tmp34_2);
  VecDestroy(&b1);
  VecDestroy(&b2);
  VecDestroy(&b34);
  VecDestroy(&b5);
  VecDestroy(&x1);
  VecDestroy(&x2);
  VecDestroy(&x34);
  VecDestroy(&x5);
118

119
    PetscFunctionReturn(0);
Sebastian Aland's avatar
Sebastian Aland committed
120
  }
121
122
123



Sebastian Aland's avatar
Sebastian Aland committed
124
  PetscSolverNSCH::PetscSolverNSCH(string name)
125
  : PetscSolverGlobalMatrix(name),
126
  pressureNullSpace(true),
127
  useOldInitialGuess(false),
128
129
  velocitySolutionMode(0),
  massSolutionMode(0),
Sebastian Aland's avatar
Sebastian Aland committed
130
  laplaceSolutionMode(0),
131
  regularizeLaplace(0),
132
  massMatrixSolverCH(NULL),
133
  laplaceMatrixSolverCH(NULL),
134
135
136
137
138
139
140
141
  deltaKMatrixSolver(NULL),
  massMatrixSolver(NULL),
  laplaceMatrixSolver(NULL),
  conDifMatrixSolver(NULL),
  nu(NULL),
  invTau(NULL),
  solution(NULL),
  phase(NULL)
142
143
  {
    Parameters::get(initFileStr + "->use old initial guess",
Sebastian Aland's avatar
Sebastian Aland committed
144
		    useOldInitialGuess);
145
146

    Parameters::get(initFileStr + "->navierstokes->velocity solver",
Sebastian Aland's avatar
Sebastian Aland committed
147
		    velocitySolutionMode);
148
149

    Parameters::get(initFileStr + "->navierstokes->mass solver",
Sebastian Aland's avatar
Sebastian Aland committed
150
		    massSolutionMode);
151
152

    Parameters::get(initFileStr + "->navierstokes->laplace solver",
Sebastian Aland's avatar
Sebastian Aland committed
153
		    laplaceSolutionMode);
154
    Parameters::get(initFileStr + "->navierstokes->regularize laplace",
Sebastian Aland's avatar
Sebastian Aland committed
155
		    regularizeLaplace);
156
157


Sebastian Aland's avatar
Sebastian Aland committed
158
  }
159
160
161


  void PetscSolverNSCH::solvePetscMatrix(SystemVector &vec,
Sebastian Aland's avatar
Sebastian Aland committed
162
163
164
						  AdaptInfo *adaptInfo)
  {
    FUNCNAME("PetscSolverNSCH::solvePetscMatrix()");
165
166

       if (useOldInitialGuess) {
Sebastian Aland's avatar
Sebastian Aland committed
167
          VecSet(getVecSolInterior(), 0.0);
168

Sebastian Aland's avatar
Sebastian Aland committed
169
170
          for (int i = 0; i < solution->getSize(); i++)
    	setDofVector(getVecSolInterior(), solution->getDOFVector(i), i, true);
171

Sebastian Aland's avatar
Sebastian Aland committed
172
173
174
          vecSolAssembly();
          KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
  }
175

Sebastian Aland's avatar
Sebastian Aland committed
176
177
    PetscSolverGlobalMatrix::solvePetscMatrix(vec, adaptInfo);
  }
178
179


Sebastian Aland's avatar
Sebastian Aland committed
180
181
182
  void PetscSolverNSCH::initSolver(KSP &ksp)
  {
    FUNCNAME("PetscSolverNSCH::initSolver()");
183
184

    // Create FGMRES based outer solver
Sebastian Aland's avatar
Sebastian Aland committed
185
186
    MSG("CREATE POS 1: %p\n", &ksp);
    KSPCreate(domainComm, &ksp);
187
188
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
    petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
189
    petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCSHELL, getRelative(), getTolerance(), getMaxIterations());
Sebastian Aland's avatar
Sebastian Aland committed
190
191
    setConstantNullSpace(ksp, componentSpaces[0]->getMesh()->getDim() , true);
  }
192
193


Sebastian Aland's avatar
Sebastian Aland committed
194
195
196
  void PetscSolverNSCH::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverNSCH::initPreconditioner()");
197

Sebastian Aland's avatar
Sebastian Aland committed
198
199
    MPI::COMM_WORLD.Barrier();
    double wtime = MPI::Wtime();
200
    int dim = componentSpaces[0]->getMesh()->getDim();
Sebastian Aland's avatar
Sebastian Aland committed
201
    pressureComponent=dim;
202
    const FiniteElemSpace *cahnHilliardFeSpace = componentSpaces[dim+1];
203
//     const FiniteElemSpace *velocityFeSpace= componentSpaces[0];
Sebastian Aland's avatar
Sebastian Aland committed
204
    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
205

Sebastian Aland's avatar
Sebastian Aland committed
206
207
208
209
210
    PCSetType(pc, PCSHELL);
    PCShellSetApply(pc, pcShell);
    PCShellSetContext(pc, &matShellContext);
    matShellContext.globalMatrixSolver = (this);
    matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0));
211
    matShellContext.dim = dim;
212
213


Sebastian Aland's avatar
Sebastian Aland committed
214
    /// Init Cahn-Hilliard part
215

Sebastian Aland's avatar
Sebastian Aland committed
216
217
    DOFMatrix laplaceMatrixCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
    Operator laplaceOpCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
218

Sebastian Aland's avatar
Sebastian Aland committed
219
220
    DOFMatrix massMatrixCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
    Operator massOpCH(cahnHilliardFeSpace, cahnHilliardFeSpace);
221

Sebastian Aland's avatar
Sebastian Aland committed
222
223
    DOFMatrix deltaKMatrix(cahnHilliardFeSpace, cahnHilliardFeSpace);
    Operator laplaceOp2(cahnHilliardFeSpace, cahnHilliardFeSpace);
224
225

    {
Sebastian Aland's avatar
Sebastian Aland committed
226
227
228
229
230
      Simple_ZOT zot;
      massOpCH.addTerm(&zot);
      laplaceOpCH.addTerm(&zot); // M
      Simple_SOT sot2((*eps)*sqrt(*delta));
      laplaceOpCH.addTerm(&sot2); // eps*sqrt(delta)*K
231
      Simple_SOT sot(-(*delta));
Sebastian Aland's avatar
Sebastian Aland committed
232
      laplaceOp2.addTerm(&sot); // -delta*K
233
234
235
236

      massMatrixCH.assembleOperator(massOpCH);
      massMatrixSolverCH = createSubSolver(dim+1, "mass_");
      massMatrixSolverCH->fillPetscMatrix(&massMatrixCH);
Sebastian Aland's avatar
Sebastian Aland committed
237
238
      // === matrix (M + eps*sqrt(delta)*K) ===
      laplaceMatrixCH.assembleOperator(laplaceOpCH);
239
      laplaceMatrixSolverCH = createSubSolver(dim+1, "laplace_");
240
241
      laplaceMatrixSolverCH->fillPetscMatrix(&laplaceMatrixCH);

Sebastian Aland's avatar
Sebastian Aland committed
242
243
      // === matrix (-delta*K) ===
      deltaKMatrix.assembleOperator(laplaceOp2);
244
      deltaKMatrixSolver = createSubSolver(dim+1, "laplace2_");
245
      deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);
Sebastian Aland's avatar
Sebastian Aland committed
246
    }
247

Sebastian Aland's avatar
Sebastian Aland committed
248
249
250
    // === Setup solver ===
    matShellContext.kspMassCH = massMatrixSolverCH->getSolver();
    matShellContext.kspLaplaceCH = laplaceMatrixSolverCH->getSolver();
251
252
    matShellContext.matMassCH = massMatrixSolverCH->getMatInterior();
    matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();
Sebastian Aland's avatar
Sebastian Aland committed
253
254
    matShellContext.eps = eps;
    matShellContext.delta = delta;
255

Sebastian Aland's avatar
Sebastian Aland committed
256
257
    petsc_helper::setSolver(matShellContext.kspMassCH, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
    petsc_helper::setSolver(matShellContext.kspLaplaceCH, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
258
259


Sebastian Aland's avatar
Sebastian Aland committed
260
    /// Init Navier-Stokes part
261

Sebastian Aland's avatar
Sebastian Aland committed
262
263
264
265
266
267
268
269
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
      Operator massOp(pressureFeSpace, pressureFeSpace);
      ZeroOrderTerm *massTerm = new Simple_ZOT;
      massOp.addTerm(massTerm);
      massMatrix.assembleOperator(massOp);
      delete massTerm;
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);
270
271
272
    matShellContext.kspMass = massMatrixSolver->getSolver();

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
Sebastian Aland's avatar
Sebastian Aland committed
273
274
275
276
      Operator laplaceOp(pressureFeSpace, pressureFeSpace);
      SecondOrderTerm *laplaceTerm = new Simple_SOT;
      laplaceOp.addTerm(laplaceTerm);
      laplaceMatrix.assembleOperator(laplaceOp);
277
      delete laplaceTerm;
Sebastian Aland's avatar
Sebastian Aland committed
278
279
    laplaceMatrixSolver = createSubSolver(pressureComponent, string("laplace_"));
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
280

Sebastian Aland's avatar
Sebastian Aland committed
281
    // === Create convection-diffusion operator ===
282

Sebastian Aland's avatar
Sebastian Aland committed
283
284
285
286
287
288
289
290
291
292
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
    DOFVector<double> vz(pressureFeSpace, "vz");
    DOFVector<double> vp(pressureFeSpace, "vp");
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
    if (dim == 3) vz.interpol(solution->getDOFVector(2));
    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    {
      Operator conDifOp(pressureFeSpace, pressureFeSpace);
293
294
295
      ZeroOrderTerm *conDif0 = NULL;
      SecondOrderTerm *conDif1 = NULL;
      FirstOrderTerm *conDif2 = NULL, *conDif3 = NULL, *conDif4 = NULL;
Sebastian Aland's avatar
Sebastian Aland committed
296
      vp.interpol(solution->getDOFVector(dim+2));
297

Sebastian Aland's avatar
Sebastian Aland committed
298
299
      densityFunctionTau = new LinearInterpolation(*rho1,*rho2,*invTau);
      viscosityFunction = new LinearInterpolation(*nu1,*nu2);
300
301
      densityFunction = new LinearInterpolation2(*rho1,*rho2);

Sebastian Aland's avatar
Sebastian Aland committed
302
      conDif0 = new VecAtQP_ZOT(&vp, densityFunctionTau);
303
      conDifOp.addTerm(conDif0);
Sebastian Aland's avatar
Sebastian Aland committed
304
      conDif1 = new VecAtQP_SOT(&vp, viscosityFunction );
305
      conDifOp.addTerm(conDif1);
Sebastian Aland's avatar
Sebastian Aland committed
306
307
      conDif2 = new Vec2AtQP_FOT(&vx, &vp, densityFunction, 0);
      conDifOp.addTerm(conDif2, GRD_PHI);
308

Sebastian Aland's avatar
Sebastian Aland committed
309
310
      conDif3 = new Vec2AtQP_FOT(&vy, &vp, densityFunction, 1);
      conDifOp.addTerm(conDif3, GRD_PHI);
311

Sebastian Aland's avatar
Sebastian Aland committed
312
313
314
315
      if (dim == 3) {
	conDif4 = new Vec2AtQP_FOT(&vz, &vp, densityFunction, 2);
	conDifOp.addTerm(conDif4, GRD_PHI);
      }
316

Sebastian Aland's avatar
Sebastian Aland committed
317
      conDifMatrix.assembleOperator(conDifOp);
318

Sebastian Aland's avatar
Sebastian Aland committed
319
320
321
322
      delete conDif0;
      delete conDif1;
      delete conDif2;
      delete conDif3;
323
      if (dim == 3) delete conDif4;
Sebastian Aland's avatar
Sebastian Aland committed
324
    }
325

Sebastian Aland's avatar
Sebastian Aland committed
326
327
    conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
    conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);
328
329
    matShellContext.matConDif = conDifMatrixSolver->getMatInterior();

330
331
332
    extractMatrixComponent(mat[0][0], dim, 1, 0, dim, &(matShellContext.matDiv));
    extractMatrixComponent(mat[0][0], 0, dim, dim, 1, &(matShellContext.matGrad));
    extractMatrixComponent(mat[0][0], 0, dim, 0, dim, &(matShellContext.velocityMat));
333

Sebastian Aland's avatar
Sebastian Aland committed
334
335
    ///erstelle kspVelocity
    KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspVelocity));
336
    petsc::ksp_set_operators(matShellContext.kspVelocity, matShellContext.velocityMat, matShellContext.velocityMat);
337

Sebastian Aland's avatar
Sebastian Aland committed
338
339
340
341
342
343
344
    ///regularisiere LaplaceMatrix
  if (regularizeLaplace)
  {
    PetscInt rows[1];
    rows[0]=0;
    MatZeroRows(laplaceMatrixSolver->getMatInterior(), 1, rows, 0, PETSC_NULL, PETSC_NULL);
    KSPCreate((meshDistributor->getMpiComm(0)), &(matShellContext.kspLaplace));
345
    petsc::ksp_set_operators(matShellContext.kspLaplace, laplaceMatrixSolver->getMatInterior(), laplaceMatrixSolver->getMatInterior());
Sebastian Aland's avatar
Sebastian Aland committed
346
347
348
349
350
  }
   else
   {  matShellContext.kspLaplace=laplaceMatrixSolver->getSolver();
     setConstantNullSpace(matShellContext.kspLaplace);
   }
351
352
353



Sebastian Aland's avatar
Sebastian Aland committed
354
    // === Setup solver ===
355

Sebastian Aland's avatar
Sebastian Aland committed
356
357
358
359
360
361
362
363
364
365
    switch (massSolutionMode) {
      case 0:
	petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
	break;
      case 1:
	petsc_helper::setSolverWithLu(matShellContext.kspMass, "mass_", KSPRICHARDSON,     PCLU, MATSOLVERMUMPS, 0.0, 1e-14, 1);
	break;
      default:
	ERROR_EXIT("No mass solution mode %d available!\n", massSolutionMode);
    }
366

Sebastian Aland's avatar
Sebastian Aland committed
367
368
369
370
371
372
373
374
375
376
    switch (laplaceSolutionMode) {
      case 0:
	petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
	break;
      case 1:
	petsc_helper::setSolverWithLu(matShellContext.kspLaplace, "mass_",   KSPRICHARDSON, PCLU, MATSOLVERMUMPS,   0.0, 1e-14, 1);
	break;
      default:
	ERROR_EXIT("No laplace solution mode %d available!\n", laplaceSolutionMode);
    }
377

Sebastian Aland's avatar
Sebastian Aland committed
378
    switch (velocitySolutionMode) {
379
      case 0:
Sebastian Aland's avatar
Sebastian Aland committed
380
381
382
383
384
385
386
387
	petsc_helper::setSolver(matShellContext.kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
	break;
      case 1:
	petsc_helper::setSolverWithLu(matShellContext.kspVelocity, "", KSPPREONLY,   PCLU, MATSOLVERMUMPS , 0.0, 1e-14, 1);
	break;
      default:
	ERROR_EXIT("No velocity solution mode %d available!\n", velocitySolutionMode);
    }
388
389


Sebastian Aland's avatar
Sebastian Aland committed
390
    PCSetFromOptions(pc);
391
392

    MSG("Setup of Cahn-Hilliard preconditioner needed %.5f seconds\n",
Sebastian Aland's avatar
Sebastian Aland committed
393
394
	MPI::Wtime() - wtime);
     }
395
396


Sebastian Aland's avatar
Sebastian Aland committed
397
398
399
     void PetscSolverNSCH::exitPreconditioner(PC pc)
     {
       FUNCNAME("PetscSolverNSCH::exitPreconditioner()");
400

Sebastian Aland's avatar
Sebastian Aland committed
401
402
403
       massMatrixSolverCH->destroyMatrixData();
       massMatrixSolverCH->destroyVectorData();
       laplaceMatrixSolverCH->destroyMatrixData();
404
       laplaceMatrixSolverCH->destroyVectorData();
Sebastian Aland's avatar
Sebastian Aland committed
405
406
       deltaKMatrixSolver->destroyMatrixData();
       deltaKMatrixSolver->destroyVectorData();
407
408


Sebastian Aland's avatar
Sebastian Aland committed
409
410
411
       massMatrixSolverCH->destroyVectorData();
       laplaceMatrixSolverCH->destroyVectorData();
       deltaKMatrixSolver->destroyVectorData();
412

Sebastian Aland's avatar
Sebastian Aland committed
413
       delete massMatrixSolverCH;
414
       massMatrixSolverCH = NULL;
415

Sebastian Aland's avatar
Sebastian Aland committed
416
       delete laplaceMatrixSolverCH;
417
       laplaceMatrixSolverCH = NULL;
418

Sebastian Aland's avatar
Sebastian Aland committed
419
       delete deltaKMatrixSolver;
420
       deltaKMatrixSolver = NULL;
421

Sebastian Aland's avatar
Sebastian Aland committed
422
423
424
425
426
427
       massMatrixSolver->destroyMatrixData();
       massMatrixSolver->destroyVectorData();
       laplaceMatrixSolver->destroyMatrixData();
       laplaceMatrixSolver->destroyVectorData();
       conDifMatrixSolver->destroyMatrixData();
       conDifMatrixSolver->destroyVectorData();
428

Sebastian Aland's avatar
Sebastian Aland committed
429
430
431
       massMatrixSolver->destroyVectorData();
       laplaceMatrixSolver->destroyVectorData();
       conDifMatrixSolver->destroyVectorData();
432
433


Sebastian Aland's avatar
Sebastian Aland committed
434
       delete massMatrixSolver;
435
       massMatrixSolver = NULL;
436

Sebastian Aland's avatar
Sebastian Aland committed
437
       delete laplaceMatrixSolver;
438
       laplaceMatrixSolver = NULL;
439

Sebastian Aland's avatar
Sebastian Aland committed
440
       delete conDifMatrixSolver;
441
       conDifMatrixSolver = NULL;
442

Sebastian Aland's avatar
Sebastian Aland committed
443
444
445
       KSPDestroy(&(matShellContext.kspVelocity));
       if (regularizeLaplace)
        KSPDestroy(&(matShellContext.kspLaplace));
446

Sebastian Aland's avatar
Sebastian Aland committed
447
448
449
       MatDestroy(&(matShellContext.matGrad));
       MatDestroy(&(matShellContext.matDiv));
       MatDestroy(&(matShellContext.velocityMat));
450

Sebastian Aland's avatar
Sebastian Aland committed
451
452
453
454
       delete densityFunction;
       delete densityFunctionTau;
       delete viscosityFunction;
     }
455
  } }
456