[Gnucap-devel] incmode, orderings, convergence

[Felix Salfelder]

On Fri, Nov 02, 2018 at 09:22:35PM -0400, al davis wrote:
> On Sun, 28 Oct 2018 11:49:19 +0100
> Felix Salfelder <address@hidden> wrote:
> > the incmode hack was meant to make the solver more robust overall.
> > my other idea was to set non-incmode only in TRLOW and also in TRHIGH,
> > which is usually not reached. this curiously helped produce results in
> > this example, but took more iterations than 6 12 18..
> 
> TRHIGH almost makes sense.  Perhaps it should be iteration number >
> itl[TRHIGH].
> 
> This would be impossible in transient, except for initial DC.  In
> transient, a step rejection sets incmode to "BAD", which means to
> revert to non-incmode on the next step.   So this is already done.

> In DC (which includes the initial DC of a transient analysis, an
> iteration count this high often means trouble, so it sort of makes
> sense.  It might also make sense to revert to non-incmode when in
> homotopy mode, but also it might not.  It might also be interesting to
> play with damping.  (options dampmax, dampmin, dampstrategy, vmin,
> vmax).

I can see the mark_inc_mode_bad in s_tr_swp.cc +473.  somehow the
additional non-incmode in itl[TRHIGH] still made a difference. will
revisit this.

> > I kind of agree that this is the wrong approach. But it's not just the
> > switch as you suggest. if i remove S3 and replace S4 with a resistor,
> > r=1e-10, it does not converge (including wrong results) in the initial
> > dc. not even with r=1e-5 (OPT::short). but finally with r=1e-3.

> OK .. then it's the resistor!!  Really, as you verified, it is a matter
> of numerics.  With a wide spread of admittances, values, whatever,
> numeric precision becomes an issue.  When you are getting strange
> results that you cannot explain, consider that numeric precision is
> likely to be the problem.

> Most of the time, failure to converge means the circuit does not work.
> In DC, it usually means there is no stable operating point.  In this
> case, failure to converge is the correct result.  Any result that looks
> stable is incorrect.

the circuit is too simple for this.

> OPT::short is only used to substitute for a resistance of 0.  It is not
> an "Rmin".  If you explicitly specify R= (smaller than OPT::short) it
> will be left as you specify.

this is not happening in the switch. it has a hardwired value for
OPT::short. is it intentional or an oversight?

> > Also interesting: the node order makes a difference. i remember some
> > experiments with alternative static orders in -uf, nothing really worked
> > perfectly yet. i will have to port it to know better. obviously,
> > changing the component order by hand gets more tricky when using
> > schematic drawings eg. gnucsator.
> 
> yes it does .. numeric precision again.  The problem is figuring out up
> front what order works better, in an efficient way.  The cost of
> finding a better order is likely to exceed the savings of having the
> better order, and the order that looks better in one way may be worse
> in another way.
> 
> Example ...
> suppose I want to solve ...  x = a+b-c;
> That could be (a+b)-c  or a+(b-c) .
> In theory, both should give the same results, but depending on the
> values of a, b, and c one or the other will be better, and I cannot say
> in general which one.

yes. it's a tradeoff between numerical stability and speed. for once it
takes time to compute the ordering. then, depending on when this
happens, the matrix would have to be refilled, or shuffled accordingly.

> Originally, I had intended to make the ordering algorithm pluggable.

but it's a bit tricky. currently, the ordering is computed before the
components have announced their needs (iwant_matrix calls). also, the
incidence graph is not stored, only a lower bound for each spike

a good ordering could be computed (pretty quickly) from the incidence
graph and a loaded matrix. but if it is done only after load, then the
matrix rebuild will take some extra time.

all this needs to be done in plugin space. and it needs to support
solutions in between the above.

> I considered making the matrix solver pluggable, but there is too much
> speed penalty to do that.  The speed penalty of plugins is that
> associated with virtual functions, which are real function calls that
> cannot be inline.  That would add a real function call to every
> external matrix access.  If you compile with -O0 (no optimization) you
> get real function calls, and a 2x slower run time on average.

doing *_load through raw double pointers instead of function calls with
matrix coordinates could make pluggable solvers feasible. there would be
some additional overhead in storing these pointers associated to
branches alongside the matrix coordinates.

a real use for this would be to compare matrix solvers easily. but there
seems to be none avaliable. KLU (the one most hyped for circuit
simulation) does not implement partial updates..

ignoring speed penalty for a moment, i have hacked s__solve.cc to pass
the numerics to UMFPACK. a corresponding transient command s_tr_umf.cc
in [1] now helps with attributing convergence issues.

cheers
felix

[1] https://git.savannah.gnu.org/cgit/gnucap/gnucap-plugins.git/log/?h=sim