<p dir="ltr">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.</p>
<p dir="ltr">   Matt</p>
<div class="gmail_quote">On Mar 20, 2015 2:52 PM, "Andrew Spott" <<a href="mailto:ansp6066@colorado.edu">ansp6066@colorado.edu</a>> wrote:<br type="attribution"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><u></u>
<div>
<span><div>I have a fairly simple problem that I’m trying to timestep:</div>
<div><br></div>
<div>u’ = A(t) u</div>
<div><br></div>
<div>I’m using the crank-nicholson method, which I understand (for this problem) to be:</div>
<div><br></div>
<div>u(t + h) = u(t) + h/2[A(t+h)*u(t+h) + A(t)*u(t)]</div>
<div>or</div>
<div>[1 - h/2 * A(t+1)] u(t+1) = [1 + h/2 * A(t)] u(t)</div>
<div><br></div>
<div>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.</div>
<div><br></div>
<div>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.</div>
<div><br></div>
<div>I’m not sure what I’m doing wrong.</div>
<div><br></div>
<div>PETSc code is taken out of the manual and is pretty simple:</div>
<div><br></div>
<div>
<div>        TSCreate( comm, &ts );</div>
<div>        TSSetProblemType( ts, TS_LINEAR);</div>
<div>        TSSetType( ts, TSCN );</div>
<div>        TSSetInitialTimeStep( ts, 0, 0.01 );</div>
<div>        TSSetDuration( ts, 5, 0.03 );</div>
<div>        TSSetFromOptions( ts );</div>
<div>        TSSetRHSFunction( ts, NULL, TSComputeRHSFunctionLinear, NULL );</div>
<div>        TSSetRHSJacobian( ts, A, A, func, &cntx );</div>
<div>        TSSolve( ts, psi0 );</div>
<div><br></div>
<div>`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.</div>
<div><br></div>
<div>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?).</div>
<div><br></div>
<div>-Andrew</div>
</div></span><div><br></div>
</div>
</blockquote></div>