So, it doesn’t seem that zeroing the given vector in the function passed to TSSetRHSJacobian is the problem. When I do that, it just zeros out the solution.
The function that is passed to TSSetRHSJacobian has only one responsibility — to create the jacobian — correct? In my case this is A(t). The solution vector is given for when you are solving nonlinear problems (A(t) also depends on U(t)). In my case, I don’t even look at the solution vector (because my A(t) doesn’t depend on it). Is this the case? or is there some other responsibility of said function? -Andrew >Ah ha! > >The function passed to TSSetRHSJacobian needs to zero the solution vector? > >As a point, this isn’t mentioned in any documentation that I can find. > >-Andrew On Friday, Mar 20, 2015 at 2:17 PM, Matthew Knepley <[email protected]>, wrote: 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
