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

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

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

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

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