<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN" "http://www.w3.org/TR/REC-html40/loose.dtd">
<html><body>
<span id="mailbox-conversation"><div>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.</div>
<div><br></div>
<div>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).</div>
<div><br></div>
<div>Is this the case? or is there some other responsibility of said function?</div>
<div><br></div>
<div>-Andrew</div>
<div><br></div>
<div>>Ah ha!</div>
<div>></div>
<div>>The function passed to TSSetRHSJacobian needs to zero the solution vector?</div>
<div>></div>
<div>>As a point, this isn’t mentioned in any documentation that I can find.</div>
<div>></div>
<div>>-Andrew</div>
<br><span id="orc-full-body-initial-text" style="display:inline;">On Friday, Mar 20, 2015 at 2:17 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>, wrote:<br></span></span><span class="none"><blockquote class="gmail_quote"><div>
<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><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>
</div></blockquote></span>
</body></html>