Actual source code: reconstruct.c

      1: #include <petsc/private/tshistoryimpl.h>
      2: #include <petscts.h>

      4: /* these two functions have been stolen from bdf.c */
      5: static inline void LagrangeBasisVals(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar L[])
      6: {
      7:   PetscInt k, j;
      8:   for (k = 0; k < n; k++) {
      9:     for (L[k] = 1, j = 0; j < n; j++) {
     10:       if (j != k) L[k] *= (t - T[j]) / (T[k] - T[j]);
     11:     }
     12:   }
     13: }

     15: static inline void LagrangeBasisDers(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar dL[])
     16: {
     17:   PetscInt k, j, i;
     18:   for (k = 0; k < n; k++) {
     19:     for (dL[k] = 0, j = 0; j < n; j++) {
     20:       if (j != k) {
     21:         PetscReal L = 1 / (T[k] - T[j]);
     22:         for (i = 0; i < n; i++) {
     23:           if (i != j && i != k) L *= (t - T[i]) / (T[k] - T[i]);
     24:         }
     25:         dL[k] += L;
     26:       }
     27:     }
     28:   }
     29: }

     31: static inline PetscInt LagrangeGetId(PetscReal t, PetscInt n, const PetscReal T[], const PetscBool Taken[])
     32: {
     33:   PetscInt _tid = 0;
     34:   while (_tid < n && PetscAbsReal(t - T[_tid]) > PETSC_SMALL) _tid++;
     35:   if (_tid < n && !Taken[_tid]) {
     36:     return _tid;
     37:   } else { /* we get back a negative id, where the maximum time is stored, since we use usually reconstruct backward in time */
     38:     PetscReal max    = PETSC_MIN_REAL;
     39:     PetscInt  maxloc = n;
     40:     _tid             = 0;
     41:     while (_tid < n) {
     42:       maxloc = (max < T[_tid] && !Taken[_tid]) ? (max = T[_tid], _tid) : maxloc;
     43:       _tid++;
     44:     }
     45:     return -maxloc - 1;
     46:   }
     47: }

     49: PetscErrorCode TSTrajectoryReconstruct_Private(TSTrajectory tj, TS ts, PetscReal t, Vec U, Vec Udot)
     50: {
     51:   TSHistory        tsh = tj->tsh;
     52:   const PetscReal *tshhist;
     53:   const PetscInt  *tshhist_id;
     54:   PetscInt         id, cnt, i, tshn;

     56:   PetscFunctionBegin;
     57:   PetscCall(TSHistoryGetLocFromTime(tsh, t, &id));
     58:   PetscCall(TSHistoryGetHistory(tsh, &tshn, &tshhist, &tshhist_id, NULL));
     59:   if (id == -1 || id == -tshn - 1) {
     60:     PetscReal t0 = tshn ? tshhist[0] : 0.0;
     61:     PetscReal tf = tshn ? tshhist[tshn - 1] : 0.0;
     62:     SETERRQ(PetscObjectComm((PetscObject)tj), PETSC_ERR_PLIB, "Requested time %g is outside the history interval [%g, %g] (%" PetscInt_FMT ")", (double)t, (double)t0, (double)tf, tshn);
     63:   }
     64:   if (tj->monitor) PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Reconstructing at time %g, order %" PetscInt_FMT "\n", (double)t, tj->lag.order));
     65:   if (!tj->lag.T) {
     66:     PetscInt o = tj->lag.order + 1;
     67:     PetscCall(PetscMalloc5(o, &tj->lag.L, o, &tj->lag.T, o, &tj->lag.WW, 2 * o, &tj->lag.TT, o, &tj->lag.TW));
     68:     for (i = 0; i < o; i++) tj->lag.T[i] = PETSC_MAX_REAL;
     69:     PetscCall(VecDuplicateVecs(U ? U : Udot, o, &tj->lag.W));
     70:   }
     71:   cnt = 0;
     72:   PetscCall(PetscArrayzero(tj->lag.TT, 2 * (tj->lag.order + 1)));
     73:   if (id < 0 || Udot) { /* populate snapshots for interpolation */
     74:     PetscInt s, nid = id < 0 ? -(id + 1) : id;

     76:     PetscInt up  = PetscMin(nid + tj->lag.order / 2 + 1, tshn);
     77:     PetscInt low = PetscMax(up - tj->lag.order - 1, 0);
     78:     up           = PetscMin(PetscMax(low + tj->lag.order + 1, up), tshn);
     79:     if (tj->monitor) PetscCall(PetscViewerASCIIPushTab(tj->monitor));

     81:     /* first see if we can reuse any */
     82:     for (s = up - 1; s >= low; s--) {
     83:       PetscReal t   = tshhist[s];
     84:       PetscInt  tid = LagrangeGetId(t, tj->lag.order + 1, tj->lag.T, tj->lag.TT);
     85:       if (tid < 0) continue;
     86:       if (tj->monitor) PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Reusing snapshot %" PetscInt_FMT ", step %" PetscInt_FMT ", time %g\n", tid, tshhist_id[s], (double)t));
     87:       tj->lag.TT[tid]                         = PETSC_TRUE;
     88:       tj->lag.WW[cnt]                         = tj->lag.W[tid];
     89:       tj->lag.TW[cnt]                         = t;
     90:       tj->lag.TT[tj->lag.order + 1 + s - low] = PETSC_TRUE; /* tell the next loop to skip it */
     91:       cnt++;
     92:     }

     94:     /* now load the missing ones */
     95:     for (s = up - 1; s >= low; s--) {
     96:       PetscReal t = tshhist[s];
     97:       PetscInt  tid;

     99:       if (tj->lag.TT[tj->lag.order + 1 + s - low]) continue;
    100:       tid = LagrangeGetId(t, tj->lag.order + 1, tj->lag.T, tj->lag.TT);
    101:       PetscCheck(tid < 0, PetscObjectComm((PetscObject)tj), PETSC_ERR_PLIB, "This should not happen");
    102:       tid = -tid - 1;
    103:       if (tj->monitor) {
    104:         if (tj->lag.T[tid] < PETSC_MAX_REAL) {
    105:           PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Discarding snapshot %" PetscInt_FMT " at time %g\n", tid, (double)tj->lag.T[tid]));
    106:         } else {
    107:           PetscCall(PetscViewerASCIIPrintf(tj->monitor, "New snapshot %" PetscInt_FMT "\n", tid));
    108:         }
    109:         PetscCall(PetscViewerASCIIPushTab(tj->monitor));
    110:       }
    111:       PetscCall(TSTrajectoryGetVecs(tj, ts, tshhist_id[s], &t, tj->lag.W[tid], NULL));
    112:       tj->lag.T[tid] = t;
    113:       if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor));
    114:       tj->lag.TT[tid]                         = PETSC_TRUE;
    115:       tj->lag.WW[cnt]                         = tj->lag.W[tid];
    116:       tj->lag.TW[cnt]                         = t;
    117:       tj->lag.TT[tj->lag.order + 1 + s - low] = PETSC_TRUE;
    118:       cnt++;
    119:     }
    120:     if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor));
    121:   }
    122:   PetscCall(PetscArrayzero(tj->lag.TT, tj->lag.order + 1));
    123:   if (id >= 0 && U) { /* requested time match */
    124:     PetscInt tid = LagrangeGetId(t, tj->lag.order + 1, tj->lag.T, tj->lag.TT);
    125:     if (tj->monitor) {
    126:       PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Retrieving solution from exact step\n"));
    127:       PetscCall(PetscViewerASCIIPushTab(tj->monitor));
    128:     }
    129:     if (tid < 0) {
    130:       tid = -tid - 1;
    131:       if (tj->monitor) {
    132:         if (tj->lag.T[tid] < PETSC_MAX_REAL) {
    133:           PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Discarding snapshot %" PetscInt_FMT " at time %g\n", tid, (double)tj->lag.T[tid]));
    134:         } else {
    135:           PetscCall(PetscViewerASCIIPrintf(tj->monitor, "New snapshot %" PetscInt_FMT "\n", tid));
    136:         }
    137:         PetscCall(PetscViewerASCIIPushTab(tj->monitor));
    138:       }
    139:       PetscCall(TSTrajectoryGetVecs(tj, ts, tshhist_id[id], &t, tj->lag.W[tid], NULL));
    140:       if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor));
    141:       tj->lag.T[tid] = t;
    142:     } else if (tj->monitor) {
    143:       PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Reusing snapshot %" PetscInt_FMT " step %" PetscInt_FMT ", time %g\n", tid, tshhist_id[id], (double)t));
    144:     }
    145:     PetscCall(VecCopy(tj->lag.W[tid], U));
    146:     PetscCall(PetscObjectStateGet((PetscObject)U, &tj->lag.Ucached.state));
    147:     PetscCall(PetscObjectGetId((PetscObject)U, &tj->lag.Ucached.id));
    148:     tj->lag.Ucached.time = t;
    149:     tj->lag.Ucached.step = tshhist_id[id];
    150:     if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor));
    151:   }
    152:   if (id < 0 && U) {
    153:     if (tj->monitor) PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Interpolating solution with %" PetscInt_FMT " snapshots\n", cnt));
    154:     LagrangeBasisVals(cnt, t, tj->lag.TW, tj->lag.L);
    155:     PetscCall(VecZeroEntries(U));
    156:     PetscCall(VecMAXPY(U, cnt, tj->lag.L, tj->lag.WW));
    157:     PetscCall(PetscObjectStateGet((PetscObject)U, &tj->lag.Ucached.state));
    158:     PetscCall(PetscObjectGetId((PetscObject)U, &tj->lag.Ucached.id));
    159:     tj->lag.Ucached.time = t;
    160:     tj->lag.Ucached.step = PETSC_INT_MIN;
    161:   }
    162:   if (Udot) {
    163:     if (tj->monitor) PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Interpolating derivative with %" PetscInt_FMT " snapshots\n", cnt));
    164:     LagrangeBasisDers(cnt, t, tj->lag.TW, tj->lag.L);
    165:     PetscCall(VecZeroEntries(Udot));
    166:     PetscCall(VecMAXPY(Udot, cnt, tj->lag.L, tj->lag.WW));
    167:     PetscCall(PetscObjectStateGet((PetscObject)Udot, &tj->lag.Udotcached.state));
    168:     PetscCall(PetscObjectGetId((PetscObject)Udot, &tj->lag.Udotcached.id));
    169:     tj->lag.Udotcached.time = t;
    170:     tj->lag.Udotcached.step = PETSC_INT_MIN;
    171:   }
    172:   PetscFunctionReturn(PETSC_SUCCESS);
    173: }