[petsc-users] TimeStepper norm problems.

[Andrew Spott]

Sorry it took so long, I wanted to create a “reduced” case (without all my parameter handling and other stuff…)





https://gist.github.com/spott/aea8070f35e79e7249e6




The first section does it using the time stepper.  The second section does it by explicitly doing the steps.  The output is:




//first section, using TimeStepper:


t: 0 step: 0 norm-1: 0

t: 0.01 step: 1 norm-1: 0

t: 0.02 step: 2 norm-1: 0.999995

t: 0.03 step: 3 norm-1: 2.99998




//Second section, using explicit code.

t: 0.01 norm-1: 0

t: 0.02 norm-1: 0

t: 0.02 norm-1: 2.22045e-16

On Fri, Mar 20, 2015 at 4:45 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:

> Andrew,
>   Send your entire code. It will be easier and faster than talking past each other.
>   Barry
>> On Mar 20, 2015, at 5:00 PM, Andrew Spott <ansp6066 at colorado.edu> wrote:
>> 
>> I’m sorry, I’m not trying to be difficult, but I’m not following.
>> 
>> The manual states (for my special case):
>> 	• u ̇ = A(t)u. Use
>> 
>> TSSetProblemType(ts,TS LINEAR); TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL); TSSetRHSJacobian(ts,A,A,YourComputeRHSJacobian,&appctx);
>> 
>> where YourComputeRHSJacobian() is a function you provide that computes A as a function of time. Or use ...
>> My `func` does this.  It is 7 lines:
>> 
>>     context* c = static_cast<context*>( G_u );
>>     PetscScalar e = c->E( t_ );
>>     MatCopy( c->D, A, SAME_NONZERO_PATTERN );
>>     MatShift( A, e );
>>     MatDiagonalSet( A, c->H0, INSERT_VALUES);
>>     MatShift( A, std::complex<double>( 0, -1 ) );
>>     return 0;
>> 
>> 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).
>> 
>> -Andrew
>> 
>> 
>> 
>> On Fri, Mar 20, 2015 at 3:26 PM, Matthew Knepley <knepley at gmail.com> wrote:
>> 
>> On Fri, Mar 20, 2015 at 3:09 PM, Andrew Spott <ansp6066 at 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.
>> 
>> 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.
>>  
>> 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).
>> 
>> Are you initializing the Jacobian to 0 first?
>> 
>>   Thanks,
>> 
>>      Matt
>>  
>> 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 <knepley at 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.
>> 
>>    Matt
>> 
>> On Mar 20, 2015 2:52 PM, "Andrew Spott" <ansp6066 at colorado.edu> 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
>> 
>> 
>> 
>> 
>> -- 
>> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
>> -- Norbert Wiener
>> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150320/02a9a356/attachment-0001.html>