[petsc-users] TimeStepper norm problems. EMIL Please read this
Barry Smith
bsmith at mcs.anl.gov
Sat Mar 21 11:26:00 CDT 2015
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 <emconsta at mcs.anl.gov> 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 <ansp6066 at 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 <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
>>>>>
>>>>
>>>>
>>>>
>>>
>>>
>>>
>>
More information about the petsc-users
mailing list