I haven't been able to compile and run. But here are a few quick notes.
The problem appears to be very stiff. Theta and theta_endpoint are defining different methods: 1) -ts_type beuler OR -ts_theta_theta 1.0: is backward Euler u(t + h) = u(t) + h*A(t+h)*u(t+h) 2) -ts_theta_theta 0.5: is the implicit midpoint rule u(t + h) = u(t) + h*[A(t+h/2)*(u(t+h)+u(t))/2]3) -ts_type cn OR -ts_theta_theta 0.5 -ts_theta_endpoint: is Crank-Nicholson/trapezoidal
u(t + h) = u(t) + h/2[A(t+h)*u(t+h) + A(t)*u(t)]Note that the last two are different. -ts_type theta -ts_theta_theta .5 is different from -ts_type cn. They the same linear stability properties if A(t)=A; but not if A depends on t.
When -ts_theta_adapt is used, then it detects the instability as an error and reduces the step by a lot! wlte=1.24e+03 which means that the reduction should be severe but the controller tries 0.1*dt and that 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. Emil On 3/20/15 10:57 PM, Barry Smith 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-debugOn 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. BarryOn 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. BarryOn 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? -AndrewAh 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. -AndrewOn 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
