1 #include <petscconvest.h> /*I "petscconvest.h" I*/ 2 #include <petscts.h> 3 #include <petscdmplex.h> 4 5 #include <petsc/private/petscconvestimpl.h> 6 7 static PetscErrorCode PetscConvEstSetTS_Private(PetscConvEst ce, PetscObject solver) { 8 PetscClassId id; 9 10 PetscFunctionBegin; 11 PetscCall(PetscObjectGetClassId(ce->solver, &id)); 12 PetscCheck(id == TS_CLASSID, PetscObjectComm((PetscObject)ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS"); 13 PetscCall(TSGetDM((TS)ce->solver, &ce->idm)); 14 PetscFunctionReturn(0); 15 } 16 17 static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u) { 18 PetscFunctionBegin; 19 PetscCall(TSComputeInitialCondition((TS)ce->solver, u)); 20 PetscFunctionReturn(0); 21 } 22 23 static PetscErrorCode PetscConvEstComputeErrorTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u, PetscReal errors[]) { 24 TS ts = (TS)ce->solver; 25 PetscErrorCode (*exactError)(TS, Vec, Vec); 26 27 PetscFunctionBegin; 28 PetscCall(TSGetComputeExactError(ts, &exactError)); 29 if (exactError) { 30 Vec e; 31 PetscInt f; 32 33 PetscCall(VecDuplicate(u, &e)); 34 PetscCall(TSComputeExactError(ts, u, e)); 35 PetscCall(VecNorm(e, NORM_2, errors)); 36 for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0]; 37 PetscCall(VecDestroy(&e)); 38 } else { 39 PetscReal t; 40 41 PetscCall(TSGetSolveTime(ts, &t)); 42 PetscCall(DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors)); 43 } 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode PetscConvEstGetConvRateTS_Temporal_Private(PetscConvEst ce, PetscReal alpha[]) { 48 TS ts = (TS)ce->solver; 49 Vec u, u0; 50 PetscReal *dt, *x, *y, slope, intercept; 51 PetscInt Ns, oNs, Nf = ce->Nf, f, Nr = ce->Nr, r; 52 53 PetscFunctionBegin; 54 PetscCall(PetscMalloc1(Nr + 1, &dt)); 55 PetscCall(TSGetTimeStep(ts, &dt[0])); 56 PetscCall(TSGetMaxSteps(ts, &oNs)); 57 PetscCall(TSGetSolution(ts, &u0)); 58 PetscCall(PetscObjectReference((PetscObject)u0)); 59 Ns = oNs; 60 for (r = 0; r <= Nr; ++r) { 61 if (r > 0) { 62 dt[r] = dt[r - 1] / ce->r; 63 Ns = PetscCeilReal(Ns * ce->r); 64 } 65 PetscCall(TSSetTime(ts, 0.0)); 66 PetscCall(TSSetStepNumber(ts, 0)); 67 PetscCall(TSSetTimeStep(ts, dt[r])); 68 PetscCall(TSSetMaxSteps(ts, Ns)); 69 PetscCall(TSGetSolution(ts, &u)); 70 PetscCall(PetscConvEstComputeInitialGuess(ce, r, NULL, u)); 71 PetscCall(TSSolve(ts, NULL)); 72 PetscCall(TSGetSolution(ts, &u)); 73 PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0)); 74 PetscCall(PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r * Nf])); 75 PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0)); 76 for (f = 0; f < Nf; ++f) { 77 ce->dofs[r * Nf + f] = 1.0 / dt[r]; 78 PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f])); 79 PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f])); 80 } 81 /* Monitor */ 82 PetscCall(PetscConvEstMonitorDefault(ce, r)); 83 } 84 /* Fit convergence rate */ 85 if (Nr) { 86 PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y)); 87 for (f = 0; f < Nf; ++f) { 88 for (r = 0; r <= Nr; ++r) { 89 x[r] = PetscLog10Real(dt[r]); 90 y[r] = PetscLog10Real(ce->errors[r * Nf + f]); 91 } 92 PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept)); 93 /* Since lg err = s lg dt + b */ 94 alpha[f] = slope; 95 } 96 PetscCall(PetscFree2(x, y)); 97 } 98 /* Reset solver */ 99 PetscCall(TSReset(ts)); 100 PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING)); 101 PetscCall(TSSetTime(ts, 0.0)); 102 PetscCall(TSSetStepNumber(ts, 0)); 103 PetscCall(TSSetTimeStep(ts, dt[0])); 104 PetscCall(TSSetMaxSteps(ts, oNs)); 105 PetscCall(TSSetSolution(ts, u0)); 106 PetscCall(PetscConvEstComputeInitialGuess(ce, 0, NULL, u0)); 107 PetscCall(VecDestroy(&u0)); 108 PetscCall(PetscFree(dt)); 109 PetscFunctionReturn(0); 110 } 111 112 static PetscErrorCode PetscConvEstGetConvRateTS_Spatial_Private(PetscConvEst ce, PetscReal alpha[]) { 113 TS ts = (TS)ce->solver; 114 Vec uInitial; 115 DM *dm; 116 PetscObject disc; 117 PetscReal *x, *y, slope, intercept; 118 PetscInt Nr = ce->Nr, r, Nf = ce->Nf, f, dim, oldlevel, oldnlev; 119 PetscErrorCode (*ifunc)(DM, PetscReal, Vec, Vec, Vec, void *); 120 PetscErrorCode (*ijac)(DM, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *); 121 PetscErrorCode (*rhsfunc)(DM, PetscReal, Vec, Vec, void *); 122 void *fctx, *jctx, *rctx; 123 void *ctx; 124 125 PetscFunctionBegin; 126 PetscCheck(ce->r == 2.0, PetscObjectComm((PetscObject)ce), PETSC_ERR_SUP, "Only refinement factor 2 is currently supported (not %g)", (double)ce->r); 127 PetscCall(DMGetDimension(ce->idm, &dim)); 128 PetscCall(DMGetApplicationContext(ce->idm, &ctx)); 129 PetscCall(DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE)); 130 PetscCall(DMGetRefineLevel(ce->idm, &oldlevel)); 131 PetscCall(PetscMalloc1((Nr + 1), &dm)); 132 PetscCall(TSGetSolution(ts, &uInitial)); 133 /* Loop over meshes */ 134 dm[0] = ce->idm; 135 for (r = 0; r <= Nr; ++r) { 136 Vec u; 137 #if defined(PETSC_USE_LOG) 138 PetscLogStage stage; 139 #endif 140 char stageName[PETSC_MAX_PATH_LEN]; 141 const char *dmname, *uname; 142 143 PetscCall(PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN - 1, "ConvEst Refinement Level %" PetscInt_FMT, r)); 144 #if defined(PETSC_USE_LOG) 145 PetscCall(PetscLogStageGetId(stageName, &stage)); 146 if (stage < 0) PetscCall(PetscLogStageRegister(stageName, &stage)); 147 #endif 148 PetscCall(PetscLogStagePush(stage)); 149 if (r > 0) { 150 if (!ce->noRefine) { 151 PetscCall(DMRefine(dm[r - 1], MPI_COMM_NULL, &dm[r])); 152 PetscCall(DMSetCoarseDM(dm[r], dm[r - 1])); 153 } else { 154 DM cdm, rcdm; 155 156 PetscCall(DMClone(dm[r - 1], &dm[r])); 157 PetscCall(DMCopyDisc(dm[r - 1], dm[r])); 158 PetscCall(DMGetCoordinateDM(dm[r - 1], &cdm)); 159 PetscCall(DMGetCoordinateDM(dm[r], &rcdm)); 160 PetscCall(DMCopyDisc(cdm, rcdm)); 161 } 162 PetscCall(DMCopyTransform(ce->idm, dm[r])); 163 PetscCall(PetscObjectGetName((PetscObject)dm[r - 1], &dmname)); 164 PetscCall(PetscObjectSetName((PetscObject)dm[r], dmname)); 165 for (f = 0; f <= Nf; ++f) { 166 PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *); 167 168 PetscCall(DMGetNullSpaceConstructor(dm[r - 1], f, &nspconstr)); 169 PetscCall(DMSetNullSpaceConstructor(dm[r], f, nspconstr)); 170 } 171 } 172 PetscCall(DMViewFromOptions(dm[r], NULL, "-conv_dm_view")); 173 /* Create solution */ 174 PetscCall(DMCreateGlobalVector(dm[r], &u)); 175 PetscCall(DMGetField(dm[r], 0, NULL, &disc)); 176 PetscCall(PetscObjectGetName(disc, &uname)); 177 PetscCall(PetscObjectSetName((PetscObject)u, uname)); 178 /* Setup solver */ 179 PetscCall(TSReset(ts)); 180 PetscCall(TSSetDM(ts, dm[r])); 181 PetscCall(DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx)); 182 if (r > 0) { 183 PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx)); 184 PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx)); 185 PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx)); 186 if (ifunc) PetscCall(DMTSSetIFunctionLocal(dm[r], ifunc, fctx)); 187 if (ijac) PetscCall(DMTSSetIJacobianLocal(dm[r], ijac, jctx)); 188 if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(dm[r], rhsfunc, rctx)); 189 } 190 PetscCall(TSSetTime(ts, 0.0)); 191 PetscCall(TSSetStepNumber(ts, 0)); 192 PetscCall(TSSetFromOptions(ts)); 193 PetscCall(TSSetSolution(ts, u)); 194 PetscCall(VecDestroy(&u)); 195 /* Create initial guess */ 196 PetscCall(TSGetSolution(ts, &u)); 197 PetscCall(PetscConvEstComputeInitialGuess(ce, r, dm[r], u)); 198 PetscCall(TSSolve(ts, NULL)); 199 PetscCall(TSGetSolution(ts, &u)); 200 PetscCall(PetscLogEventBegin(ce->event, ce, 0, 0, 0)); 201 PetscCall(PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r * Nf])); 202 PetscCall(PetscLogEventEnd(ce->event, ce, 0, 0, 0)); 203 for (f = 0; f < Nf; ++f) { 204 PetscSection s, fs; 205 PetscInt lsize; 206 207 /* Could use DMGetOutputDM() to add in Dirichlet dofs */ 208 PetscCall(DMGetLocalSection(dm[r], &s)); 209 PetscCall(PetscSectionGetField(s, f, &fs)); 210 PetscCall(PetscSectionGetConstrainedStorageSize(fs, &lsize)); 211 PetscCallMPI(MPI_Allreduce(&lsize, &ce->dofs[r * Nf + f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)ts))); 212 PetscCall(PetscLogEventSetDof(ce->event, f, ce->dofs[r * Nf + f])); 213 PetscCall(PetscLogEventSetError(ce->event, f, ce->errors[r * Nf + f])); 214 } 215 /* Monitor */ 216 PetscCall(PetscConvEstMonitorDefault(ce, r)); 217 if (!r) { 218 /* PCReset() does not wipe out the level structure */ 219 SNES snes; 220 KSP ksp; 221 PC pc; 222 223 PetscCall(TSGetSNES(ts, &snes)); 224 PetscCall(SNESGetKSP(snes, &ksp)); 225 PetscCall(KSPGetPC(ksp, &pc)); 226 PetscCall(PCMGGetLevels(pc, &oldnlev)); 227 } 228 /* Cleanup */ 229 PetscCall(PetscLogStagePop()); 230 } 231 PetscCall(DMTSGetIFunctionLocal(dm[r - 1], &ifunc, &fctx)); 232 PetscCall(DMTSGetIJacobianLocal(dm[r - 1], &ijac, &jctx)); 233 PetscCall(DMTSGetRHSFunctionLocal(dm[r - 1], &rhsfunc, &rctx)); 234 for (r = 1; r <= Nr; ++r) PetscCall(DMDestroy(&dm[r])); 235 /* Fit convergence rate */ 236 PetscCall(PetscMalloc2(Nr + 1, &x, Nr + 1, &y)); 237 for (f = 0; f < Nf; ++f) { 238 for (r = 0; r <= Nr; ++r) { 239 x[r] = PetscLog10Real(ce->dofs[r * Nf + f]); 240 y[r] = PetscLog10Real(ce->errors[r * Nf + f]); 241 } 242 PetscCall(PetscLinearRegression(Nr + 1, x, y, &slope, &intercept)); 243 /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */ 244 alpha[f] = -slope * dim; 245 } 246 PetscCall(PetscFree2(x, y)); 247 PetscCall(PetscFree(dm)); 248 /* Restore solver */ 249 PetscCall(TSReset(ts)); 250 { 251 /* PCReset() does not wipe out the level structure */ 252 SNES snes; 253 KSP ksp; 254 PC pc; 255 256 PetscCall(TSGetSNES(ts, &snes)); 257 PetscCall(SNESGetKSP(snes, &ksp)); 258 PetscCall(KSPGetPC(ksp, &pc)); 259 PetscCall(PCMGSetLevels(pc, oldnlev, NULL)); 260 PetscCall(DMSetRefineLevel(ce->idm, oldlevel)); /* The damn DMCoarsen() calls in PCMG can reset this */ 261 } 262 PetscCall(TSSetDM(ts, ce->idm)); 263 PetscCall(DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx)); 264 if (ifunc) PetscCall(DMTSSetIFunctionLocal(ce->idm, ifunc, fctx)); 265 if (ijac) PetscCall(DMTSSetIJacobianLocal(ce->idm, ijac, jctx)); 266 if (rhsfunc) PetscCall(DMTSSetRHSFunctionLocal(ce->idm, rhsfunc, rctx)); 267 PetscCall(TSSetConvergedReason(ts, TS_CONVERGED_ITERATING)); 268 PetscCall(TSSetTime(ts, 0.0)); 269 PetscCall(TSSetStepNumber(ts, 0)); 270 PetscCall(TSSetFromOptions(ts)); 271 PetscCall(TSSetSolution(ts, uInitial)); 272 PetscCall(PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial)); 273 PetscFunctionReturn(0); 274 } 275 276 PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal) { 277 PetscFunctionBegin; 278 PetscValidHeaderSpecific(ce, PETSC_OBJECT_CLASSID, 1); 279 ce->ops->setsolver = PetscConvEstSetTS_Private; 280 ce->ops->initguess = PetscConvEstInitGuessTS_Private; 281 ce->ops->computeerror = PetscConvEstComputeErrorTS_Private; 282 if (checkTemporal) { 283 ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private; 284 } else { 285 ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private; 286 } 287 PetscFunctionReturn(0); 288 } 289