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