1 #include <petsc/private/tshistoryimpl.h>
2 #include <petscts.h>
3
4 /* these two functions have been stolen from bdf.c */
LagrangeBasisVals(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar L[])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 }
14
LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[])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 }
30
LagrangeGetId(PetscReal t,PetscInt n,const PetscReal T[],const PetscBool Taken[])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 }
48
TSTrajectoryReconstruct_Private(TSTrajectory tj,TS ts,PetscReal t,Vec U,Vec Udot)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;
55
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;
75
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));
80
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 }
93
94 /* now load the missing ones */
95 for (s = up - 1; s >= low; s--) {
96 PetscReal t = tshhist[s];
97 PetscInt tid;
98
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 }
174