Actual source code: reconstruct.c

  1: #include <petsc/private/tshistoryimpl.h>
  2: #include <petscts.h>

  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: }

 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: }

 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: }

 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;

 53:   TSHistoryGetLocFromTime(tsh,t,&id);
 54:   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] (%d)",(double)t,(double)t0,(double)tf,tshn);
 59:   }
 60:   if (tj->monitor) {
 61:     PetscViewerASCIIPrintf(tj->monitor,"Reconstructing at time %g, order %D\n",(double)t,tj->lag.order);
 62:   }
 63:   if (!tj->lag.T) {
 64:     PetscInt o = tj->lag.order+1;
 65:     PetscMalloc5(o,&tj->lag.L,o,&tj->lag.T,o,&tj->lag.WW,2*o,&tj->lag.TT,o,&tj->lag.TW);
 66:     for (i = 0; i < o; i++) tj->lag.T[i] = PETSC_MAX_REAL;
 67:     VecDuplicateVecs(U ? U : Udot,o,&tj->lag.W);
 68:   }
 69:   cnt = 0;
 70:   PetscArrayzero(tj->lag.TT,2*(tj->lag.order+1));
 71:   if (id < 0 || Udot) { /* populate snapshots for interpolation */
 72:     PetscInt s,nid = id < 0 ? -(id+1) : id;

 74:     PetscInt up = PetscMin(nid + tj->lag.order/2+1,tshn);
 75:     PetscInt low = PetscMax(up-tj->lag.order-1,0);
 76:     up = PetscMin(PetscMax(low + tj->lag.order + 1,up),tshn);
 77:     if (tj->monitor) {
 78:       PetscViewerASCIIPushTab(tj->monitor);
 79:     }

 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) {
 87:         PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %D, step %D, time %g\n",tid,tshhist_id[s],(double)t);
 88:       }
 89:       tj->lag.TT[tid] = PETSC_TRUE;
 90:       tj->lag.WW[cnt] = tj->lag.W[tid];
 91:       tj->lag.TW[cnt] = t;
 92:       tj->lag.TT[tj->lag.order+1 + s-low] = PETSC_TRUE; /* tell the next loop to skip it */
 93:       cnt++;
 94:     }

 96:     /* now load the missing ones */
 97:     for (s = up-1; s >= low; s--) {
 98:       PetscReal t = tshhist[s];
 99:       PetscInt tid;

101:       if (tj->lag.TT[tj->lag.order+1 + s-low]) continue;
102:       tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT);
104:       tid = -tid-1;
105:       if (tj->monitor) {
106:         if (tj->lag.T[tid] < PETSC_MAX_REAL) {
107:           PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %D at time %g\n",tid,(double)tj->lag.T[tid]);
108:         } else {
109:           PetscViewerASCIIPrintf(tj->monitor,"New snapshot %D\n",tid);
110:         }
111:         PetscViewerASCIIPushTab(tj->monitor);
112:       }
113:       TSTrajectoryGetVecs(tj,ts,tshhist_id[s],&t,tj->lag.W[tid],NULL);
114:       tj->lag.T[tid] = t;
115:       if (tj->monitor) {
116:         PetscViewerASCIIPopTab(tj->monitor);
117:       }
118:       tj->lag.TT[tid] = PETSC_TRUE;
119:       tj->lag.WW[cnt] = tj->lag.W[tid];
120:       tj->lag.TW[cnt] = t;
121:       tj->lag.TT[tj->lag.order+1 + s-low] = PETSC_TRUE;
122:       cnt++;
123:     }
124:     if (tj->monitor) {
125:       PetscViewerASCIIPopTab(tj->monitor);
126:     }
127:   }
128:   PetscArrayzero(tj->lag.TT,tj->lag.order+1);
129:   if (id >=0 && U) { /* requested time match */
130:     PetscInt tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT);
131:     if (tj->monitor) {
132:       PetscViewerASCIIPrintf(tj->monitor,"Retrieving solution from exact step\n");
133:       PetscViewerASCIIPushTab(tj->monitor);
134:     }
135:     if (tid < 0) {
136:       tid = -tid-1;
137:       if (tj->monitor) {
138:         if (tj->lag.T[tid] < PETSC_MAX_REAL) {
139:           PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %D at time %g\n",tid,(double)tj->lag.T[tid]);
140:         } else {
141:           PetscViewerASCIIPrintf(tj->monitor,"New snapshot %D\n",tid);
142:         }
143:         PetscViewerASCIIPushTab(tj->monitor);
144:       }
145:       TSTrajectoryGetVecs(tj,ts,tshhist_id[id],&t,tj->lag.W[tid],NULL);
146:       if (tj->monitor) {
147:         PetscViewerASCIIPopTab(tj->monitor);
148:       }
149:       tj->lag.T[tid] = t;
150:     } else if (tj->monitor) {
151:       PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %D step %D, time %g\n",tid,tshhist_id[id],(double)t);
152:     }
153:     VecCopy(tj->lag.W[tid],U);
154:     PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);
155:     PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);
156:     tj->lag.Ucached.time = t;
157:     tj->lag.Ucached.step = tshhist_id[id];
158:     if (tj->monitor) {
159:       PetscViewerASCIIPopTab(tj->monitor);
160:     }
161:   }
162:   if (id < 0 && U) {
163:     if (tj->monitor) {
164:       PetscViewerASCIIPrintf(tj->monitor,"Interpolating solution with %D snapshots\n",cnt);
165:     }
166:     LagrangeBasisVals(cnt,t,tj->lag.TW,tj->lag.L);
167:     VecZeroEntries(U);
168:     VecMAXPY(U,cnt,tj->lag.L,tj->lag.WW);
169:     PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);
170:     PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);
171:     tj->lag.Ucached.time = t;
172:     tj->lag.Ucached.step = PETSC_MIN_INT;
173:   }
174:   if (Udot) {
175:     if (tj->monitor) {
176:       PetscViewerASCIIPrintf(tj->monitor,"Interpolating derivative with %D snapshots\n",cnt);
177:     }
178:     LagrangeBasisDers(cnt,t,tj->lag.TW,tj->lag.L);
179:     VecZeroEntries(Udot);
180:     VecMAXPY(Udot,cnt,tj->lag.L,tj->lag.WW);
181:     PetscObjectStateGet((PetscObject)Udot,&tj->lag.Udotcached.state);
182:     PetscObjectGetId((PetscObject)Udot,&tj->lag.Udotcached.id);
183:     tj->lag.Udotcached.time = t;
184:     tj->lag.Udotcached.step = PETSC_MIN_INT;
185:   }
186:   return 0;
187: }