Looking at the code for TSInterpolate_Theta (X is the output vector):
1 ierr = VecCopy(ts->vec_sol,th->X);CHKERRQ(ierr); 2 if (th->endpoint) alpha *= th->Theta; 3 ierr = VecWAXPY(X,alpha,th->Xdot,th->X);CHKERRQ(ierr); I can’t figure out why line 1 exists. And, since this appears to just be using a simple taylor series: U(t + dt) ~= U(t) + Udot(t) * dt, I’m not sure why th->endpoint has any effect on this. Is th->Xdot scaled by th->Theta at some point? Attempting to see if something I didn’t understand was going on, I also took a look at TSInterpolate_Euler (X is still the output vector). PetscReal alpha = (ts->ptime - t)/ts->time_step; ierr = VecAXPBY(ts->vec_sol,1.0-alpha,alpha,X);CHKERRQ(ierr); Which doesn’t even modify X (and does modify the solution…). I assume that this is an attempt to find U(t-δt) ~= U(t) (1-δt/Δt) + U(t-Δt)(δt/Δt). However, I’m not sure where U(t-Δt) is supposed to come from… TSInterpolate_Theta is easy enough to fix, but TSInterpolate_Euler seems to need a retained solution, which isn’t supplied. If a retained solution WERE supplied, then the TSInterpolate_Euler version could work for every time step method, and should be a default. I’ll fix these if I’m right, but since I’m still learning the internals of Petsc, I wanted to check and make sure these are problems first. -Andrew
