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