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) { maxloc = (max < T[_tid] && !Taken[_tid]) ? (max = T[_tid],_tid) : maxloc; _tid++; } 42 return -maxloc-1; 43 } 44 } 45 46 PetscErrorCode TSTrajectoryReconstruct_Private(TSTrajectory tj,TS ts,PetscReal t,Vec U,Vec Udot) 47 { 48 TSHistory tsh = tj->tsh; 49 const PetscReal *tshhist; 50 const PetscInt *tshhist_id; 51 PetscInt id, cnt, i, tshn; 52 53 PetscFunctionBegin; 54 PetscCall(TSHistoryGetLocFromTime(tsh,t,&id)); 55 PetscCall(TSHistoryGetHistory(tsh,&tshn,&tshhist,&tshhist_id,NULL)); 56 if (id == -1 || id == -tshn - 1) { 57 PetscReal t0 = tshn ? tshhist[0] : 0.0; 58 PetscReal tf = tshn ? tshhist[tshn-1] : 0.0; 59 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); 60 } 61 if (tj->monitor) { 62 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Reconstructing at time %g, order %" PetscInt_FMT "\n",(double)t,tj->lag.order)); 63 } 64 if (!tj->lag.T) { 65 PetscInt o = tj->lag.order+1; 66 PetscCall(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 PetscCall(VecDuplicateVecs(U ? U : Udot,o,&tj->lag.W)); 69 } 70 cnt = 0; 71 PetscCall(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; 74 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) PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 79 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) { 86 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %" PetscInt_FMT ", step %" PetscInt_FMT ", time %g\n",tid,tshhist_id[s],(double)t)); 87 } 88 tj->lag.TT[tid] = PETSC_TRUE; 89 tj->lag.WW[cnt] = tj->lag.W[tid]; 90 tj->lag.TW[cnt] = t; 91 tj->lag.TT[tj->lag.order+1 + s-low] = PETSC_TRUE; /* tell the next loop to skip it */ 92 cnt++; 93 } 94 95 /* now load the missing ones */ 96 for (s = up-1; s >= low; s--) { 97 PetscReal t = tshhist[s]; 98 PetscInt tid; 99 100 if (tj->lag.TT[tj->lag.order+1 + s-low]) continue; 101 tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT); 102 PetscCheck(tid < 0,PetscObjectComm((PetscObject)tj),PETSC_ERR_PLIB,"This should not happen"); 103 tid = -tid-1; 104 if (tj->monitor) { 105 if (tj->lag.T[tid] < PETSC_MAX_REAL) { 106 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %" PetscInt_FMT " at time %g\n",tid,(double)tj->lag.T[tid])); 107 } else { 108 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"New snapshot %" PetscInt_FMT "\n",tid)); 109 } 110 PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 111 } 112 PetscCall(TSTrajectoryGetVecs(tj,ts,tshhist_id[s],&t,tj->lag.W[tid],NULL)); 113 tj->lag.T[tid] = t; 114 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 115 tj->lag.TT[tid] = PETSC_TRUE; 116 tj->lag.WW[cnt] = tj->lag.W[tid]; 117 tj->lag.TW[cnt] = t; 118 tj->lag.TT[tj->lag.order+1 + s-low] = PETSC_TRUE; 119 cnt++; 120 } 121 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 122 } 123 PetscCall(PetscArrayzero(tj->lag.TT,tj->lag.order+1)); 124 if (id >=0 && U) { /* requested time match */ 125 PetscInt tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT); 126 if (tj->monitor) { 127 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Retrieving solution from exact step\n")); 128 PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 129 } 130 if (tid < 0) { 131 tid = -tid-1; 132 if (tj->monitor) { 133 if (tj->lag.T[tid] < PETSC_MAX_REAL) { 134 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %" PetscInt_FMT " at time %g\n",tid,(double)tj->lag.T[tid])); 135 } else { 136 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"New snapshot %" PetscInt_FMT "\n",tid)); 137 } 138 PetscCall(PetscViewerASCIIPushTab(tj->monitor)); 139 } 140 PetscCall(TSTrajectoryGetVecs(tj,ts,tshhist_id[id],&t,tj->lag.W[tid],NULL)); 141 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 142 tj->lag.T[tid] = t; 143 } else if (tj->monitor) { 144 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %" PetscInt_FMT " step %" PetscInt_FMT ", time %g\n",tid,tshhist_id[id],(double)t)); 145 } 146 PetscCall(VecCopy(tj->lag.W[tid],U)); 147 PetscCall(PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state)); 148 PetscCall(PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id)); 149 tj->lag.Ucached.time = t; 150 tj->lag.Ucached.step = tshhist_id[id]; 151 if (tj->monitor) PetscCall(PetscViewerASCIIPopTab(tj->monitor)); 152 } 153 if (id < 0 && U) { 154 if (tj->monitor) { 155 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Interpolating solution with %" PetscInt_FMT " snapshots\n",cnt)); 156 } 157 LagrangeBasisVals(cnt,t,tj->lag.TW,tj->lag.L); 158 PetscCall(VecZeroEntries(U)); 159 PetscCall(VecMAXPY(U,cnt,tj->lag.L,tj->lag.WW)); 160 PetscCall(PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state)); 161 PetscCall(PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id)); 162 tj->lag.Ucached.time = t; 163 tj->lag.Ucached.step = PETSC_MIN_INT; 164 } 165 if (Udot) { 166 if (tj->monitor) { 167 PetscCall(PetscViewerASCIIPrintf(tj->monitor,"Interpolating derivative with %" PetscInt_FMT " snapshots\n",cnt)); 168 } 169 LagrangeBasisDers(cnt,t,tj->lag.TW,tj->lag.L); 170 PetscCall(VecZeroEntries(Udot)); 171 PetscCall(VecMAXPY(Udot,cnt,tj->lag.L,tj->lag.WW)); 172 PetscCall(PetscObjectStateGet((PetscObject)Udot,&tj->lag.Udotcached.state)); 173 PetscCall(PetscObjectGetId((PetscObject)Udot,&tj->lag.Udotcached.id)); 174 tj->lag.Udotcached.time = t; 175 tj->lag.Udotcached.step = PETSC_MIN_INT; 176 } 177 PetscFunctionReturn(0); 178 } 179