[petsc-users] TimeStepper norm problems.

Andrew Spott ansp6066 at colorado.edu
Fri Mar 20 22:18:55 CDT 2015 

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 <bsmith at 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 <ansp6066 at 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 <bsmith at 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 <ansp6066 at 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 <knepley at gmail.com> wrote: 
>> > 
>> > On Fri, Mar 20, 2015 at 3:09 PM, Andrew Spott <ansp6066 at 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 <knepley at 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" <ansp6066 at 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 
>> > 
>> 
>> 
>> 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150320/ef0fae6d/attachment.html>

More information about the petsc-users mailing list