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: }