Actual source code: reconstruct.c
petsc-3.11.4 2019-09-28
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)
11: L[k] *= (t - T[j])/(T[k] - T[j]);
12: }
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: }
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: }
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;
52: TSHistoryGetLocFromTime(tsh,t,&id);
53: TSHistoryGetHistory(tsh,&tshn,&tshhist,&tshhist_id,NULL);
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: PetscViewerASCIIPrintf(tj->monitor,"Reconstructing at time %g, order %D\n",(double)t,tj->lag.order);
61: }
62: if (!tj->lag.T) {
63: PetscInt o = tj->lag.order+1;
64: PetscMalloc5(o,&tj->lag.L,o,&tj->lag.T,o,&tj->lag.WW,2*o,&tj->lag.TT,o,&tj->lag.TW);
65: for (i = 0; i < o; i++) tj->lag.T[i] = PETSC_MAX_REAL;
66: VecDuplicateVecs(U ? U : Udot,o,&tj->lag.W);
67: }
68: cnt = 0;
69: PetscMemzero(tj->lag.TT,2*(tj->lag.order+1)*sizeof(PetscBool));
70: if (id < 0 || Udot) { /* populate snapshots for interpolation */
71: PetscInt s,nid = id < 0 ? -(id+1) : id;
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: PetscViewerASCIIPushTab(tj->monitor);
78: }
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: PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %D, step %D, 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: }
95: /* now load the missing ones */
96: for (s = up-1; s >= low; s--) {
97: PetscReal t = tshhist[s];
98: PetscInt tid;
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: if (tid >= 0) SETERRQ(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: PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %D at time %g\n",tid,(double)tj->lag.T[tid]);
107: } else {
108: PetscViewerASCIIPrintf(tj->monitor,"New snapshot %D\n",tid);
109: }
110: PetscViewerASCIIPushTab(tj->monitor);
111: }
112: TSTrajectoryGetVecs(tj,ts,tshhist_id[s],&t,tj->lag.W[tid],NULL);
113: tj->lag.T[tid] = t;
114: if (tj->monitor) {
115: PetscViewerASCIIPopTab(tj->monitor);
116: }
117: tj->lag.TT[tid] = PETSC_TRUE;
118: tj->lag.WW[cnt] = tj->lag.W[tid];
119: tj->lag.TW[cnt] = t;
120: tj->lag.TT[tj->lag.order+1 + s-low] = PETSC_TRUE;
121: cnt++;
122: }
123: if (tj->monitor) {
124: PetscViewerASCIIPopTab(tj->monitor);
125: }
126: }
127: PetscMemzero(tj->lag.TT,(tj->lag.order+1)*sizeof(PetscBool));
128: if (id >=0 && U) { /* requested time match */
129: PetscInt tid = LagrangeGetId(t,tj->lag.order+1,tj->lag.T,tj->lag.TT);
130: if (tj->monitor) {
131: PetscViewerASCIIPrintf(tj->monitor,"Retrieving solution from exact step\n");
132: PetscViewerASCIIPushTab(tj->monitor);
133: }
134: if (tid < 0) {
135: tid = -tid-1;
136: if (tj->monitor) {
137: if (tj->lag.T[tid] < PETSC_MAX_REAL) {
138: PetscViewerASCIIPrintf(tj->monitor,"Discarding snapshot %D at time %g\n",tid,(double)tj->lag.T[tid]);
139: } else {
140: PetscViewerASCIIPrintf(tj->monitor,"New snapshot %D\n",tid);
141: }
142: PetscViewerASCIIPushTab(tj->monitor);
143: }
144: TSTrajectoryGetVecs(tj,ts,tshhist_id[id],&t,tj->lag.W[tid],NULL);
145: if (tj->monitor) {
146: PetscViewerASCIIPopTab(tj->monitor);
147: }
148: tj->lag.T[tid] = t;
149: } else if (tj->monitor) {
150: PetscViewerASCIIPrintf(tj->monitor,"Reusing snapshot %D step %D, time %g\n",tid,tshhist_id[id],(double)t);
151: }
152: VecCopy(tj->lag.W[tid],U);
153: PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);
154: PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);
155: tj->lag.Ucached.time = t;
156: tj->lag.Ucached.step = tshhist_id[id];
157: if (tj->monitor) {
158: PetscViewerASCIIPopTab(tj->monitor);
159: }
160: }
161: if (id < 0 && U) {
162: if (tj->monitor) {
163: PetscViewerASCIIPrintf(tj->monitor,"Interpolating solution with %D snapshots\n",cnt);
164: }
165: LagrangeBasisVals(cnt,t,tj->lag.TW,tj->lag.L);
166: VecZeroEntries(U);
167: VecMAXPY(U,cnt,tj->lag.L,tj->lag.WW);
168: PetscObjectStateGet((PetscObject)U,&tj->lag.Ucached.state);
169: PetscObjectGetId((PetscObject)U,&tj->lag.Ucached.id);
170: tj->lag.Ucached.time = t;
171: tj->lag.Ucached.step = PETSC_MIN_INT;
172: }
173: if (Udot) {
174: if (tj->monitor) {
175: PetscViewerASCIIPrintf(tj->monitor,"Interpolating derivative with %D snapshots\n",cnt);
176: }
177: LagrangeBasisDers(cnt,t,tj->lag.TW,tj->lag.L);
178: VecZeroEntries(Udot);
179: VecMAXPY(Udot,cnt,tj->lag.L,tj->lag.WW);
180: PetscObjectStateGet((PetscObject)Udot,&tj->lag.Udotcached.state);
181: PetscObjectGetId((PetscObject)Udot,&tj->lag.Udotcached.id);
182: tj->lag.Udotcached.time = t;
183: tj->lag.Udotcached.step = PETSC_MIN_INT;
184: }
185: return(0);
186: }