<div>Thanks! Let me know if I can do anything to help.</div>
<div><br></div>
<div>-Andrew</div>
<div class="mailbox_signature">
<br>—<br>Andrew</div>
<br><br><div class="gmail_quote"><p>On Sun, Mar 22, 2015 at 8:21 PM, Emil Constantinescu <span dir="ltr"><<a href="mailto:emconsta@mcs.anl.gov" target="_blank">emconsta@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>Hi Andrew,
<br><br>I can reproduce this issue and I agree that there is something wrong. 
<br>I'll look into it.
<br><br>Emil
<br><br><br>On 3/22/15 3:29 PM, Andrew Spott wrote:
<br>> So, I’m now even more confused.
<br>>
<br>> I’m attempting to solve an equation that looks like this:
<br>>
<br>> u’ = -i(H0 + e(t) D) u
<br>>
<br>> Where H0 is a purely real diagonal matrix, D is an off-diagonal block
<br>> matrix, and e(t) is a function of time (The schrödinger equation in the
<br>> energy basis).
<br>>
<br>> I’ve rewritten the e(t) function in my code to just return 0.0.  So the
<br>> new equation is just u’ = -iH0 u.  The matrix is time independent and
<br>> diagonal (I’ve checked this).  H0[0] ~= -.5 (with no imaginary
<br>> component).  and u(t=0) = [1,0,0,0,..]
<br>>
<br>> This problem SHOULD be incredibly simple: u’ = i (0.5) u.
<br>>
<br>> However, I’m still getting the same blowup with the TS.:
<br>>
<br>> //with e(t) == 0
<br>> //TS
<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.9999953125635765
<br>> t: 0.03 step: 3 norm-1: 2.999981250276277
<br>> //Hand rolled
<br>> t: 0.01 norm-1: 0 ef 0
<br>> t: 0.02 norm-1: 0 ef 0
<br>> t: 0.03 norm-1: -1.110223024625157e-16 ef 0
<br>> ——————————————————————————————
<br>> //with e(t) != 0
<br>> //TS
<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.9999953125635765
<br>> t: 0.03 step: 3 norm-1: 2.999981250276277
<br>> //Hand rolled
<br>> t: 0.01 norm-1: 0 ef 9.474814647559372e-11
<br>> t: 0.02 norm-1: 0 ef 7.57983838406065e-10
<br>> t: 0.03 norm-1: -1.110223024625157e-16 ef 2.558187954267552e-09
<br>>
<br>> I’ve updated the gist.
<br>>
<br>> -Andrew
<br>>
<br>>
<br>>
<br>> On Fri, Mar 20, 2015 at 9:57 PM, Barry Smith <bsmith@mcs.anl.gov
<br>> <mailto:bsmith@mcs.anl.gov>> wrote:
<br>>
<br>>
<br>>     Andrew,
<br>>
<br>>     I'm afraid Emil will have to take a look at this and explain it. The
<br>>     -ts_type beuler and -ts_type theta -ts_theta_theta .5 are stable but
<br>>     the -ts_type cn is not stable. It turns out that -ts_type cn is
<br>>     equivalent to -ts_type theta -ts_theta_theta .5 -ts_theta_endpoint
<br>>     and somehow this endpoint business (which I don't understand) is
<br>>     causing a problem. Meanwhile if I add -ts_theta_adapt to the
<br>>     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
<br>>     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
<br>>     -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
<br>>     -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
<br>>     -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
<br>>     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
<br>>     family='theta' scheme=0:'(null)' dt=1.000e-02
<br>>     TSAdapt 'basic': step 1 accepted t=0.01 + 1.000e-02 wlte= 0
<br>>     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
<br>>     family='theta' scheme=0:'(null)' dt=1.000e-02
<br>>     TSAdapt 'basic': step 2 accepted t=0.02 + 1.000e-02 wlte= 0
<br>>     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
<br>>     <ansp6066@colorado.edu> wrote:
<br>>      >
<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>>      >
<br>>     https://www.dropbox.com/s/sb59q38vqvjoypk/energy_eigenvalues_vector.dat?dl=0
<br>>
<br>>      >
<br>>      > -Andrew
<br>>      >
<br>>      >
<br>>      >
<br>>      > On Fri, Mar 20, 2015 at 7:25 PM, Barry Smith <bsmith@mcs.anl.gov>
<br>>     wrote:
<br>>      >
<br>>      > Data files are needed
<br>>      >
<br>>      > PetscViewerBinaryOpen( PETSC_COMM_WORLD,
<br>>     "hamiltonian/energy_eigenvalues_vector.dat", FILE_MODE_READ, &view );
<br>>      > VecLoad( H0, view );
<br>>      > PetscViewerBinaryOpen( PETSC_COMM_WORLD,
<br>>     "hamiltonian/dipole_matrix.dat", FILE_MODE_READ, &view );
<br>>      >
<br>>      > BTW: You do not need to call Mat/VecAssembly on Mats and Vecs
<br>>     after they have been loaded.
<br>>      >
<br>>      > Barry
<br>>      >
<br>>      >
<br>>      > > On Mar 20, 2015, at 6:39 PM, Andrew Spott
<br>>     <ansp6066@colorado.edu> wrote:
<br>>      > >
<br>>      > > Sorry it took so long, I wanted to create a “reduced” case
<br>>     (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
<br>>     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>>      > >
<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
<br>>     <bsmith@mcs.anl.gov> wrote:
<br>>      > >
<br>>      > > Andrew,
<br>>      > >
<br>>      > > Send your entire code. It will be easier and faster than
<br>>     talking past each other.
<br>>      > >
<br>>      > > Barry
<br>>      > >
<br>>      > > > On Mar 20, 2015, at 5:00 PM, Andrew Spott
<br>>     <ansp6066@colorado.edu> wrote:
<br>>      > > >
<br>>      > > > I’m sorry, I’m not trying to be difficult, but I’m not
<br>>     following.
<br>>      > > >
<br>>      > > > The manual states (for my special case):
<br>>      > > > • u ̇ = A(t)u. Use
<br>>      > > >
<br>>      > > > TSSetProblemType(ts,TS LINEAR);
<br>>     TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,NULL);
<br>>     TSSetRHSJacobian(ts,A,A,YourComputeRHSJacobian,&appctx);
<br>>      > > >
<br>>      > > > where YourComputeRHSJacobian() is a function you provide that
<br>>     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>>      > > >
<br>>      > > > SHOULD `func` touch U? If so, what should `func` do to U? I
<br>>     thought that the RHSJacobian function was only meant to create A,
<br>>     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
<br>>     <knepley@gmail.com> wrote:
<br>>      > > >
<br>>      > > > On Fri, Mar 20, 2015 at 3:09 PM, Andrew Spott
<br>>     <ansp6066@colorado.edu> wrote:
<br>>      > > > So, it doesn’t seem that zeroing the given vector in the
<br>>     function passed to TSSetRHSJacobian is the problem. When I do that,
<br>>     it just zeros out the solution.
<br>>      > > >
<br>>      > > > I would think you would zero the residual vector (if you add
<br>>     to it to construct the residual, as in FEM methods), not the solution.
<br>>      > > >
<br>>      > > > The function that is passed to TSSetRHSJacobian has only one
<br>>     responsibility — to create the jacobian — correct? In my case this
<br>>     is A(t). The solution vector is given for when you are solving
<br>>     nonlinear problems (A(t) also depends on U(t)). In my case, I don’t
<br>>     even look at the solution vector (because my A(t) doesn’t depend on
<br>>     it).
<br>>      > > >
<br>>      > > > Are you initializing the Jacobian to 0 first?
<br>>      > > >
<br>>      > > > Thanks,
<br>>      > > >
<br>>      > > > Matt
<br>>      > > >
<br>>      > > > Is this the case? or is there some other responsibility of
<br>>     said function?
<br>>      > > >
<br>>      > > > -Andrew
<br>>      > > >
<br>>      > > > >Ah ha!
<br>>      > > > >
<br>>      > > > >The function passed to TSSetRHSJacobian needs to zero the
<br>>     solution vector?
<br>>      > > > >
<br>>      > > > >As a point, this isn’t mentioned in any documentation that I
<br>>     can find.
<br>>      > > > >
<br>>      > > > >-Andrew
<br>>      > > >
<br>>      > > > On Friday, Mar 20, 2015 at 2:17 PM, Matthew Knepley
<br>>     <knepley@gmail.com>, wrote:
<br>>      > > > This sounds like a problem in your calculation function where
<br>>     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"
<br>>     <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
<br>>     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>>      > > >
<br>>      > > > When I attempt to timestep using PETSc, the norm of `u` blows
<br>>     up. When I do it directly (using the above), the norm of `u` doesn’t
<br>>     blow up.
<br>>      > > >
<br>>      > > > It is important to note that the solution generated after the
<br>>     first step is identical for both, but the second step for Petsc has
<br>>     a norm of ~2, while for the directly calculated version it is ~1.
<br>>     The third step for petsc has a norm of ~4, while the directly
<br>>     calculated version it is still ~1.
<br>>      > > >
<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>>      > > >
<br>>      > > > `func` just constructs A(t) at the time given. The same code
<br>>     for calculating A(t) is used in both calculations, along with the
<br>>     same initial vector psi0, and the same time steps.
<br>>      > > >
<br>>      > > > Let me know what other information is needed. I’m not sure
<br>>     what could be the problem. `func` doesn’t touch U at all (should it?).
<br>>      > > >
<br>>      > > > -Andrew
<br>>      > > >
<br>>      > > >
<br>>      > > >
<br>>      > > >
<br>>      > > > --
<br>>      > > > What most experimenters take for granted before they begin
<br>>     their experiments is infinitely more interesting than any results to
<br>>     which their experiments lead.
<br>>      > > > -- Norbert Wiener
<br>>      > > >
<br>>      > >
<br>>      > >
<br>>      > >
<br>>      >
<br>>      >
<br>>      >
<br>>
<br>>
<br></p></blockquote></div><br>