Re: [Gnucap-devel] Integration methods

[al davis]

On Saturday 23 December 2006 09:54, nbs_r wrote:
> Which integration methods are currently implemented in
> gnucap's transient analysis? As far as I can see, only
> trapezoid and backward-euler is present (e_storag.cc), and
> gear method falls back to euler.

That is correct.

> Is there any plan to extend 
> this set?

I think it is obvious from the code that the answer is "yes".  
With the long list of things to do, it has not been a priority.

The plan is (and has been) to add Gear-2 and heuristics to pick 
between the 3 provided methods, by device, similar to 
what "Spectre" does.

At first glance, adding Gear-2 involves putting the correct 
number in 3 tables in e_storag.h, and adding the formula (one 
line of code) in e_storag.cc.  That will seem to work most of 
the time.  You also need to keep one more time step, and add 
the code to keep it updated.  (one more line, a few places).  
It seems pretty simple.

It looks easy so far, but there are a few issues that will need 
to be addressed.

How do you deal with variable time step?  BE and trap use one 
time step, so the formula is simple.  Gear-2 uses 2 time steps.  
The classic formula assumes they are equal, but they are not 
guaranteed to be equal.  If they are unequal, most likely they 
will be very unequal.  One possibility is to switch to a single 
step method for that one time step, but they you lose the 
advantage of Gear-2.  Switching to trap would not be a good 
idea because time step changes tend to occur when the ringing 
problem with trap is most pronounced.  Switching to euler would 
not be a good idea because time step changes tend to occur when 
additional accuracy is suddenly needed.  An effective value for 
the misplaced step could be generated by interpolation, but 
this has similar effects to a method change.  That leaves 
generating a new formula for unequal steps.

How do you control the time step?  This may seem easy, because 
in theory the same algorithm as for trap works, with only a 
coefficient change.  In practice, it doesn't work that way.  
For the same truncation error, you need smaller steps with Gear 
than with Trap.  If the step control algorithm is working 
correctly, Trap ringing will not be a problem, so you might as 
well just use Trap.  Considering that step control might not 
always work correctly, that is where Gear pays off, but often 
just using smaller steps with Trap works as well, usually with 
more accuracy.

The real issue with time step control is that you don't want to 
use the same algorithm.  The reason for a stiffly stable method 
is to knowingly ignore the stiff poles.  The best step control 
with Gear may be "none", as with Euler.  Actually, there still 
is step control based on events.

To take this a bit farther, it might make sense to use 
truncation error to select the method.  If the error is low, 
use Trap.  If the error is borderline, use Gear-2.  If the 
error is high use Euler.  It might be better to continuously 
vary between the 3 methods.  In all cases, leave the step size 
constant, or maybe tweek it to keep the leading coefficient 
constant.  I am not sure how effective this would be.  It is 
worth thinking about.  


> Also, I would like to hack some other algorithms. How to test
> the correctness and performance of the integration method?
> Is there any testbench available reliably exercising
> integration accuracy and stability? If not, how to isolate
> this part of analysis? I.e. how to kill time-step control and
> save all transient simulation data.

That is worth doing.  I have never seen a truly good writeup of 
this, but there are some that will get you going.

Testing will take more of your time than algorithm development.  
(Remember ... it is just a few lines.)  I have a few good test 
cases, but you will need more.

Gnucap does not have a command to turn off time-step control 
completely, but you can get the effect in several ways.  One 
way is to set "trtol" and "chgtol" to a high number.  In the 
code, you can make STORAGE::tr_review set  _time_future to 
NEVER, and make it return (NEVER, NEVER).  This will not 
disable step control from other devices, but it will probably 
do what you want.

Saving all transient data is easy ...  add "trace all" to the 
transient command.