PetscSolverNavierStokes.cc 6.67 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
    KSPSolve(data->kspLaplace, x, y);
    MatMult(data->matConDif, y, x);
    KSPSolve(data->kspMass, x, y);
32
  }
33
  
34 35 36 37 38

  PetscSolverNavierStokes::PetscSolverNavierStokes(string name)
    : PetscSolverGlobalMatrix(name),
      pressureComponent(-1),
      massMatrixSolver(NULL),
39 40 41
      laplaceMatrixSolver(NULL),
      nu(NULL),
      invTau(NULL),
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
42 43
      solution(NULL),
      phase(NULL)
44
  {
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
45
    Parameters::get(initFileStr + "->navierstokes->pressure component", 
46 47 48 49 50 51 52 53 54
		    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");
  }


55
  void PetscSolverNavierStokes::initSolver(KSP &ksp)
56 57 58
  {
    FUNCNAME("PetscSolverNavierStokes::initSolver()");

59
    // Create FGMRES based outer solver
60
    KSPCreate(mpiCommGlobal, &ksp);
61
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN);
62
    KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL);
63
    petsc_helper::setSolver(ksp, "ns_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 100);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
64
    
65
    // Create null space information.
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
66
    setConstantNullSpace(ksp, pressureComponent, true);
67 68 69 70 71 72 73
  }


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

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
74 75 76 77
    TEST_EXIT(nu)("nu pointer not set!\n");
    TEST_EXIT(invTau)("invtau pointer not set!\n");
    TEST_EXIT(solution)("solution pointer not set!\n");

78 79 80
    vector<int> velocityComponents;
    velocityComponents.push_back(0);
    velocityComponents.push_back(1);
81

82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
    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);

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
101
    petsc_helper::setSolver(kspVelocity, "velocity_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
102 103

    KSPSetType(kspSchur, KSPPREONLY);
104
    PC pcSub;
105 106 107 108
    KSPGetPC(kspSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcSchurShell);
    PCShellSetContext(pcSub, &matShellContext);
109

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
110
    setConstantNullSpace(kspSchur);
111

112 113

    // === Mass matrix solver ===
114 115 116 117

    const FiniteElemSpace *pressureFeSpace = componentSpaces[pressureComponent];
    DOFMatrix massMatrix(pressureFeSpace, pressureFeSpace);
    Operator massOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
118 119 120 121
    //    if (!phase)
      massOp.addTerm(new Simple_ZOT);
    // else
    //   massOp.addTerm(new VecAtQP_ZOT(phase, &idFct));
122 123 124 125 126
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(pressureComponent, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);


Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
127 128
    //    VtkWriter::writeFile(phase, "phase.vtu");

129 130 131 132
    // === Laplace matrix solver ===

    DOFMatrix laplaceMatrix(pressureFeSpace, pressureFeSpace);
    Operator laplaceOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
133 134 135 136
    if (!phase)
      laplaceOp.addTerm(new Simple_SOT);
    else
      laplaceOp.addTerm(new VecAtQP_SOT(phase, &idFct));
137 138 139 140 141 142 143
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(pressureComponent, "laplace_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);


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

144 145
    DOFVector<double> vx(pressureFeSpace, "vx");
    DOFVector<double> vy(pressureFeSpace, "vy");
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
146
    DOFVector<double> vp(pressureFeSpace, "vp");
147 148
    vx.interpol(solution->getDOFVector(0));
    vy.interpol(solution->getDOFVector(1));
149 150 151

    DOFMatrix conDifMatrix(pressureFeSpace, pressureFeSpace);
    Operator conDifOp(pressureFeSpace, pressureFeSpace);
Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178

    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);
    } 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);
    }
179 180 181 182 183 184

    conDifMatrix.assembleOperator(conDifOp);

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

185 186

    // === Setup solver ===
187 188 189 190 191

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

Thomas Witkowski's avatar
blub  
Thomas Witkowski committed
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
    petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
    //    petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);

    KSPSetType(matShellContext.kspLaplace, KSPRICHARDSON);
    KSPSetTolerances(matShellContext.kspLaplace, 0.0, 1e-14, PETSC_DEFAULT, 1);
    KSPSetOptionsPrefix(matShellContext.kspLaplace, "laplace_");
    KSPSetFromOptions(matShellContext.kspLaplace);

    MatNullSpace matNullSpace;
    MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullSpace);
    KSPSetNullSpace(matShellContext.kspLaplace, matNullSpace);
    MatNullSpaceDestroy(&matNullSpace);

    {
      PC pc;
      KSPGetPC(matShellContext.kspLaplace, &pc);
      PCSetType(pc, PCHYPRE);
      PCSetFromOptions(pc);
    }



    //    setConstantNullSpace(matShellContext.kspLaplace);

216 217
  }

218

219
}