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 { 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 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 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 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_MIN_INT; 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_MIN_INT; 171 } 172 PetscFunctionReturn(PETSC_SUCCESS); 173 } 174