[petsc-users] TimeStepper norm problems. EMIL Please read this

Andrew Spott ansp6066 at colorado.edu
Sat Mar 21 11:26:35 CDT 2015 

>I haven’t been able to compile and run.




You might need -std=c++11 or -std=c++1y.  The code uses a lambda.  If you want, I can refactor that away. 




>The problem appears to be very stiff.





So, one of the things that is weird is that when I explicitly calculate the first few steps using the crank-nicholson method, I don’t have the exponential blowup.  However, when I use the TimeStepper, I DO get an exponential blowup.




Since, theoretically, the two methods should be equivalent and I’m using the same recurrence relationship for the Crank Nicholson method that you are, where does the difference come from?

On Sat, Mar 21, 2015 at 8:32 AM, Emil Constantinescu <emconsta at mcs.anl.gov>
wrote:

> I haven't been able to compile and run. But here are a few quick notes.
> The problem appears to be very stiff.
> Theta and theta_endpoint are defining different methods:
> 1) -ts_type beuler OR -ts_theta_theta 1.0: is backward Euler
> u(t + h) = u(t) + h*A(t+h)*u(t+h)
> 2) -ts_theta_theta 0.5: is the implicit midpoint rule
> u(t + h) = u(t) + h*[A(t+h/2)*(u(t+h)+u(t))/2]
> 3) -ts_type cn OR -ts_theta_theta 0.5 -ts_theta_endpoint: is 
> Crank-Nicholson/trapezoidal
> u(t + h) = u(t) + h/2[A(t+h)*u(t+h) + A(t)*u(t)]
> Note that the last two are different. -ts_type theta -ts_theta_theta .5 
> is different from -ts_type cn. They the same linear stability properties 
> if A(t)=A; but not if A depends on t.
> When -ts_theta_adapt is used, then it detects the instability as an 
> error and reduces the step by a lot! wlte=1.24e+03 which means that the 
> reduction should be severe but the controller tries 0.1*dt and that 
> seems to pass but it "jig-saws" (take a look at the next attempted 
> step), which means that it is likely unstable.
> I'll try to build the example to get more insight.
> Emil
> On 3/20/15 10:57 PM, Barry Smith wrote:
>>
>>    Andrew,
>>
>>     I'm afraid Emil will have to take a look at this and explain it.  The -ts_type beuler and -ts_type theta -ts_theta_theta .5 are stable but the -ts_type cn is not stable. It turns out that -ts_type cn is equivalent to -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint and somehow this endpoint business (which I don't understand) is causing a problem. Meanwhile if I add -ts_theta_adapt to the endpoint one it becomes stable ? Anyways all cases are displayed below.
>>
>>    Emil,
>>
>>      What's up with this? Does the endpoint business have a bug or can it not be used for this problem (the matrix A is a function of t.)
>>
>>    Barry
>>
>>
>> $ ./ex2  -ts_type cn
>> t: 0 step: 0 norm-1: 0
>> t: 0.01 step: 1 norm-1: 0
>> t: 0.02 step: 2 norm-1: 1
>> t: 0.03 step: 3 norm-1: 3
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type theta
>> t: 0 step: 0 norm-1: 0
>> t: 0.01 step: 1 norm-1: 0
>> t: 0.02 step: 2 norm-1: 0
>> t: 0.03 step: 3 norm-1: 0
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type theta -ts_theta_theta .5
>> t: 0 step: 0 norm-1: 0
>> t: 0.01 step: 1 norm-1: 0
>> t: 0.02 step: 2 norm-1: 0
>> t: 0.03 step: 3 norm-1: 0
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint
>> t: 0 step: 0 norm-1: 0
>> t: 0.01 step: 1 norm-1: 0
>> t: 0.02 step: 2 norm-1: 1
>> t: 0.03 step: 3 norm-1: 3
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint -ts_theta_adapt
>> t: 0 step: 0 norm-1: 0
>> t: 0.01 step: 1 norm-1: 0
>> t: 0.02 step: 2 norm-1: 0
>> t: 0.03 step: 3 norm-1: 0
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint -ts_theta_adapt -ts_monitor
>> 0 TS dt 0.01 time 0
>> t: 0 step: 0 norm-1: 0
>> 0 TS dt 0.01 time 0
>> 1 TS dt 0.1 time 0.01
>> t: 0.01 step: 1 norm-1: 0
>> 1 TS dt 0.1 time 0.01
>> 2 TS dt 0.1 time 0.02
>> t: 0.02 step: 2 norm-1: 0
>> 2 TS dt 0.1 time 0.02
>> 3 TS dt 0.1 time 0.03
>> t: 0.03 step: 3 norm-1: 0
>> 3 TS dt 0.1 time 0.03
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint -ts_theta_adapt -ts_monitor -ts_adapt_monitor
>> 0 TS dt 0.01 time 0
>> t: 0 step: 0 norm-1: 0
>> 0 TS dt 0.01 time 0
>>        TSAdapt 'basic': step   0 accepted t=0          + 1.000e-02 wlte=    0 family='theta' scheme=0:'(null)' dt=1.000e-01
>> 1 TS dt 0.1 time 0.01
>> t: 0.01 step: 1 norm-1: 0
>> 1 TS dt 0.1 time 0.01
>>        TSAdapt 'basic': step   1 rejected t=0.01       + 1.000e-01 wlte=1.24e+03 family='theta' scheme=0:'(null)' dt=1.000e-02
>>        TSAdapt 'basic': step   1 accepted t=0.01       + 1.000e-02 wlte=    0 family='theta' scheme=0:'(null)' dt=1.000e-01
>> 2 TS dt 0.1 time 0.02
>> t: 0.02 step: 2 norm-1: 0
>> 2 TS dt 0.1 time 0.02
>>        TSAdapt 'basic': step   2 rejected t=0.02       + 1.000e-01 wlte=1.24e+03 family='theta' scheme=0:'(null)' dt=1.000e-02
>>        TSAdapt 'basic': step   2 accepted t=0.02       + 1.000e-02 wlte=    0 family='theta' scheme=0:'(null)' dt=1.000e-01
>> 3 TS dt 0.1 time 0.03
>> t: 0.03 step: 3 norm-1: 0
>> 3 TS dt 0.1 time 0.03
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type beuler
>> t: 0 step: 0 norm-1: 0
>> t: 0.01 step: 1 norm-1: 0
>> t: 0.02 step: 2 norm-1: 0
>> t: 0.03 step: 3 norm-1: 0
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>> $ ./ex2  -ts_type euler
>> t: 0 step: 0 norm-1: 0
>> t: 0.01 step: 1 norm-1: 0
>> t: 0.02 step: 2 norm-1: 0
>> t: 0.03 step: 3 norm-1: 0
>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug
>>
>>
>>> On Mar 20, 2015, at 10:18 PM, Andrew Spott <ansp6066 at colorado.edu> wrote:
>>>
>>> here are the data files.
>>>
>>> dipole_matrix.dat:
>>> https://www.dropbox.com/s/2ahkljzt6oo9bdr/dipole_matrix.dat?dl=0
>>>
>>> energy_eigenvalues_vector.dat
>>> https://www.dropbox.com/s/sb59q38vqvjoypk/energy_eigenvalues_vector.dat?dl=0
>>>
>>> -Andrew
>>>
>>>
>>>
>>> On Fri, Mar 20, 2015 at 7:25 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>>>
>>> Data files are needed
>>>
>>> PetscViewerBinaryOpen( PETSC_COMM_WORLD, "hamiltonian/energy_eigenvalues_vector.dat", FILE_MODE_READ, &view );
>>> VecLoad( H0, view );
>>> PetscViewerBinaryOpen( PETSC_COMM_WORLD, "hamiltonian/dipole_matrix.dat", FILE_MODE_READ, &view );
>>>
>>> BTW: You do not need to call Mat/VecAssembly on Mats and Vecs after they have been loaded.
>>>
>>> Barry
>>>
>>>
>>>> On Mar 20, 2015, at 6:39 PM, Andrew Spott <ansp6066 at colorado.edu> wrote:
>>>>
>>>> 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/20150321/39111395/attachment.html>

More information about the petsc-users mailing list