<span id="mailbox-conversation"><div>Sorry it took so long, I wanted to create a “reduced” case (without all my parameter handling and other stuff…)</div>

<div><br></div>

<div>

<div id="mb-reply">https://gist.github.com/spott/aea8070f35e79e7249e6</div>

<div id="mb-reply"><br></div>

<div id="mb-reply">The first section does it using the time stepper.  The second section does it by explicitly doing the steps.  The output is:</div>

<div id="mb-reply"><br></div>

<div id="mb-reply">//first section, using TimeStepper:</div>

<div id="mb-reply">

<div>t: 0 step: 0 norm-1: 0</div>

<div>t: 0.01 step: 1 norm-1: 0</div>

<div>t: 0.02 step: 2 norm-1: 0.999995</div>

<div id="mb-reply">t: 0.03 step: 3 norm-1: 2.99998</div>

<div id="mb-reply"><br></div>

<div id="mb-reply">//Second section, using explicit code.</div>

<div>t: 0.01 norm-1: 0</div>

<div>t: 0.02 norm-1: 0</div>

<div>t: 0.02 norm-1: 2.22045e-16</div>

</div>

</div></span><div class="mailbox_signature"><br></div>

<br><br><div class="gmail_quote"><p>On Fri, Mar 20, 2015 at 4:45 PM, Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span> wrote:<br></p><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;"><p>Andrew,

Send your entire code. It will be easier and faster than talking past each other.

<br><br>  Barry

<br><br>> On Mar 20, 2015, at 5:00 PM, Andrew Spott <ansp6066@colorado.edu> wrote:

<br>> 

> I’m sorry, I’m not trying to be difficult, but I’m not following.

<br>> 

<br>> The manual states (for my special case):

<br>>  • u ̇ = A(t)u. Use

<br>> 

<br>> TSSetProblemType(ts,TS LINEAR); TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL); TSSetRHSJacobian(ts,A,A,YourComputeRHSJacobian,&appctx);

<br>> 

> where YourComputeRHSJacobian() is a function you provide that computes A as a function of time. Or use ...

<br>> My `func` does this.  It is 7 lines:

<br>> 

<br>>     context* c = static_cast<context*>( G_u );

<br>>     PetscScalar e = c->E( t_ );

<br>>     MatCopy( c->D, A, SAME_NONZERO_PATTERN );

<br>>     MatShift( A, e );

<br>>     MatDiagonalSet( A, c->H0, INSERT_VALUES);

<br>>     MatShift( A, std::complex<double>( 0, -1 ) );

<br>>     return 0;

<br>> 

> SHOULD `func` touch U?  If so, what should `func` do to U?  I thought that the RHSJacobian function was only meant to create A, since dG/du = A(t) (for this special case).

<br>> 

<br>> -Andrew

<br>> 

<br>> 

<br>> 

<br>> On Fri, Mar 20, 2015 at 3:26 PM, Matthew Knepley <knepley@gmail.com> wrote:

<br>> 

<br>> On Fri, Mar 20, 2015 at 3:09 PM, Andrew Spott <ansp6066@colorado.edu> wrote:

> 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.

<br>> 

> 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.

<br>>  

> 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).

<br>> 

> Are you initializing the Jacobian to 0 first?

<br>> 

<br>>   Thanks,

<br>> 

<br>>      Matt

<br>>  

> Is this the case? or is there some other responsibility of said function?

<br>> 

<br>> -Andrew

<br>> 

<br>> >Ah ha!

<br>> >

> >The function passed to TSSetRHSJacobian needs to zero the solution vector?

<br>> >

> >As a point, this isn’t mentioned in any documentation that I can find.

<br>> >

<br>> >-Andrew

<br>> 

<br>> On Friday, Mar 20, 2015 at 2:17 PM, Matthew Knepley <knepley@gmail.com>, 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.

<br>> 

<br>>    Matt

<br>> 

<br>> On Mar 20, 2015 2:52 PM, "Andrew Spott" <ansp6066@colorado.edu> wrote:

<br>> I have a fairly simple problem that I’m trying to timestep:

<br>> 

<br>> u’ = A(t) u

<br>> 

<br>> I’m using the crank-nicholson method, which I understand (for this problem) to be:

<br>> 

<br>> u(t + h) = u(t) + h/2[A(t+h)*u(t+h) + A(t)*u(t)]

<br>> or

<br>> [1 - h/2 * A(t+1)] u(t+1) = [1 + h/2 * A(t)] u(t)

<br>> 

> 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.

<br>> 

> 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.

<br>> 

> I’m not sure what I’m doing wrong.

<br>> 

<br>> PETSc code is taken out of the manual and is pretty simple:

<br>> 

<br>>         TSCreate( comm, &ts );

<br>>         TSSetProblemType( ts, TS_LINEAR);

<br>>         TSSetType( ts, TSCN );

<br>>         TSSetInitialTimeStep( ts, 0, 0.01 );

<br>>         TSSetDuration( ts, 5, 0.03 );

<br>>         TSSetFromOptions( ts );

<br>>         TSSetRHSFunction( ts, NULL, TSComputeRHSFunctionLinear, NULL );

<br>>         TSSetRHSJacobian( ts, A, A, func, &cntx );

<br>>         TSSolve( ts, psi0 );

<br>> 

> `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.

<br>> 

> 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?).

<br>> 

<br>> -Andrew

<br>> 

<br>> 

<br>> 

<br>> 

<br>> -- 

> 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

<br>> 

<br><br></p></blockquote></div><br>