<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Mar 20, 2015 at 3:09 PM, Andrew Spott <span dir="ltr"><<a href="mailto:ansp6066@colorado.edu" target="_blank">ansp6066@colorado.edu</a>></span> 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>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></span></div></blockquote><div><br></div><div>I would think you would zero the residual vector (if you add to it to construct the residual, as in FEM methods), not the solution.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><span>
<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></span></div></blockquote><div><br></div><div>Are you initializing the Jacobian to 0 first?</div><div><br></div><div>  Thanks,</div><div><br></div><div>     Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><span>
<div>Is this the case? or is there some other responsibility of said function?</div><span class="HOEnZb"><font color="#888888">
<div><br></div>
<div>-Andrew</div></font></span><span class="">
<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><div><div class="h5"><span 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></div></div></span><div><div class="h5"><span><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" target="_blank">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>
</div></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>