Actual source code: reconstruct.c
petsc-3.14.6 2021-03-30
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: }