Actual source code: tsconvest.c

      1: #include <petscconvest.h>
      2: #include <petscts.h>
      3: #include <petscdmplex.h>

      5: #include <petsc/private/petscconvestimpl.h>

      7: static PetscErrorCode PetscConvEstSetTS_Private(PetscConvEst ce, PetscObject solver)
      8: {
      9:   PetscClassId   id;

     13:   PetscObjectGetClassId(ce->solver, &id);
     14:   if (id != TS_CLASSID) SETERRQ(PetscObjectComm((PetscObject) ce), PETSC_ERR_ARG_WRONG, "Solver was not a TS");
     15:   TSGetDM((TS) ce->solver, &ce->idm);
     16:   return(0);
     17: }

     19: static PetscErrorCode PetscConvEstInitGuessTS_Private(PetscConvEst ce, PetscInt r, DM dm, Vec u)
     20: {

     24:   TSComputeInitialCondition((TS) ce->solver, u);
     25:   return(0);
     26: }

     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;

     35:   TSGetComputeExactError(ts, &exactError);
     36:   if (exactError) {
     37:     Vec      e;
     38:     PetscInt f;

     40:     VecDuplicate(u, &e);
     41:     TSComputeExactError(ts, u, e);
     42:     VecNorm(e, NORM_2, errors);
     43:     for (f = 1; f < ce->Nf; ++f) errors[f] = errors[0];
     44:     VecDestroy(&e);
     45:   } else {
     46:     PetscReal t;

     48:     TSGetSolveTime(ts, &t);
     49:     DMComputeL2FieldDiff(dm, t, ce->exactSol, ce->ctxs, u, errors);
     50:   }
     51:   return(0);
     52: }

     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;

     63:   TSGetSolution(ts, &u);
     64:   PetscMalloc1(Nr+1, &dt);
     65:   TSGetTimeStep(ts, &dt[0]);
     66:   TSGetMaxSteps(ts, &oNs);
     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:     TSSetTime(ts, 0.0);
     74:     TSSetStepNumber(ts, 0);
     75:     TSSetTimeStep(ts, dt[r]);
     76:     TSSetMaxSteps(ts, Ns);
     77:     PetscConvEstComputeInitialGuess(ce, r, NULL, u);
     78:     TSSolve(ts, u);
     79:     PetscLogEventBegin(ce->event, ce, 0, 0, 0);
     80:     PetscConvEstComputeError(ce, r, ce->idm, u, &ce->errors[r*Nf]);
     81:     PetscLogEventEnd(ce->event, ce, 0, 0, 0);
     82:     for (f = 0; f < Nf; ++f) {
     83:       PetscLogEventSetDof(ce->event, f, 1.0/dt[r]);
     84:       PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
     85:     }
     86:     /* Monitor */
     87:     PetscConvEstMonitorDefault(ce, r);
     88:   }
     89:   /* Fit convergence rate */
     90:   if (Nr) {
     91:     PetscMalloc2(Nr+1, &x, Nr+1, &y);
     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:       PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
     98:       /* Since lg err = s lg dt + b */
     99:       alpha[f] = slope;
    100:     }
    101:     PetscFree2(x, y);
    102:   }
    103:   /* Reset solver */
    104:   TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
    105:   TSSetTime(ts, 0.0);
    106:   TSSetStepNumber(ts, 0);
    107:   TSSetTimeStep(ts, dt[0]);
    108:   TSSetMaxSteps(ts, oNs);
    109:   PetscConvEstComputeInitialGuess(ce, 0, NULL, u);
    110:   PetscFree(dt);
    111:   return(0);
    112: }

    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;

    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:   DMGetDimension(ce->idm, &dim);
    128:   DMGetApplicationContext(ce->idm, &ctx);
    129:   DMPlexSetRefinementUniform(ce->idm, PETSC_TRUE);
    130:   DMGetRefineLevel(ce->idm, &oldlevel);
    131:   PetscMalloc1((Nr+1), &dm);
    132:   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;

    143:     PetscSNPrintf(stageName, PETSC_MAX_PATH_LEN-1, "ConvEst Refinement Level %D", r);
    144:     PetscLogStageRegister(stageName, &stage);
    145:     PetscLogStagePush(stage);
    146:     if (r > 0) {
    147:       DMRefine(dm[r-1], MPI_COMM_NULL, &dm[r]);
    148:       DMSetCoarseDM(dm[r], dm[r-1]);
    149:       DMCopyTransform(ce->idm, dm[r]);
    150:       PetscObjectGetName((PetscObject) dm[r-1], &dmname);
    151:       PetscObjectSetName((PetscObject) dm[r], dmname);
    152:       for (f = 0; f <= Nf; ++f) {
    153:         PetscErrorCode (*nspconstr)(DM, PetscInt, PetscInt, MatNullSpace *);

    155:         DMGetNullSpaceConstructor(dm[r-1], f, &nspconstr);
    156:         DMSetNullSpaceConstructor(dm[r],   f,  nspconstr);
    157:       }
    158:     }
    159:     DMViewFromOptions(dm[r], NULL, "-conv_dm_view");
    160:     /* Create solution */
    161:     DMCreateGlobalVector(dm[r], &u);
    162:     DMGetField(dm[r], 0, NULL, &disc);
    163:     PetscObjectGetName(disc, &uname);
    164:     PetscObjectSetName((PetscObject) u, uname);
    165:     /* Setup solver */
    166:     TSReset(ts);
    167:     TSSetDM(ts, dm[r]);
    168:     DMTSSetBoundaryLocal(dm[r], DMPlexTSComputeBoundary, ctx);
    169:     DMTSSetIFunctionLocal(dm[r], DMPlexTSComputeIFunctionFEM, ctx);
    170:     DMTSSetIJacobianLocal(dm[r], DMPlexTSComputeIJacobianFEM, ctx);
    171:     TSSetTime(ts, 0.0);
    172:     TSSetStepNumber(ts, 0);
    173:     TSSetFromOptions(ts);
    174:     /* Create initial guess */
    175:     PetscConvEstComputeInitialGuess(ce, r, dm[r], u);
    176:     TSSolve(ts, u);
    177:     PetscLogEventBegin(ce->event, ce, 0, 0, 0);
    178:     PetscConvEstComputeError(ce, r, dm[r], u, &ce->errors[r*Nf]);
    179:     PetscLogEventEnd(ce->event, ce, 0, 0, 0);
    180:     for (f = 0; f < Nf; ++f) {
    181:       PetscSection s, fs;
    182:       PetscInt     lsize;

    184:       /* Could use DMGetOutputDM() to add in Dirichlet dofs */
    185:       DMGetLocalSection(dm[r], &s);
    186:       PetscSectionGetField(s, f, &fs);
    187:       PetscSectionGetConstrainedStorageSize(fs, &lsize);
    188:       MPI_Allreduce(&lsize, &ce->dofs[r*Nf+f], 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject) ts));
    189:       PetscLogEventSetDof(ce->event, f, ce->dofs[r*Nf+f]);
    190:       PetscLogEventSetError(ce->event, f, ce->errors[r*Nf+f]);
    191:     }
    192:     /* Monitor */
    193:     PetscConvEstMonitorDefault(ce, r);
    194:     if (!r) {
    195:       /* PCReset() does not wipe out the level structure */
    196:       SNES snes;
    197:       KSP  ksp;
    198:       PC   pc;

    200:       TSGetSNES(ts, &snes);
    201:       SNESGetKSP(snes, &ksp);
    202:       KSPGetPC(ksp, &pc);
    203:       PCMGGetLevels(pc, &oldnlev);
    204:     }
    205:     /* Cleanup */
    206:     VecDestroy(&u);
    207:     PetscLogStagePop();
    208:   }
    209:   for (r = 1; r <= Nr; ++r) {
    210:     DMDestroy(&dm[r]);
    211:   }
    212:   /* Fit convergence rate */
    213:   PetscMalloc2(Nr+1, &x, Nr+1, &y);
    214:   for (f = 0; f < Nf; ++f) {
    215:     for (r = 0; r <= Nr; ++r) {
    216:       x[r] = PetscLog10Real(ce->dofs[r*Nf+f]);
    217:       y[r] = PetscLog10Real(ce->errors[r*Nf+f]);
    218:     }
    219:     PetscLinearRegression(Nr+1, x, y, &slope, &intercept);
    220:     /* Since h^{-dim} = N, lg err = s lg N + b = -s dim lg h + b */
    221:     alpha[f] = -slope * dim;
    222:   }
    223:   PetscFree2(x, y);
    224:   PetscFree(dm);
    225:   /* Restore solver */
    226:   TSReset(ts);
    227:   {
    228:     /* PCReset() does not wipe out the level structure */
    229:     SNES snes;
    230:     KSP  ksp;
    231:     PC   pc;

    233:     TSGetSNES(ts, &snes);
    234:     SNESGetKSP(snes, &ksp);
    235:     KSPGetPC(ksp, &pc);
    236:     PCMGSetLevels(pc, oldnlev, NULL);
    237:     DMSetRefineLevel(ce->idm, oldlevel); /* The damn DMCoarsen() calls in PCMG can reset this */
    238:   }
    239:   TSSetDM(ts, ce->idm);
    240:   DMTSSetBoundaryLocal(ce->idm, DMPlexTSComputeBoundary, ctx);
    241:   DMTSSetIFunctionLocal(ce->idm, DMPlexTSComputeIFunctionFEM, ctx);
    242:   DMTSSetIJacobianLocal(ce->idm, DMPlexTSComputeIJacobianFEM, ctx);
    243:   TSSetConvergedReason(ts, TS_CONVERGED_ITERATING);
    244:   TSSetTime(ts, 0.0);
    245:   TSSetStepNumber(ts, 0);
    246:   TSSetFromOptions(ts);
    247:   TSSetSolution(ts, uInitial);
    248:   PetscConvEstComputeInitialGuess(ce, 0, NULL, uInitial);
    249:   return(0);
    250: }

    252: PetscErrorCode PetscConvEstUseTS(PetscConvEst ce, PetscBool checkTemporal)
    253: {
    256:   ce->ops->setsolver     = PetscConvEstSetTS_Private;
    257:   ce->ops->initguess     = PetscConvEstInitGuessTS_Private;
    258:   ce->ops->computeerror  = PetscConvEstComputeErrorTS_Private;
    259:   if (checkTemporal) {
    260:     ce->ops->getconvrate = PetscConvEstGetConvRateTS_Temporal_Private;
    261:   } else {
    262:     ce->ops->getconvrate = PetscConvEstGetConvRateTS_Spatial_Private;
    263:   }
    264:   return(0);
    265: }