Thanks! Let me know if I can do anything to help.
-Andrew — Andrew On Sun, Mar 22, 2015 at 8:21 PM, Emil Constantinescu <[email protected]> wrote: > Hi Andrew, > I can reproduce this issue and I agree that there is something wrong. > I'll look into it. > Emil > On 3/22/15 3:29 PM, Andrew Spott wrote: >> So, I’m now even more confused. >> >> I’m attempting to solve an equation that looks like this: >> >> u’ = -i(H0 + e(t) D) u >> >> Where H0 is a purely real diagonal matrix, D is an off-diagonal block >> matrix, and e(t) is a function of time (The schrödinger equation in the >> energy basis). >> >> I’ve rewritten the e(t) function in my code to just return 0.0. So the >> new equation is just u’ = -iH0 u. The matrix is time independent and >> diagonal (I’ve checked this). H0[0] ~= -.5 (with no imaginary >> component). and u(t=0) = [1,0,0,0,..] >> >> This problem SHOULD be incredibly simple: u’ = i (0.5) u. >> >> However, I’m still getting the same blowup with the TS.: >> >> //with e(t) == 0 >> //TS >> t: 0 step: 0 norm-1: 0 >> t: 0.01 step: 1 norm-1: 0 >> t: 0.02 step: 2 norm-1: 0.9999953125635765 >> t: 0.03 step: 3 norm-1: 2.999981250276277 >> //Hand rolled >> t: 0.01 norm-1: 0 ef 0 >> t: 0.02 norm-1: 0 ef 0 >> t: 0.03 norm-1: -1.110223024625157e-16 ef 0 >> —————————————————————————————— >> //with e(t) != 0 >> //TS >> t: 0 step: 0 norm-1: 0 >> t: 0.01 step: 1 norm-1: 0 >> t: 0.02 step: 2 norm-1: 0.9999953125635765 >> t: 0.03 step: 3 norm-1: 2.999981250276277 >> //Hand rolled >> t: 0.01 norm-1: 0 ef 9.474814647559372e-11 >> t: 0.02 norm-1: 0 ef 7.57983838406065e-10 >> t: 0.03 norm-1: -1.110223024625157e-16 ef 2.558187954267552e-09 >> >> I’ve updated the gist. >> >> -Andrew >> >> >> >> On Fri, Mar 20, 2015 at 9:57 PM, Barry Smith <[email protected] >> <mailto:[email protected]>> 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 >> <[email protected]> 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 <[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 >> > > > >> > > >> > > >> > > >> > >> > >> > >> >>
