<span id="mailbox-conversation"><div>>I haven’t been able to compile and run.</div>

<div><br></div>

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

<div><br></div>

<div><div>>The problem appears to be very stiff.</div></div>

<div><br></div>

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

<div><br></div>

<div>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?</div>

<div><br></div>

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

<div class="gmail_quote">On Sat, Mar 21, 2015 at 8:32 AM, Emil Constantinescu <span dir="ltr"><<a href="mailto:emconsta@mcs.anl.gov" target="_blank">emconsta@mcs.anl.gov</a>> wrote: <blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex;">I haven't been able to compile and run. But here are a few quick notes.

The problem appears to be very stiff.

<br><br>Theta and theta_endpoint are defining different methods:

<br><br>1) -ts_type beuler OR -ts_theta_theta 1.0: is backward Euler

<br><br>u(t + h) = u(t) + h*A(t+h)*u(t+h)

<br><br>2) -ts_theta_theta 0.5: is the implicit midpoint rule

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

<br><br>3) -ts_type cn OR -ts_theta_theta 0.5 -ts_theta_endpoint: is 

<br>Crank-Nicholson/trapezoidal

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

<br><br>Note that the last two are different. -ts_type theta -ts_theta_theta .5 

<br>is different from -ts_type cn. They the same linear stability properties 

if A(t)=A; but not if A depends on t.

<br><br>When -ts_theta_adapt is used, then it detects the instability as an 

<br>error and reduces the step by a lot! wlte=1.24e+03 which means that the 

<br>reduction should be severe but the controller tries 0.1*dt and that 

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

<br><br>Emil

<br><br>On 3/20/15 10:57 PM, Barry Smith wrote:

<br>>

<br>>    Andrew,

<br>>

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

<br>>

<br>>    Emil,

<br>>

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

<br>>

<br>>    Barry

<br>>

<br>>

<br>> $ ./ex2  -ts_type cn

<br>> t: 0 step: 0 norm-1: 0

<br>> t: 0.01 step: 1 norm-1: 0

<br>> t: 0.02 step: 2 norm-1: 1

<br>> t: 0.03 step: 3 norm-1: 3

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type theta

<br>> t: 0 step: 0 norm-1: 0

<br>> t: 0.01 step: 1 norm-1: 0

<br>> t: 0.02 step: 2 norm-1: 0

<br>> t: 0.03 step: 3 norm-1: 0

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type theta -ts_theta_theta .5

<br>> t: 0 step: 0 norm-1: 0

<br>> t: 0.01 step: 1 norm-1: 0

<br>> t: 0.02 step: 2 norm-1: 0

<br>> t: 0.03 step: 3 norm-1: 0

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint

<br>> t: 0 step: 0 norm-1: 0

<br>> t: 0.01 step: 1 norm-1: 0

<br>> t: 0.02 step: 2 norm-1: 1

<br>> t: 0.03 step: 3 norm-1: 3

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint -ts_theta_adapt

<br>> t: 0 step: 0 norm-1: 0

<br>> t: 0.01 step: 1 norm-1: 0

<br>> t: 0.02 step: 2 norm-1: 0

<br>> t: 0.03 step: 3 norm-1: 0

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint -ts_theta_adapt -ts_monitor

<br>> 0 TS dt 0.01 time 0

<br>> t: 0 step: 0 norm-1: 0

<br>> 0 TS dt 0.01 time 0

<br>> 1 TS dt 0.1 time 0.01

<br>> t: 0.01 step: 1 norm-1: 0

<br>> 1 TS dt 0.1 time 0.01

<br>> 2 TS dt 0.1 time 0.02

<br>> t: 0.02 step: 2 norm-1: 0

<br>> 2 TS dt 0.1 time 0.02

<br>> 3 TS dt 0.1 time 0.03

<br>> t: 0.03 step: 3 norm-1: 0

<br>> 3 TS dt 0.1 time 0.03

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type theta -ts_theta_theta .5  -ts_theta_endpoint -ts_theta_adapt -ts_monitor -ts_adapt_monitor

<br>> 0 TS dt 0.01 time 0

<br>> t: 0 step: 0 norm-1: 0

<br>> 0 TS dt 0.01 time 0

<br>>        TSAdapt 'basic': step   0 accepted t=0          + 1.000e-02 wlte=    0 family='theta' scheme=0:'(null)' dt=1.000e-01

<br>> 1 TS dt 0.1 time 0.01

<br>> t: 0.01 step: 1 norm-1: 0

<br>> 1 TS dt 0.1 time 0.01

<br>>        TSAdapt 'basic': step   1 rejected t=0.01       + 1.000e-01 wlte=1.24e+03 family='theta' scheme=0:'(null)' dt=1.000e-02

<br>>        TSAdapt 'basic': step   1 accepted t=0.01       + 1.000e-02 wlte=    0 family='theta' scheme=0:'(null)' dt=1.000e-01

<br>> 2 TS dt 0.1 time 0.02

<br>> t: 0.02 step: 2 norm-1: 0

<br>> 2 TS dt 0.1 time 0.02

<br>>        TSAdapt 'basic': step   2 rejected t=0.02       + 1.000e-01 wlte=1.24e+03 family='theta' scheme=0:'(null)' dt=1.000e-02

<br>>        TSAdapt 'basic': step   2 accepted t=0.02       + 1.000e-02 wlte=    0 family='theta' scheme=0:'(null)' dt=1.000e-01

<br>> 3 TS dt 0.1 time 0.03

<br>> t: 0.03 step: 3 norm-1: 0

<br>> 3 TS dt 0.1 time 0.03

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type beuler

<br>> t: 0 step: 0 norm-1: 0

<br>> t: 0.01 step: 1 norm-1: 0

<br>> t: 0.02 step: 2 norm-1: 0

<br>> t: 0.03 step: 3 norm-1: 0

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>> $ ./ex2  -ts_type euler

<br>> t: 0 step: 0 norm-1: 0

<br>> t: 0.01 step: 1 norm-1: 0

<br>> t: 0.02 step: 2 norm-1: 0

<br>> t: 0.03 step: 3 norm-1: 0

<br>> ~/Src/petsc/test-dir (barry/more-tchem-work *=) arch-debug

<br>>

<br>>

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

<br>>>

>> here are the data files.

<br>>>

<br>>> dipole_matrix.dat:

<br>>> https://www.dropbox.com/s/2ahkljzt6oo9bdr/dipole_matrix.dat?dl=0

<br>>>

<br>>> energy_eigenvalues_vector.dat

<br>>> https://www.dropbox.com/s/sb59q38vqvjoypk/energy_eigenvalues_vector.dat?dl=0

<br>>>

<br>>> -Andrew

<br>>>

<br>>>

<br>>>

<br>>> On Fri, Mar 20, 2015 at 7:25 PM, Barry Smith <bsmith@mcs.anl.gov> wrote:

<br>>>

<br>>> Data files are needed

<br>>>

<br>>> PetscViewerBinaryOpen( PETSC_COMM_WORLD, "hamiltonian/energy_eigenvalues_vector.dat", FILE_MODE_READ, &view );

<br>>> VecLoad( H0, view );

<br>>> PetscViewerBinaryOpen( PETSC_COMM_WORLD, "hamiltonian/dipole_matrix.dat", FILE_MODE_READ, &view );

<br>>>

>> BTW: You do not need to call Mat/VecAssembly on Mats and Vecs after they have been loaded.

<br>>>

<br>>> Barry

<br>>>

<br>>>

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

<br>>>>

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

<br>>>>

<br>>>> https://gist.github.com/spott/aea8070f35e79e7249e6

<br>>>>

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

<br>>>>

<br>>>> //first section, using TimeStepper:

<br>>>> t: 0 step: 0 norm-1: 0

<br>>>> t: 0.01 step: 1 norm-1: 0

<br>>>> t: 0.02 step: 2 norm-1: 0.999995

<br>>>> t: 0.03 step: 3 norm-1: 2.99998

<br>>>>

>>> //Second section, using explicit code.

<br>>>> t: 0.01 norm-1: 0

<br>>>> t: 0.02 norm-1: 0

<br>>>> t: 0.02 norm-1: 2.22045e-16

<br>>>>

<br>>>>

<br>>>>

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

<br>>>>

<br>>>> Andrew,

<br>>>>

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

<br>>>>

<br>>>

<br>>>

<br>>>

<br>>

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