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