PetscSolverCahnHilliard_DD.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 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
 ******************************************************************************/



#include "PetscSolverCahnHilliard_DD.h"
#include "parallel/PetscHelper.h"
#include "TransformDOF.h"
#include "Expressions.h"

namespace AMDiS { namespace Parallel {


  using namespace std;


  PetscErrorCode pcChShell_dd_b(PC pc, Vec b, Vec x) // solve Px=b
  {FUNCNAME("PCApply()");
    void *ctx;
    PCShellGetContext(pc, &ctx);
    CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
39

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

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

48 49 50 51 52 53 54 55 56
    PetscFunctionReturn(0);
  }

  /// solve Cahn-Hilliard Preconditioner
  PetscErrorCode pcChSchurShell_dd(PC pc, Vec b, Vec x) // solve Px=b
  {FUNCNAME("PCApply()");
    void *ctx;
    PCShellGetContext(pc, &ctx);
    CahnHilliardData2* data = static_cast<CahnHilliardData2*>(ctx);
57

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

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

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

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

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

86 87
    PetscFunctionReturn(0);
  }
88 89


90 91

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


  void PetscSolverCahnHilliard_DD::initSolver(KSP &ksp)
  {
    FUNCNAME("PetscSolverCahnHilliard_DD::initSolver()");

    // Create FGMRES based outer solver
114

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

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

127 128 129 130 131 132 133 134 135
  }


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

    MPI::COMM_WORLD.Barrier();
    double wtime = MPI::Wtime();
136 137


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

143 144 145 146 147 148 149 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 177
    vector<int> chPotentialComponent;
    chPotentialComponent.push_back(0);
    vector<int> chSchurComponent;
    chSchurComponent.push_back(1);

    PCSetType(pc, PCFIELDSPLIT);
    PCFieldSplitSetType(pc, PC_COMPOSITE_SCHUR);
    PCFieldSplitSetSchurFactType(pc, PC_FIELDSPLIT_SCHUR_FACT_FULL);

    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_dd);
    PCShellSetContext(pcSub, &matShellContext);

    const FiniteElemSpace *feSpace = componentSpaces[0];
178

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

187 188 189 190 191 192 193 194
    // === Laplace matrix solver ===
    DOFMatrix laplaceMatrix(feSpace, feSpace);
    Operator laplaceOp(feSpace, feSpace);
    addSOT(laplaceOp, ((*eps)*std::sqrt(*delta)) * valueOf(phase));
    addSOT(laplaceOp, valueOf(phase));
    laplaceMatrix.assembleOperator(laplaceOp);
    laplaceMatrixSolver = createSubSolver(0, "MpK_");
    laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix);
195

196 197 198 199 200 201 202 203
    // === matrix (-delta*K) ===
    DOFMatrix deltaKMatrix(feSpace, feSpace);
    Operator laplaceOp2(feSpace, feSpace);
    addSOT(laplaceOp2, -(*delta) * valueOf(phase));
    deltaKMatrix.assembleOperator(laplaceOp2);
    deltaKMatrixSolver = createSubSolver(0, "laplace_");
    deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix);

204

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

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

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


  void PetscSolverCahnHilliard_DD::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverCahnHilliard_DD::exitPreconditioner()");

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

234 235 236 237

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

239 240 241 242 243
    delete massMatrixSolver;
    massMatrixSolver = NULL;

    delete laplaceMatrixSolver;
    laplaceMatrixSolver = NULL;
244

245 246 247
    delete deltaKMatrixSolver;
    deltaKMatrixSolver = NULL;

248

249 250 251 252 253
    if (tau) {
      delete delta;
    }
  }
} }