PetscSolverNavierStokes.cc 7.14 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



#include "parallel/PetscSolverNavierStokes.h"

namespace AMDiS {

  using namespace std;

20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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

  int petscMultNavierStokesSchur(Mat mat, Vec x, Vec y)
  {
    void *ctx;
    MatShellGetContext(mat, &ctx);
    MatShellNavierStokesSchur* data = static_cast<MatShellNavierStokesSchur*>(ctx);

    switch (data->solverMode) {
    case 0:
      {
	Vec vec0, vec1;
	VecDuplicate(x, &vec0);
	MatGetVecs(data->A00, &vec1, PETSC_NULL);
	
	MatMult(data->A11, x, y);
	MatMult(data->A01, x, vec1);
	KSPSolve(data->kspVelocity, vec1, vec1);
	MatMult(data->A10, vec1, y);
	VecAYPX(y, -1.0, vec0);
	
	VecDestroy(&vec0);
	VecDestroy(&vec1);
      }
      break;
    case 1:
      {
	Vec vec0, vec1;
	VecDuplicate(x, &vec0);
	VecDuplicate(x, &vec1);

	KSPSolve(data->kspLaplace, x, vec0);
	MatMult(data->matConDif, vec0, vec1);
	KSPSolve(data->kspMass, vec1, y);

	VecDestroy(&vec0);
	VecDestroy(&vec1);
      }
      break;
    default:
      ERROR_EXIT("Wrong solver mode %d\n", data->solverMode);
    }
  }
62
  
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
      massMatrixSolver(NULL),
      laplaceMatrixSolver(NULL)
  {
    Parameters::get(initFileStr + "->stokes->pressure component", 
		    pressureComponent);

    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\n");

    TEST_EXIT(pressureComponent == 2)("Fixed for pressure component = 2\n");
  }


80
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
81
82
83
84
85
86
87
88
89
90
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

    KSPCreate(mpiCommGlobal, &ksp);
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), 
		    SAME_NONZERO_PATTERN); 
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPGMRES);
    KSPSetOptionsPrefix(ksp, "ns_");
    KSPSetFromOptions(ksp);
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);

    
    // === Create null space information. ===

    Vec nullSpaceBasis;
    VecDuplicate(getVecSolInterior(), &nullSpaceBasis);

    SystemVector basisVec("tmp", componentSpaces, componentSpaces.size(), true);
    basisVec.set(0.0);
    basisVec.getDOFVector(pressureComponent)->set(1.0);

    setDofVector(nullSpaceBasis, basisVec, true);
    VecAssemblyBegin(nullSpaceBasis);
    VecAssemblyEnd(nullSpaceBasis);
    VecNormalize(nullSpaceBasis, PETSC_NULL);

    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_FALSE, 1, &nullSpaceBasis, &matNullSpace);
110
111
112
113
114
115
116
117
118
119
120
121
  }


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

    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);

    PCSetType(pc, PCFIELDSPLIT);
122
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
123
124
125
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);
126
127
128
129
130
131
132

    KSPSetUp(kspInterior);

    KSP *subKsp;
    int nSubKsp;
    PCFieldSplitGetSubKSP(pc, &nSubKsp, &subKsp);

133
134
    TEST_EXIT(nSubKsp == 2)
      ("Wrong numer of KSPs inside of the fieldsplit preconditioner!\n");
135

136
137
    KSP kspVelocity = subKsp[0];
    KSP kspSchur = subKsp[1];
138
139
140
141
    PetscFree(subKsp);

    Mat A00, A01, A10, A11;
    PCFieldSplitGetSchurBlocks(pc, &A00, &A01, &A10, &A11);
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
    matShellContext.A00 = A00;
    matShellContext.A01 = A01;
    matShellContext.A10 = A10;
    matShellContext.A11 = A11;
    matShellContext.kspVelocity = kspVelocity;

    Mat matShell;
    MatDuplicate(A11, MAT_DO_NOT_COPY_VALUES, &matShell);
    MatSetType(matShell, MATSHELL);
    MatShellSetContext(matShell, &matShellContext);
    MatShellSetOperation(matShell, MATOP_MULT, (void(*)(void))petscMultNavierStokesSchur);
    MatSetUp(matShell);

    KSPSetType(kspVelocity, KSPPREONLY);
    PC pcSub;
    KSPGetPC(kspVelocity, &pcSub);
    PCSetType(pcSub, PCLU);
    PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);
160
161


162
163
164
165
166
    KSPSetType(kspSchur, KSPGMRES);
    KSPSetTolerances(kspSchur, 0, 1e-14, 1e+3, 1000);
    KSPSetOperators(kspSchur, matShell, matShell, SAME_NONZERO_PATTERN);
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCNONE);
167
168


169
    // === Mass matrix solver ===
170

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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
    MSG("Initialize mass matrix solver ...\n");

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
    Simple_ZOT zot;
    massOp.addTerm(&zot);
    massMatrix.assembleOperator(massOp);

    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);

    MSG("... OK!\n");


    // === Laplace matrix solver ===

    MSG("Initialize laplace matrix solver ...\n");

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
    Simple_SOT sot;
    laplaceOp.addTerm(&sot);
    laplaceMatrix.assembleOperator(laplaceOp);

    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);

    MSG("... OK!\n");


    // === Create convection-diffusion operator ===

    MSG("Initialize pressure convection-diffusion operator ...\n");

    double timestep = 1.0;
    Parameters::get("navierstokes->adapt->timestep", timestep);
    timestep = 1.0 / timestep;

    double nu = 0.01;
    Parameters::get("nu", nu);

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
    Simple_ZOT conDif0(timestep);
    conDifOp.addTerm(&conDif0);
    Simple_SOT conDif1(nu);
    conDifOp.addTerm(&conDif1);
    Vector_FOT conDif2(0);
    //    conDifOp.addTerm(&conDif2, GRD_PHI);
    Vector_FOT conDif3(1);
    //    conDifOp.addTerm(&conDif3, GRD_PHI);

    conDifMatrix.assembleOperator(conDifOp);

    conDifMatrixSolver = createSubSolver(pressureComponent, "condif_");
    conDifMatrixSolver->fillPetscMatrix(&conDifMatrix);

    MSG("... OK!\n");

    matShellContext.kspMass = massMatrixSolver->getSolver();
    matShellContext.kspLaplace = laplaceMatrixSolver->getSolver();
    matShellContext.matConDif = conDifMatrixSolver->getMatInterior();    

    KSPSetType(matShellContext.kspMass, KSPRICHARDSON);
    KSPSetTolerances(matShellContext.kspMass, 0, 1e-13, 1e+3, 1);
    KSPGetPC(matShellContext.kspMass, &pcSub);
238
239
    PCSetType(pcSub, PCLU);
    PCFactorSetMatSolverPackage(pcSub, MATSOLVERMUMPS);
240
241
242
243
244
245
246
247
248
    KSPSetFromOptions(matShellContext.kspMass);

    KSPSetType(matShellContext.kspLaplace, KSPCG);
    KSPSetTolerances(matShellContext.kspLaplace, 0, 1e-13, 1e+3, 10000);
    KSPGetPC(matShellContext.kspLaplace, &pcSub);
    PCSetType(pcSub, PCJACOBI);
    KSPSetFromOptions(matShellContext.kspLaplace);

    matShellContext.solverMode = 1;
249
250


251
252
    KSPSetFromOptions(kspVelocity);
    KSPSetFromOptions(kspSchur);
253
254
255
  }

}