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: 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) L[k] *= (t - T[j])/(T[k] - T[j]);
     11:     }
     12:   }
     13: }

     15: PETSC_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: PETSC_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;
     52:   PetscErrorCode  ierr;

     55:   TSHistoryGetLocFromTime(tsh,t,&id);
     56:   TSHistoryGetHistory(tsh,&tshn,&tshhist,&tshhist_id,NULL);
     57:   if (id == -1 || id == -tshn - 1) {
     58:     PetscReal t0 = tshn ? tshhist[0]      : 0.0;
     59:     PetscReal tf = tshn ? tshhist[tshn-1] : 0.0;
     60:     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);
     61:   }
     62:   if (tj->monitor) {
     63:     PetscViewerASCIIPrintf(tj->monitor,"Reconstructing at time %g, order %D\n",(double)t,tj->lag.order);
     64:   }
     65:   if (!tj->lag.T) {
     66:     PetscInt o = tj->lag.order+1;
     67:     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:     VecDuplicateVecs(U ? U : Udot,o,&tj->lag.W);
     70:   }
     71:   cnt = 0;
     72:   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;

     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) {
     80:       PetscViewerASCIIPushTab(tj->monitor);
     81:     }

     83:     /* first see if we can reuse any */
     84:     for (s = up-1; s >= low; s--) {
     85:       PetscReal t = tshhist[s];
     86:       PetscInt tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT);
     87:       if (tid < 0) continue;
     88:       if (tj->monitor) {
     89:         PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %D, step %D, time %g\n",tid,tshhist_id[s],(double)t);
     90:       }
     91:       tj->lag.TT[tid] = PETSC_TRUE;
     92:       tj->lag.WW[cnt] = tj->lag.W[tid];
     93:       tj->lag.TW[cnt] = t;
     94:       tj->lag.TT[tj->lag.order+1 + s-low] = PETSC_TRUE; /* tell the next loop to skip it */
     95:       cnt++;
     96:     }

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

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