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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> 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 <[email protected]> >> > wrote: >> > >> > On Fri, Mar 20, 2015 at 3:09 PM, Andrew Spott <[email protected]> >> > 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 <[email protected]>, >> > 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" <[email protected]> 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 >> > >> >> >>
