This sounds like a problem in your calculation function where a Vec or Mat does not get reset to 0, but it does in your by hand code.
Matt On Mar 20, 2015 2:52 PM, "Andrew Spott" <[email protected]> wrote: > I have a fairly simple problem that I’m trying to timestep: > > u’ = A(t) u > > I’m using the crank-nicholson method, which I understand (for this > problem) to be: > > u(t + h) = u(t) + h/2[A(t+h)*u(t+h) + A(t)*u(t)] > or > [1 - h/2 * A(t+1)] u(t+1) = [1 + h/2 * A(t)] u(t) > > When I attempt to timestep using PETSc, the norm of `u` blows up. When I > do it directly (using the above), the norm of `u` doesn’t blow up. > > It is important to note that the solution generated after the first step > is identical for both, but the second step for Petsc has a norm of ~2, > while for the directly calculated version it is ~1. The third step for > petsc has a norm of ~4, while the directly calculated version it is still > ~1. > > I’m not sure what I’m doing wrong. > > PETSc code is taken out of the manual and is pretty simple: > > TSCreate( comm, &ts ); > TSSetProblemType( ts, TS_LINEAR); > TSSetType( ts, TSCN ); > TSSetInitialTimeStep( ts, 0, 0.01 ); > TSSetDuration( ts, 5, 0.03 ); > TSSetFromOptions( ts ); > TSSetRHSFunction( ts, NULL, TSComputeRHSFunctionLinear, NULL ); > TSSetRHSJacobian( ts, A, A, func, &cntx ); > TSSolve( ts, psi0 ); > > `func` just constructs A(t) at the time given. The same code for > calculating A(t) is used in both calculations, along with the same initial > vector psi0, and the same time steps. > > Let me know what other information is needed. I’m not sure what could be > the problem. `func` doesn’t touch U at all (should it?). > > -Andrew > >
