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