Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

PetscSolverNavierStokes.cc 6.92 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
//
// 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"
15
#include "parallel/PetscHelper.h"
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
16
#include "TransformDOF.h"
17
18
19
20
21

namespace AMDiS {

  using namespace std;

22

23
  PetscErrorCode pcSchurShell(PC pc, Vec x, Vec y)
Thomas Witkowski's avatar
Thomas Witkowski committed
24
25
26
  {
    void *ctx;
    PCShellGetContext(pc, &ctx);
27
    NavierStokesSchurData* data = static_cast<NavierStokesSchurData*>(ctx);
Thomas Witkowski's avatar
Thomas Witkowski committed
28

29
30
31
32
33
34
35
36
37
38
    // Project out constant null space
    {
      int vecSize;
      VecGetSize(x, &vecSize);
      PetscScalar vecSum;
      VecSum(x, &vecSum);
      vecSum = vecSum / static_cast<PetscScalar>(-1.0 * vecSize);
      VecShift(x, vecSum); 
    }

39
40
41
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
42
  }
43
  
44
45
46
47

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
48
      pressureNullSpace(true),
49
      massMatrixSolver(NULL),
50
51
52
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
53
54
      solution(NULL),
      phase(NULL)
55
  {
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
56
    Parameters::get(initFileStr + "->navierstokes->pressure component", 
57
58
59
		    pressureComponent);
    TEST_EXIT(pressureComponent >= 0)
      ("For using PETSc stokes solver you must define a pressure component!\n");
60
61
62
63

    Parameters::get(initFileStr + "->navierstokes->pressure null space",
		    pressureNullSpace);
    TEST_EXIT(pressureNullSpace)("This is not yet tested, may be wrong!\n");
64
65
66
  }


67
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
68
69
70
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

71
    // Create FGMRES based outer solver
72
    KSPCreate(mpiCommGlobal, &ksp);
73
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
74
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
75
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 10000);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
76
    
77
    // Create null space information.
78
79
    if (pressureNullSpace)
      setConstantNullSpace(ksp, pressureComponent, true);
80
81
82
83
84
85
86
  }


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

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
87
88
89
90
    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
    TEST_EXIT(solution)("solution pointer not set!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
91
92
    int dim = componentSpaces[pressureComponent]->getMesh()->getDim();

93
94
95
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
Thomas Witkowski's avatar
Thomas Witkowski committed
96
97
    if (dim == 3)
      velocityComponents.push_back(2);
98

99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
    createFieldSplit(pc, "velocity", velocityComponents);
    createFieldSplit(pc, "pressure", pressureComponent);

    KSPSetUp(kspInterior);

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

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

    KSP kspVelocity = subKsp[0];
    KSP kspSchur = subKsp[1];
    PetscFree(subKsp);

118
    petsc_helper::setSolver(kspVelocity, "", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
119
120

    KSPSetType(kspSchur, KSPPREONLY);
121
    PC pcSub;
122
123
124
125
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
126

127
128
    if (pressureNullSpace)
      setConstantNullSpace(kspSchur);
129

130
131

    // === Mass matrix solver ===
132
133
134
135

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
136
    if (!phase)
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
137
      massOp.addTerm(new Simple_ZOT);
138
139
    else
      massOp.addTerm(new VecAtQP_ZOT(phase, &idFct));
140
141
142
143
144
145
146
147
148
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
149
150
151
152
    if (!phase)
      laplaceOp.addTerm(new Simple_SOT);
    else
      laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct));
153
154
155
156
157
158
159
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

160
161
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
Thomas Witkowski committed
162
    DOFVector<double> vz(pressureFeSpace, "vz");
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
163
    DOFVector<double> vp(pressureFeSpace, "vp");
164
165
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
Thomas Witkowski's avatar
Thomas Witkowski committed
166
167
    if (dim == 3)
      vz.interpol(solution->getDOFVector(2));
168
169
170

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
171
172
173
174
175
176
177
178
179
180
181

    if (!phase) {
      MSG("INIT WITHOUT PHASE!\n");
      Simple_ZOT *conDif0 = new Simple_ZOT(*invTau);
      conDifOp.addTerm(conDif0);
      Simple_SOT *conDif1 = new Simple_SOT(*nu);
      conDifOp.addTerm(conDif1);
      VecAtQP_FOT *conDif2 = new VecAtQP_FOT(&vx, &idFct, 0);
      conDifOp.addTerm(conDif2, GRD_PHI);
      VecAtQP_FOT *conDif3 = new VecAtQP_FOT(&vy, &idFct, 1);
      conDifOp.addTerm(conDif3, GRD_PHI);
Thomas Witkowski's avatar
Thomas Witkowski committed
182
183
184
185
186

      if (dim == 3) {
	VecAtQP_FOT *conDif4 = new VecAtQP_FOT(&vz, &idFct, 2);	
	conDifOp.addTerm(conDif4, GRD_PHI);
      }
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
    } else {
      MSG("INIT WITH PHASE!\n");
      
      vp.interpol(phase);
      VecAtQP_ZOT *conDif0 = new VecAtQP_ZOT(&vp, new MultConstFct(*invTau));
      conDifOp.addTerm(conDif0);

      VecAtQP_SOT *conDif1 = new VecAtQP_SOT(&vp, new MultConstFct(*nu));
      conDifOp.addTerm(conDif1);

      Vec2AtQP_FOT *conDif2 = new Vec2AtQP_FOT(&vx, &vp, new Multiplier3(), 0);
      conDifOp.addTerm(conDif2, GRD_PHI);

      Vec2AtQP_FOT *conDif3 = new Vec2AtQP_FOT(&vy, &vp, new Multiplier3(), 1);
      conDifOp.addTerm(conDif3, GRD_PHI);
Thomas Witkowski's avatar
Thomas Witkowski committed
202
203
204
205
206

      if (dim == 3) {
	Vec2AtQP_FOT *conDif4 = new Vec2AtQP_FOT(&vz, &vp, new Multiplier3(), 2);
	conDifOp.addTerm(conDif4, GRD_PHI);
      }
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
207
    }
208
209
210
211
212
213

    conDifMatrix.assembleOperator(conDifOp);

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

214
215

    // === Setup solver ===
216
217
218
219
220

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

Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
221
    petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
222
    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
Thomas Witkowski's avatar
blub    
Thomas Witkowski committed
223

224
    setConstantNullSpace(matShellContext.kspLaplace);
225
226
  }

227

228
}