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 <ansp6...@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 <bsm...@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 <ansp6...@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 <bsm...@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 <ansp6...@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 <knep...@gmail.com> > > > wrote: > > > > > > On Fri, Mar 20, 2015 at 3:09 PM, Andrew Spott <ansp6...@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 <knep...@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" <ansp6...@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 > > > > > > > > > > > >