Depending on your C++ compiler you may need to build PETSc with the additional ./configure option --with-cxx-dialect=C++11
> On Mar 21, 2015, at 9:32 AM, Emil Constantinescu <[email protected]> wrote: > > 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-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 >>>>> >>>> >>>> >>>> >>> >>> >>> >>
