PetscSolverCahnHilliard2.cc 7.45 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 20 21 22 23
 ******************************************************************************/



#include "parallel/PetscSolverCahnHilliard2.h"
24
#include "parallel/PetscHelper.h"
25
#include "TransformDOF.h"
26

27
namespace AMDiS { namespace Parallel {
28

29

30
  using namespace std;
31 32 33


   PetscErrorCode pcChShell2b(PC pc, Vec b, Vec x) // solve Px=b
34 35 36 37
  {FUNCNAME("PCApply()");
    void *ctx;
    PCShellGetContext(pc, &ctx);
    CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
38

39 40 41
    Vec y1, y2;
    VecDuplicate(b, &y1);
    VecDuplicate(b, &y2);
42

43
    KSPSolve(data->kspMplusK, b, y1);
44
    MatMult(data->matMass, y1, y2);
45
    KSPSolve(data->kspMplusK, y2, x);
46

47
    PetscFunctionReturn(0);
48 49 50
  }

  /// solve Cahn-Hilliard Preconditioner
51
  PetscErrorCode pcChSchurShell(PC pc, Vec b, Vec x) // solve Px=b
52 53 54 55
  {FUNCNAME("PCApply()");
    void *ctx;
    PCShellGetContext(pc, &ctx);
    CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
56

57 58 59
    /// create S = M + eps^2*delta*K*D^(-1)*K where D=diag(M)
    Mat K, S;
    MatDuplicate(data->matMinusDeltaK, MAT_COPY_VALUES, &K);
60

61 62 63
    MatGetDiagonal(data->matMass, x);
    VecReciprocal(x);
    MatDiagonalScale(K, x, PETSC_NULL); 						// K' := -delta*D^(-1)*K
64
    MatMatMult(data->matMinusDeltaK, K, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &S); 		// S := -delta*K*K'
65
    MatAYPX(S, sqr(*data->eps)/(*data->delta), data->matMass, DIFFERENT_NONZERO_PATTERN); 	// S = eps^2/delta*S + M
66

67 68 69
    /// create new solver for S
    KSP kspS;
    KSPCreate(*data->mpiCommGlobal, &kspS);
70
    petsc::ksp_set_operators(kspS, S, S);
71 72 73 74 75 76 77
    petsc_helper::setSolver(kspS, "S_", KSPFGMRES, PCSHELL, 1e-6, 1e-8, 1);
    {
      PC pc;
      KSPGetPC(kspS, &pc);
      PCShellSetApply(pc, pcChShell2b);
      PCShellSetContext(pc, data);
    }
78

79
    KSPSolve(kspS, b, x); 				// S*x2 = x1
80

81 82 83
    MatDestroy(&S);
    MatDestroy(&K);
    KSPDestroy(&kspS);
84

85
    PetscFunctionReturn(0);
86
  }
87 88


89

90
  PetscSolverCahnHilliard2::PetscSolverCahnHilliard2(string name)
91 92
    : PetscSolverGlobalMatrix(name),
      useOldInitialGuess(false),
93
      laplaceSolutionMode(0),
94
      massMatrixSolver(NULL),
95
      laplaceMatrixSolver(NULL),
96 97 98 99 100 101
      deltaKMatrixSolver(NULL),
      eps(NULL),
      delta(NULL),
      tau(NULL),
      solution(NULL),
      phase(NULL)
102
  {
103 104
    Parameters::get(initFileStr + "->use old initial guess",
		    useOldInitialGuess);
105 106 107
  }


108 109 110 111 112
  void PetscSolverCahnHilliard2::initSolver(KSP &ksp)
  {
    FUNCNAME("PetscSolverCahnHilliard2::initSolver()");

    // Create FGMRES based outer solver
113

114
    MSG("CREATE POS 1: %p\n", &ksp);
115
    KSPCreate(domainComm, &ksp);
116
    petsc::ksp_set_operators(ksp, getMatInterior(), getMatInterior());
117
    if (getInfo() >= 10)
118
      petsc::ksp_monitor_set(ksp, KSPMonitorDefault);
119
    else if (getInfo() >= 20)
120
      petsc::ksp_monitor_set(ksp, KSPMonitorTrueResidualNorm);
121
    petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCNONE, getRelative(), getTolerance(), getMaxIterations());
122

123 124
    if (useOldInitialGuess)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
125

126 127 128 129 130 131 132
  }


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

133 134
    MPI::COMM_WORLD.Barrier();
    double wtime = MPI::Wtime();
135 136


137 138 139 140
    if (tau) {
      delta = new double;
      *delta = sqrt(*tau);
    }
141

142 143 144 145
    vector<int> chPotentialComponent;
    chPotentialComponent.push_back(0);
    vector<int> chSchurComponent;
    chSchurComponent.push_back(1);
146

147 148 149
    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);
150

151 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
    createFieldSplit(pc, "ch_potential", chPotentialComponent);
    createFieldSplit(pc, "ch_schur", chSchurComponent);
    PCSetFromOptions(pc);

    KSPSetUp(kspInterior);

    KSP *subKspCH;
    int nSubKspCH;
    PCFieldSplitGetSubKSP(pc, &nSubKspCH, &subKspCH);


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

    KSP kspChPotential = subKspCH[0];
    KSP kspChSchur = subKspCH[1];
    PetscFree(subKspCH);

    KSPSetType(kspChSchur, KSPPREONLY);
    PC pcSub;
    KSPGetPC(kspChSchur, &pcSub);
    PCSetType(pcSub, PCSHELL);
    PCShellSetApply(pcSub, pcChSchurShell);
    PCShellSetContext(pcSub, &matShellContext);

      const FiniteElemSpace *feSpace = componentSpaces[0];
177

178
    // === Mass matrix solver ===
179 180 181 182 183 184 185
    DOFMatrix massMatrix(feSpace, feSpace);
    Operator massOp(feSpace, feSpace);
    Simple_ZOT zot;
    massOp.addTerm(&zot);
    massMatrix.assembleOperator(massOp);
    massMatrixSolver = createSubSolver(0, "mass_");
    massMatrixSolver->fillPetscMatrix(&massMatrix);
186

187
    // === Laplace matrix solver ===
188 189 190 191 192 193 194 195
    DOFMatrix laplaceMatrix(feSpace, feSpace);
    Operator laplaceOp(feSpace, feSpace);
    laplaceOp.addTerm(&zot); // M
    Simple_SOT sot2((*eps)*sqrt(*delta));
    laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(0, "MpK_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
196

197 198 199 200 201 202 203 204 205
    // === matrix (-delta*K) ===
    DOFMatrix deltaKMatrix(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
    Simple_SOT sot(-(*delta));
    laplaceOp2.addTerm(&sot); // -delta*K
    deltaKMatrix.assembleOperator(laplaceOp2);
    deltaKMatrixSolver = createSubSolver(0, "laplace_");
    deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);

206

207 208
    // === Setup solver ===
    matShellContext.kspMplusK = laplaceMatrixSolver->getSolver();
209
    matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior();
210 211
    matShellContext.eps = eps;
    matShellContext.delta = delta;
212
    matShellContext.kspMass = massMatrixSolver->getSolver();
213
    matShellContext.matMass = massMatrixSolver->getMatInterior();
214 215 216 217
    matShellContext.mpiCommGlobal = &(meshDistributor->getMpiComm(0));

    petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPCG, PCJACOBI, 0.0, 1e-14, 2);
    petsc_helper::setSolver(matShellContext.kspMplusK, "MpK_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
218
    petsc_helper::setSolver(kspChPotential, "chPotential",  KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1);
219 220

    MSG("Setup of Cahn-Hilliard preconditioner needed %.5f seconds\n",
221
	MPI::Wtime() - wtime);
222 223 224 225 226
  }


  void PetscSolverCahnHilliard2::exitPreconditioner(PC pc)
  {
227 228
    FUNCNAME("PetscSolverCahnHilliard2::exitPreconditioner()");

229
    massMatrixSolver->destroyMatrixData();
230
    massMatrixSolver->destroyVectorData();
231
    laplaceMatrixSolver->destroyMatrixData();
232
    laplaceMatrixSolver->destroyVectorData();
233
    deltaKMatrixSolver->destroyMatrixData();
234
    deltaKMatrixSolver->destroyVectorData();
235

236 237 238 239

    massMatrixSolver->destroyVectorData();
    laplaceMatrixSolver->destroyVectorData();
    deltaKMatrixSolver->destroyVectorData();
240

241
    delete massMatrixSolver;
242
    massMatrixSolver = NULL;
243 244

    delete laplaceMatrixSolver;
245
    laplaceMatrixSolver = NULL;
246

247
    delete deltaKMatrixSolver;
248
    deltaKMatrixSolver = NULL;
249

250

251 252 253
    if (tau) {
      delete delta;
    }
254
  }
255
} }