1 #include <petsc/private/tshistoryimpl.h> 2 #include <petscts.h> 3 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 PetscInt k, j; 7 for (k = 0; k < n; k++) { 8 for (L[k] = 1, j = 0; j < n; j++) { 9 if (j != k) L[k] *= (t - T[j]) / (T[k] - T[j]); 10 } 11 } 12 } 13 14 static inline void LagrangeBasisDers(PetscInt n, PetscReal t, const PetscReal T[], PetscScalar dL[]) { 15 PetscInt k, j, i; 16 for (k = 0; k < n; k++) { 17 for (dL[k] = 0, j = 0; j < n; j++) { 18 if (j != k) { 19 PetscReal L = 1 / (T[k] - T[j]); 20 for (i = 0; i < n; i++) { 21 if (i != j && i != k) L *= (t - T[i]) / (T[k] - T[i]); 22 } 23 dL[k] += L; 24 } 25 } 26 } 27 } 28 29 static inline PetscInt LagrangeGetId(PetscReal t, PetscInt n, const PetscReal T[], const PetscBool Taken[]) { 30 PetscInt _tid = 0; 31 while (_tid < n && PetscAbsReal(t - T[_tid]) > PETSC_SMALL) _tid++; 32 if (_tid < n && !Taken[_tid]) { 33 return _tid; 34 } else { /* we get back a negative id, where the maximum time is stored, since we use usually reconstruct backward in time */ 35 PetscReal max = PETSC_MIN_REAL; 36 PetscInt maxloc = n; 37 _tid = 0; 38 while (_tid < n) { 39 maxloc = (max < T[_tid] && !Taken[_tid]) ? (max = T[_tid], _tid) : maxloc; 40 _tid++; 41 } 42 return -maxloc - 1; 43 } 44 } 45 46 PetscErrorCode TSTrajectoryReconstruct_Private(TSTrajectory tj, TS ts, PetscReal t, Vec U, Vec Udot) { 47 TSHistory tsh = tj->tsh; 48 const PetscReal *tshhist; 49 const PetscInt *tshhist_id; 50 PetscInt id, cnt, i, tshn; 51 52 PetscFunctionBegin; 53 PetscCall(TSHistoryGetLocFromTime(tsh, t, &id)); 54 PetscCall(TSHistoryGetHistory(tsh, &tshn, &tshhist, &tshhist_id, NULL)); 55 if (id == -1 || id == -tshn - 1) { 56 PetscReal t0 = tshn ? tshhist[0] : 0.0; 57 PetscReal tf = tshn ? tshhist[tshn - 1] : 0.0; 58 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); 59 } 60 if (tj->monitor) { PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Reconstructing at time %g, order %" PetscInt_FMT "\n", (double)t, tj->lag.order)); } 61 if (!tj->lag.T) { 62 PetscInt o = tj->lag.order + 1; 63 PetscCall(PetscMalloc5(o, &tj->lag.L, o, &tj->lag.T, o, &tj->lag.WW, 2 * o, &tj->lag.TT, o, &tj->lag.TW)); 64 for (i = 0; i < o; i++) tj->lag.T[i] = PETSC_MAX_REAL; 65 PetscCall(VecDuplicateVecs(U ? U : Udot, o, &tj->lag.W)); 66 } 67 cnt = 0; 68 PetscCall(PetscArrayzero(tj->lag.TT, 2 * (tj->lag.order + 1))); 69 if (id < 0 || Udot) { /* populate snapshots for interpolation */ 70 PetscInt s, nid = id < 0 ? -(id + 1) : id; 71 72 PetscInt up = PetscMin(nid + tj->lag.order / 2 + 1, tshn); 73 PetscInt low = PetscMax(up - tj->lag.order - 1, 0); 74 up = PetscMin(PetscMax(low + tj->lag.order + 1, up), tshn); 75 if (tj->monitor) PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 76 77 /* first see if we can reuse any */ 78 for (s = up - 1; s >= low; s--) { 79 PetscReal t = tshhist[s]; 80 PetscInt tid = LagrangeGetId(t, tj->lag.order + 1, tj->lag.T, tj->lag.TT); 81 if (tid < 0) continue; 82 if (tj->monitor) { PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Reusing snapshot %" PetscInt_FMT ", step %" PetscInt_FMT ", time %g\n", tid, tshhist_id[s], (double)t)); } 83 tj->lag.TT[tid] = PETSC_TRUE; 84 tj->lag.WW[cnt] = tj->lag.W[tid]; 85 tj->lag.TW[cnt] = t; 86 tj->lag.TT[tj->lag.order + 1 + s - low] = PETSC_TRUE; /* tell the next loop to skip it */ 87 cnt++; 88 } 89 90 /* now load the missing ones */ 91 for (s = up - 1; s >= low; s--) { 92 PetscReal t = tshhist[s]; 93 PetscInt tid; 94 95 if (tj->lag.TT[tj->lag.order + 1 + s - low]) continue; 96 tid = LagrangeGetId(t, tj->lag.order + 1, tj->lag.T, tj->lag.TT); 97 PetscCheck(tid < 0, PetscObjectComm((PetscObject)tj), PETSC_ERR_PLIB, "This should not happen"); 98 tid = -tid - 1; 99 if (tj->monitor) { 100 if (tj->lag.T[tid] < PETSC_MAX_REAL) { 101 PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Discarding snapshot %" PetscInt_FMT " at time %g\n", tid, (double)tj->lag.T[tid])); 102 } else { 103 PetscCall(PetscViewerASCIIPrintf(tj->monitor, "New snapshot %" PetscInt_FMT "\n", tid)); 104 } 105 PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 106 } 107 PetscCall(TSTrajectoryGetVecs(tj, ts, tshhist_id[s], &t, tj->lag.W[tid], NULL)); 108 tj->lag.T[tid] = t; 109 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 110 tj->lag.TT[tid] = PETSC_TRUE; 111 tj->lag.WW[cnt] = tj->lag.W[tid]; 112 tj->lag.TW[cnt] = t; 113 tj->lag.TT[tj->lag.order + 1 + s - low] = PETSC_TRUE; 114 cnt++; 115 } 116 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 117 } 118 PetscCall(PetscArrayzero(tj->lag.TT, tj->lag.order + 1)); 119 if (id >= 0 && U) { /* requested time match */ 120 PetscInt tid = LagrangeGetId(t, tj->lag.order + 1, tj->lag.T, tj->lag.TT); 121 if (tj->monitor) { 122 PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Retrieving solution from exact step\n")); 123 PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 124 } 125 if (tid < 0) { 126 tid = -tid - 1; 127 if (tj->monitor) { 128 if (tj->lag.T[tid] < PETSC_MAX_REAL) { 129 PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Discarding snapshot %" PetscInt_FMT " at time %g\n", tid, (double)tj->lag.T[tid])); 130 } else { 131 PetscCall(PetscViewerASCIIPrintf(tj->monitor, "New snapshot %" PetscInt_FMT "\n", tid)); 132 } 133 PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 134 } 135 PetscCall(TSTrajectoryGetVecs(tj, ts, tshhist_id[id], &t, tj->lag.W[tid], NULL)); 136 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 137 tj->lag.T[tid] = t; 138 } else if (tj->monitor) { 139 PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Reusing snapshot %" PetscInt_FMT " step %" PetscInt_FMT ", time %g\n", tid, tshhist_id[id], (double)t)); 140 } 141 PetscCall(VecCopy(tj->lag.W[tid], U)); 142 PetscCall(PetscObjectStateGet((PetscObject)U, &tj->lag.Ucached.state)); 143 PetscCall(PetscObjectGetId((PetscObject)U, &tj->lag.Ucached.id)); 144 tj->lag.Ucached.time = t; 145 tj->lag.Ucached.step = tshhist_id[id]; 146 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 147 } 148 if (id < 0 && U) { 149 if (tj->monitor) { PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Interpolating solution with %" PetscInt_FMT " snapshots\n", cnt)); } 150 LagrangeBasisVals(cnt, t, tj->lag.TW, tj->lag.L); 151 PetscCall(VecZeroEntries(U)); 152 PetscCall(VecMAXPY(U, cnt, tj->lag.L, tj->lag.WW)); 153 PetscCall(PetscObjectStateGet((PetscObject)U, &tj->lag.Ucached.state)); 154 PetscCall(PetscObjectGetId((PetscObject)U, &tj->lag.Ucached.id)); 155 tj->lag.Ucached.time = t; 156 tj->lag.Ucached.step = PETSC_MIN_INT; 157 } 158 if (Udot) { 159 if (tj->monitor) { PetscCall(PetscViewerASCIIPrintf(tj->monitor, "Interpolating derivative with %" PetscInt_FMT " snapshots\n", cnt)); } 160 LagrangeBasisDers(cnt, t, tj->lag.TW, tj->lag.L); 161 PetscCall(VecZeroEntries(Udot)); 162 PetscCall(VecMAXPY(Udot, cnt, tj->lag.L, tj->lag.WW)); 163 PetscCall(PetscObjectStateGet((PetscObject)Udot, &tj->lag.Udotcached.state)); 164 PetscCall(PetscObjectGetId((PetscObject)Udot, &tj->lag.Udotcached.id)); 165 tj->lag.Udotcached.time = t; 166 tj->lag.Udotcached.step = PETSC_MIN_INT; 167 } 168 PetscFunctionReturn(0); 169 } 170