PetscSolverCahnHilliard.cc 8.22 KB
 Praetorius, Simon committed Dec 19, 2012 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ``````// // 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 "PetscSolverCahnHilliard.h" #include "parallel/PetscHelper.h" #include "parallel/PetscSolverGlobalMatrix.h" namespace AMDiS { `````` Sebastian Aland committed Feb 05, 2013 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 `````` class Multiplier: public AbstractFunction { public: Multiplier(double d_) : AbstractFunction(2),d(d_) {}; double operator()(const double& x) const { return d*x; } protected: double d; }; `````` Praetorius, Simon committed Dec 19, 2012 34 `````` using namespace std; `````` Sebastian Aland committed Feb 05, 2013 35 `````` `````` Praetorius, Simon committed Dec 19, 2012 36 37 38 39 `````` /// solve Cahn-Hilliard Preconditioner PetscErrorCode pcChShell(PC pc, Vec b, Vec x) // solve Px=b {FUNCNAME("PCApply()"); `````` Sebastian Aland committed Feb 05, 2013 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 `````` void *ctx; PCShellGetContext(pc, &ctx); CahnHilliardData* data = static_cast(ctx); Vec b1, b2, x1, x2; VecNestGetSubVec(b, 0, &b1); VecNestGetSubVec(b, 1, &b2); VecNestGetSubVec(x, 0, &x1); VecNestGetSubVec(x, 1, &x2); Vec y1, y2; VecDuplicate(b1, &y1); VecDuplicate(b2, &y2); // MatGetDiagonal(data->matM, y2); // VecReciprocal(y2); // VecPointwiseMult(y1, y2, b1); KSPSolve(data->kspMass, b1, y1); // M*y1 = b1 MatMultAdd(data->matMinusDeltaK, y1, b2, x1); // -> x1 := b2-delta*K*y1 KSPSolve(data->kspLaplace, x1, y2); // (M+eps*sqrt(delta))*y2 = x1 MatMult(data->matM, y2, x1); // x1 := M*y2 KSPSolve(data->kspLaplace, x1, x2); // (M+eps*sqrt(delta))*x2 = x1 double factor = (*data->eps)/sqrt(*data->delta); VecCopy(x2, x1); // x1 := x2 VecAXPBYPCZ(x1, 1.0, factor, -factor, y1, y2); // x1 = 1*y1 + factor*y2 - factor*x1 VecDestroy(&y1); VecDestroy(&y2); `````` Praetorius, Simon committed Dec 19, 2012 72 73 `````` } `````` Sebastian Aland committed Feb 05, 2013 74 `````` `````` Praetorius, Simon committed Dec 19, 2012 75 `````` PetscSolverCahnHilliard::PetscSolverCahnHilliard(string name, double *epsPtr, double *deltaPtr) `````` Sebastian Aland committed Feb 05, 2013 76 77 78 79 80 81 82 83 84 85 86 `````` : PetscSolverGlobalBlockMatrix(name), massMatrixSolver(NULL), laplaceMatrixSolver(NULL), deltaKMatrixSolver(NULL), useOldInitialGuess(false), eps(epsPtr), delta(deltaPtr) { Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess); } `````` Praetorius, Simon committed Dec 19, 2012 87 88 89 90 `````` void PetscSolverCahnHilliard::solvePetscMatrix(SystemVector &vec, AdaptInfo *adaptInfo) { FUNCNAME("PetscSolverCahnHilliard::solvePetscMatrix()"); `````` Sebastian Aland committed Feb 05, 2013 91 92 93 94 `````` /* if (useOldInitialGuess) { if (getVecSolInterior()) {VecSet(getVecSolInterior(), 0.0); `````` Praetorius, Simon committed Dec 19, 2012 95 96 97 98 `````` for (int i = 0; i < solution->getSize(); i++) { Vec tmp; `````` Sebastian Aland committed Feb 05, 2013 99 `````` VecNestGetSubVec(getVecSolInterior(), i, &tmp); `````` Praetorius, Simon committed Dec 19, 2012 100 101 102 103 104 105 `````` setDofVector(tmp, solution->getDOFVector(i)); } vecSolAssembly(); KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE); } `````` Sebastian Aland committed Feb 05, 2013 106 107 `````` KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE); }*/ `````` Praetorius, Simon committed Dec 19, 2012 108 109 `````` PetscSolverGlobalBlockMatrix::solvePetscMatrix(vec, adaptInfo); } `````` Sebastian Aland committed Feb 05, 2013 110 `````` `````` Praetorius, Simon committed Dec 19, 2012 111 112 113 `````` void PetscSolverCahnHilliard::initSolver(KSP &ksp) { FUNCNAME("PetscSolverCahnHilliard::initSolver()"); `````` Sebastian Aland committed Feb 05, 2013 114 `````` `````` Praetorius, Simon committed Dec 19, 2012 115 116 117 118 119 `````` // Create FGMRES based outer solver KSPCreate(meshDistributor->getMpiComm(0), &ksp); KSPSetOperators(ksp, getMatInterior(), getMatInterior(), SAME_NONZERO_PATTERN); KSPMonitorSet(ksp, KSPMonitorTrueResidualNorm, PETSC_NULL, PETSC_NULL); petsc_helper::setSolver(ksp, "ch_", KSPFGMRES, PCNONE, 1e-6, 1e-8, 1000); `````` Sebastian Aland committed Feb 05, 2013 120 121 `````` KSPSetFromOptions(ksp); `````` Praetorius, Simon committed Dec 19, 2012 122 `````` } `````` Sebastian Aland committed Feb 05, 2013 123 124 `````` `````` Praetorius, Simon committed Dec 19, 2012 125 126 127 128 `````` void PetscSolverCahnHilliard::initPreconditioner(PC pc) { FUNCNAME("PetscSolverCahnHilliard::initPreconditioner()"); MSG("PetscSolverCahnHilliard::initPreconditioner()\n"); `````` Sebastian Aland committed Feb 05, 2013 129 130 131 `````` // KSPSetUp(kspInterior); `````` Praetorius, Simon committed Dec 19, 2012 132 133 134 `````` PCSetType(pc, PCSHELL); PCShellSetApply(pc, pcChShell); PCShellSetContext(pc, &matShellContext); `````` Sebastian Aland committed Feb 05, 2013 135 `````` `````` Praetorius, Simon committed Dec 19, 2012 136 137 138 139 140 141 `````` const FiniteElemSpace *feSpace = componentSpaces[0]; // === matrix M === DOFMatrix laplaceMatrix(feSpace, feSpace); Operator laplaceOp(feSpace, feSpace); `````` Sebastian Aland committed Feb 05, 2013 142 143 144 `````` DOFMatrix massMatrix(feSpace, feSpace); Operator massOp(feSpace, feSpace); `````` Praetorius, Simon committed Dec 19, 2012 145 146 `````` DOFMatrix deltaKMatrix(feSpace, feSpace); Operator laplaceOp2(feSpace, feSpace); `````` Sebastian Aland committed Feb 05, 2013 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 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 `````` if (phase) { VecAtQP_ZOT zot(phase, NULL); massOp.addTerm(&zot); laplaceOp.addTerm(&zot); // M VecAtQP_SOT sot2(phase, new Multiplier((*eps)*sqrt(*delta))); laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K VecAtQP_SOT sot(phase, new Multiplier(-(*delta))); laplaceOp2.addTerm(&sot); // -delta*K massMatrix.assembleOperator(massOp); massMatrixSolver = createSubSolver(0, "mass_"); massMatrixSolver->fillPetscMatrix(&massMatrix); // === matrix (M + eps*sqrt(delta)*K) === laplaceMatrix.assembleOperator(laplaceOp); laplaceMatrixSolver = createSubSolver(0, "laplace_"); laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix); // === matrix (-delta*K) === deltaKMatrix.assembleOperator(laplaceOp2); deltaKMatrixSolver = createSubSolver(0, "laplace2_"); deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix); } else { Simple_ZOT zot; massOp.addTerm(&zot); laplaceOp.addTerm(&zot); // M Simple_SOT sot2((*eps)*sqrt(*delta)); laplaceOp.addTerm(&sot2); // eps*sqrt(delta)*K Simple_SOT sot(-(*delta)); laplaceOp2.addTerm(&sot); // -delta*K massMatrix.assembleOperator(massOp); massMatrixSolver = createSubSolver(0, "mass_"); massMatrixSolver->fillPetscMatrix(&massMatrix); // === matrix (M + eps*sqrt(delta)*K) === laplaceMatrix.assembleOperator(laplaceOp); laplaceMatrixSolver = createSubSolver(0, "laplace_"); laplaceMatrixSolver->fillPetscMatrix(&laplaceMatrix); // === matrix (-delta*K) === deltaKMatrix.assembleOperator(laplaceOp2); deltaKMatrixSolver = createSubSolver(0, "laplace2_"); deltaKMatrixSolver->fillPetscMatrix(&deltaKMatrix); } `````` Praetorius, Simon committed Dec 19, 2012 202 203 204 205 206 207 208 209 210 211 `````` // === Setup solver === matShellContext.kspMass = massMatrixSolver->getSolver(); matShellContext.kspLaplace = laplaceMatrixSolver->getSolver(); matShellContext.matM = massMatrixSolver->getMatInterior(); matShellContext.matMinusDeltaK = deltaKMatrixSolver->getMatInterior(); matShellContext.eps = eps; matShellContext.delta = delta; matShellContext.mpiCommGlobal= &(meshDistributor->getMpiComm(0)); `````` Sebastian Aland committed Feb 05, 2013 212 `````` `````` Praetorius, Simon committed Dec 19, 2012 213 `````` petsc_helper::setSolver(matShellContext.kspMass, "", KSPCG, PCJACOBI, 0.0, 1e-14, 2); `````` Sebastian Aland committed Feb 05, 2013 214 215 216 217 218 219 220 `````` // petsc_helper::setSolver(matShellContext.kspMass, "mass_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1); // { // PC pc; // KSPGetPC(matShellContext.kspMass, &pc); // PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS); // } `````` Praetorius, Simon committed Dec 19, 2012 221 `````` petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCHYPRE, 0.0, 1e-14, 1); `````` Sebastian Aland committed Feb 05, 2013 222 223 224 225 226 227 228 `````` // petsc_helper::setSolver(matShellContext.kspLaplace, "laplace_", KSPRICHARDSON, PCLU, 0.0, 1e-14, 1); // { // PC pc; // KSPGetPC(matShellContext.kspLaplace, &pc); // PCFactorSetMatSolverPackage(pc, MATSOLVERMUMPS); // } PCSetFromOptions(pc); `````` Praetorius, Simon committed Dec 19, 2012 229 `````` } `````` Sebastian Aland committed Feb 05, 2013 230 `````` `````` Praetorius, Simon committed Dec 19, 2012 231 232 233 234 235 `````` PetscSolver* PetscSolverCahnHilliard::createSubSolver(int component, string kspPrefix) { FUNCNAME("PetscSolverCahnHilliard::createSubSolver()"); `````` Sebastian Aland committed Feb 05, 2013 236 `````` `````` Praetorius, Simon committed Dec 19, 2012 237 238 `````` vector fe; fe.push_back(componentSpaces[component]); `````` Sebastian Aland committed Feb 05, 2013 239 `````` `````` Praetorius, Simon committed Dec 19, 2012 240 241 242 243 `````` PetscSolver* subSolver = new PetscSolverGlobalMatrix(""); subSolver->setKspPrefix(kspPrefix.c_str()); subSolver->setMeshDistributor(meshDistributor, 0); subSolver->init(fe, fe); `````` Sebastian Aland committed Feb 05, 2013 244 `````` `````` Praetorius, Simon committed Dec 19, 2012 245 246 247 `````` ParallelDofMapping &subDofMap = subSolver->getDofMap(); subDofMap[0] = dofMap[component]; subDofMap.update(); `````` Sebastian Aland committed Feb 05, 2013 248 `````` `````` Praetorius, Simon committed Dec 19, 2012 249 250 251 252 253 254 255 256 257 258 `````` return subSolver; } void PetscSolverCahnHilliard::exitPreconditioner(PC pc) { FUNCNAME("PetscSolverNavierStokes::exitPreconditioner()"); massMatrixSolver->destroyMatrixData(); laplaceMatrixSolver->destroyMatrixData(); deltaKMatrixSolver->destroyMatrixData(); `````` Sebastian Aland committed Feb 05, 2013 259 `````` `````` Praetorius, Simon committed Dec 19, 2012 260 261 262 `````` massMatrixSolver->destroyVectorData(); laplaceMatrixSolver->destroyVectorData(); deltaKMatrixSolver->destroyVectorData(); `````` Sebastian Aland committed Feb 05, 2013 263 `````` `````` Praetorius, Simon committed Dec 19, 2012 264 265 266 `````` delete massMatrixSolver; massMatrixSolver = NULL; `````` Sebastian Aland committed Feb 05, 2013 267 `````` `````` Praetorius, Simon committed Dec 19, 2012 268 269 `````` delete laplaceMatrixSolver; laplaceMatrixSolver = NULL; `````` Sebastian Aland committed Feb 05, 2013 270 `````` `````` Praetorius, Simon committed Dec 19, 2012 271 272 273 `````` delete deltaKMatrixSolver; deltaKMatrixSolver = NULL; } `````` Sebastian Aland committed Feb 05, 2013 274 275 `````` `````` Praetorius, Simon committed Dec 19, 2012 276 ``}``