1 #include <petsc/private/tshistoryimpl.h> 2 #include <petscts.h> 3 4 /* these two functions have been stolen from bdf.c */ 5 PETSC_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) 11 L[k] *= (t - T[j])/(T[k] - T[j]); 12 } 13 14 PETSC_STATIC_INLINE void LagrangeBasisDers(PetscInt n,PetscReal t,const PetscReal T[],PetscScalar dL[]) 15 { 16 PetscInt k,j,i; 17 for (k=0; k<n; k++) 18 for (dL[k]=0, j=0; j<n; j++) 19 if (j != k) { 20 PetscReal L = 1/(T[k] - T[j]); 21 for (i=0; i<n; i++) 22 if (i != j && i != k) 23 L *= (t - T[i])/(T[k] - T[i]); 24 dL[k] += L; 25 } 26 } 27 28 PETSC_STATIC_INLINE PetscInt LagrangeGetId(PetscReal t, PetscInt n, const PetscReal T[], const PetscBool Taken[]) 29 { 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) { maxloc = (max < T[_tid] && !Taken[_tid]) ? (max = T[_tid],_tid) : maxloc; _tid++; } 39 return -maxloc-1; 40 } 41 } 42 43 PetscErrorCode TSTrajectoryReconstruct_Private(TSTrajectory tj,TS ts,PetscReal t,Vec U,Vec Udot) 44 { 45 TSHistory tsh = tj->tsh; 46 const PetscReal *tshhist; 47 const PetscInt *tshhist_id; 48 PetscInt id, cnt, i, tshn; 49 PetscErrorCode ierr; 50 51 PetscFunctionBegin; 52 ierr = TSHistoryGetLocFromTime(tsh,t,&id);CHKERRQ(ierr); 53 ierr = TSHistoryGetHistory(tsh,&tshn,&tshhist,&tshhist_id,NULL);CHKERRQ(ierr); 54 if (id == -1 || id == -tshn - 1) { 55 PetscReal t0 = tshn ? tshhist[0] : 0.0; 56 PetscReal tf = tshn ? tshhist[tshn-1] : 0.0; 57 SETERRQ4(PetscObjectComm((PetscObject)tj),PETSC_ERR_PLIB,"Requested time %g is outside the history interval [%g, %g] (%d)",(double)t,(double)t0,(double)tf,tshn); 58 } 59 if (tj->monitor) { 60 ierr = PetscViewerASCIIPrintf(tj->monitor,"Reconstructing at time %g, order %D\n",(double)t,tj->lag.order);CHKERRQ(ierr); 61 } 62 if (!tj->lag.T) { 63 PetscInt o = tj->lag.order+1; 64 ierr = PetscMalloc5(o,&tj->lag.L,o,&tj->lag.T,o,&tj->lag.WW,2*o,&tj->lag.TT,o,&tj->lag.TW);CHKERRQ(ierr); 65 for (i = 0; i < o; i++) tj->lag.T[i] = PETSC_MAX_REAL; 66 ierr = VecDuplicateVecs(U ? U : Udot,o,&tj->lag.W);CHKERRQ(ierr); 67 } 68 cnt = 0; 69 ierr = PetscMemzero(tj->lag.TT,2*(tj->lag.order+1)*sizeof(PetscBool));CHKERRQ(ierr); 70 if (id < 0 || Udot) { /* populate snapshots for interpolation */ 71 PetscInt s,nid = id < 0 ? -(id+1) : id; 72 73 PetscInt up = PetscMin(nid + tj->lag.order/2+1,tshn); 74 PetscInt low = PetscMax(up-tj->lag.order-1,0); 75 up = PetscMin(PetscMax(low + tj->lag.order + 1,up),tshn); 76 if (tj->monitor) { 77 ierr = PetscViewerASCIIPushTab(tj->monitor);CHKERRQ(ierr); 78 } 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 ierr = PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %D, step %D, time %g\n",tid,tshhist_id[s],(double)t);CHKERRQ(ierr); 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 if (tj->lag.TT[tj->lag.order+1 + s-low]) continue; 99 PetscInt tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT); 100 if (tid >= 0) SETERRQ(PetscObjectComm((PetscObject)tj),PETSC_ERR_PLIB,"This should not happen"); 101 tid = -tid-1; 102 if (tj->monitor) { 103 if (tj->lag.T[tid] < PETSC_MAX_REAL) { 104 ierr = PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %D at time %g\n",tid,(double)tj->lag.T[tid]);CHKERRQ(ierr); 105 } else { 106 ierr = PetscViewerASCIIPrintf(tj->monitor,"New snapshot %D\n",tid);CHKERRQ(ierr); 107 } 108 ierr = PetscViewerASCIIPushTab(tj->monitor);CHKERRQ(ierr); 109 } 110 ierr = TSTrajectoryGetVecs(tj,ts,tshhist_id[s],&t,tj->lag.W[tid],NULL);CHKERRQ(ierr); 111 tj->lag.T[tid] = t; 112 if (tj->monitor) { 113 ierr = PetscViewerASCIIPopTab(tj->monitor);CHKERRQ(ierr); 114 } 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) { 122 ierr = PetscViewerASCIIPopTab(tj->monitor);CHKERRQ(ierr); 123 } 124 } 125 ierr = PetscMemzero(tj->lag.TT,(tj->lag.order+1)*sizeof(PetscBool));CHKERRQ(ierr); 126 if (id >=0 && U) { /* requested time match */ 127 PetscInt tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT); 128 if (tj->monitor) { 129 ierr = PetscViewerASCIIPrintf(tj->monitor,"Retrieving solution from exact step\n");CHKERRQ(ierr); 130 ierr = PetscViewerASCIIPushTab(tj->monitor);CHKERRQ(ierr); 131 } 132 if (tid < 0) { 133 tid = -tid-1; 134 if (tj->monitor) { 135 if (tj->lag.T[tid] < PETSC_MAX_REAL) { 136 ierr = PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %D at time %g\n",tid,(double)tj->lag.T[tid]);CHKERRQ(ierr); 137 } else { 138 ierr = PetscViewerASCIIPrintf(tj->monitor,"New snapshot %D\n",tid);CHKERRQ(ierr); 139 } 140 ierr = PetscViewerASCIIPushTab(tj->monitor);CHKERRQ(ierr); 141 } 142 ierr = TSTrajectoryGetVecs(tj,ts,tshhist_id[id],&t,tj->lag.W[tid],NULL);CHKERRQ(ierr); 143 if (tj->monitor) { 144 ierr = PetscViewerASCIIPopTab(tj->monitor);CHKERRQ(ierr); 145 } 146 tj->lag.T[tid] = t; 147 } else if (tj->monitor) { 148 ierr = PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %D step %D, time %g\n",tid,tshhist_id[id],(double)t);CHKERRQ(ierr); 149 } 150 ierr = VecCopy(tj->lag.W[tid],U);CHKERRQ(ierr); 151 ierr = PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);CHKERRQ(ierr); 152 ierr = PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);CHKERRQ(ierr); 153 tj->lag.Ucached.time = t; 154 tj->lag.Ucached.step = tshhist_id[id]; 155 if (tj->monitor) { 156 ierr = PetscViewerASCIIPopTab(tj->monitor);CHKERRQ(ierr); 157 } 158 } 159 if (id < 0 && U) { 160 if (tj->monitor) { 161 ierr = PetscViewerASCIIPrintf(tj->monitor,"Interpolating solution with %D snapshots\n",cnt);CHKERRQ(ierr); 162 } 163 LagrangeBasisVals(cnt,t,tj->lag.TW,tj->lag.L); 164 ierr = VecZeroEntries(U);CHKERRQ(ierr); 165 ierr = VecMAXPY(U,cnt,tj->lag.L,tj->lag.WW);CHKERRQ(ierr); 166 ierr = PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);CHKERRQ(ierr); 167 ierr = PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);CHKERRQ(ierr); 168 tj->lag.Ucached.time = t; 169 tj->lag.Ucached.step = PETSC_MIN_INT; 170 } 171 if (Udot) { 172 if (tj->monitor) { 173 ierr = PetscViewerASCIIPrintf(tj->monitor,"Interpolating derivative with %D snapshots\n",cnt);CHKERRQ(ierr); 174 } 175 LagrangeBasisDers(cnt,t,tj->lag.TW,tj->lag.L); 176 ierr = VecZeroEntries(Udot);CHKERRQ(ierr); 177 ierr = VecMAXPY(Udot,cnt,tj->lag.L,tj->lag.WW);CHKERRQ(ierr); 178 ierr = PetscObjectStateGet((PetscObject)Udot,&tj->lag.Udotcached.state);CHKERRQ(ierr); 179 ierr = PetscObjectGetId((PetscObject)Udot,&tj->lag.Udotcached.id);CHKERRQ(ierr); 180 tj->lag.Udotcached.time = t; 181 tj->lag.Udotcached.step = PETSC_MIN_INT; 182 } 183 PetscFunctionReturn(0); 184 } 185