From: Soeren D. S. <soe...@gm...> - 2013-08-26 18:29:29
|
Hi again, I created a new branch called "numerical_improvements" which includes some changes to improve the convergence behavior of transient analysis. Summary of the main changes: * In the corrector and the non-linear solver, the implementation of Newton's method was changed to an "update" style rather than computing the solution directly. This leads to a number of advantages. First, the achievable accuracy is higher. That's because the updates between the integration steps are rather small, so a specific relative error (which is usually determined by the condition of the Jacobian) has a smaller impact when it's on the update rather than on the solution itself. The program is less likely to experience Newton rejections upon small step size. Additionally, the change enables the use of the so-called "simplified Newton method" which may save some computation time. * Specifically, energy-storage components are no longer modeled as "conductance + current source", but their values are directly written into the so-called "F matrix" which is later combined with the A matrix according to the step size. * Source stepping is now used as a convergence helper for the initial DC analysis. All other convergence helpers are temporarily disabled. (It's possible to re-enable them, but I currently don't see any need because Newton's method converges locally.) It was a bit surprising to me that source stepping works the way it is implemented, because the code looks erroneous: After a Newton rejection, the rejected guess is used as a starting value for the new iteration. If I "fix" it, I get into a local minimum. If I leave it this way, the method seems to generate "random" guesses which eventually converge after a few 100 steps. (Maybe it's better to use actually random starting values?) * An additional rejection criterion for the Newton iteration was introduced: The iteration is rejected if the update of the solution in one step is greater than the update in the previous Newton step. A similar criterion is thinkable for the RHS, but we probably don't need two. * The Gear predictor is used for *all* integration methods. As I explained earlier, the others can never exceed order 1 because there are no precise values for the derivative of the solution available. This boosts the Trapezoidal method quite a bit, while Euler and Gear are unaffected (they coincide for order 1 anyway). Higher-order Adams methods are always bad, so a better predictor guess doesn't do too much. * e_trsolver was temporarily disabled. What exactly does it do? The changes were done in a conservative way, so old models remain working. Their performance may vary a bit due to some changes in the convergence criteria. On the other hand, almost everything was updated to the new way, except for one problem: * Non-linear capacitors (and inductors) should to be reconsidered. I fail to find any mathematical convergence analysis for what Qucs currently does. There is some literature on BDF concerning the code DASSL, but the exact algorithm is a bit different. So while the current method may work, I think what saves a lot of trouble is switching to the "charge-oriented MNA", which introduces an extra equation and an extra variable for the charge of each non-linear capacitor. In this case, the numerical properties are much better-known, and we can also plot the charge easily if it's interesting to the user. Unfortunately, such a change is not completely transparent, and it requires changes in VerilogA-generated code. Is there anyone who knows this code well and can estimate if it's easy to add an extra column and an extra row for each non-linear capacitor? What's left on my TODO list: * Clean up the code. * Document it in the technical papers. * Speed it up. I'm currently following "premature optimization is the root of all evil". My code introduces a bit of CPU time and memory overhead. While matrix decomposition is still slow, it doesn't really matter, but it shouldn't be there. On the other hand, implementing the simplified Newton method might reduce the number of necessary decompositions by a factor of 2 immediately. * Implement the Radau5 one-step method. Just implementing the corrector isn't hard, but making the method work well is. Error estimation is based on an embedded order 4 method, the predictor extrapolates an internal polynomial ("collocation polynomial") from the previous step, and an *efficient* solver has to solve both a real-valued and a complex-valued matrix. This is a bit hard to realize using the current abstractions in "nasolver". (Without this optimization, LU decomposition is slower by a factor of 9.) In the original code, a Gustaffson controller is used for step size control. * Implement the method for improving the condition of the Jacobian that I developed in my Bachelor's thesis. If you want to help, please do some tests (I already tested everything except for the n-ary inductor myself). Tell me if the code works as I described, or if things stop working for you. As a benchmark problem, I attached the circuit of the "Blameless Amplifier" from "Douglas Self: Audio Power Amplifier Design Handbook". It doesn't work in the old code primarily because no source stepping is used in the Initial DC analysis. With the new code, Euler and Trapezoidal work better, but the difference for Gear is small. Circuits that work well with the new code and badly with the old code (or vice versa) would be very interesting for me. Sören |
From: Richard C. <r.c...@ed...> - 2013-08-27 09:34:57
|
On 26/08/2013 19:29, Soeren D. Schulze wrote: > Hi again, > > > * e_trsolver was temporarily disabled. What exactly does it do? > > e_trsolver is the nascent external interface (external trsolver) I am creating to the qucs transient solvers. It provides a means to interact with the solver between each time step, or even completely control the stepping from another program when qucs is compiled as a library (which has to be done manually at the moment). The external interface currently works in two ways, synchonous or asynchronous. In the asynchonous version it is intended that two time steps will be provided and the qucs solvers solve between these as they normally do, choosing their own minor time steps and simply reporting the values of the solution at the end. In the synchronous version however it is intended that all step control is performed by the external program, using the previous solution values reported by qucs. I have a prototype of this working using the Matlab/Octave solvers, which determine the local gradient by fitting a 3-point spline to the previous solution values. This is crude but seems to work ok for testing. Richard -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. |
From: Soeren D. S. <soe...@gm...> - 2013-08-27 10:43:34
|
Am 27.08.2013 11:34, schrieb Richard Crozier: > On 26/08/2013 19:29, Soeren D. Schulze wrote: >> Hi again, >> >> * e_trsolver was temporarily disabled. What exactly does it do? > > e_trsolver is the nascent external interface (external trsolver) I am > creating to the qucs transient solvers. It provides a means to interact > with the solver between each time step, or even completely control the > stepping from another program when qucs is compiled as a library (which > has to be done manually at the moment). > > The external interface currently works in two ways, synchonous or > asynchronous. In the asynchonous version it is intended that two time > steps will be provided and the qucs solvers solve between these as they > normally do, choosing their own minor time steps and simply reporting > the values of the solution at the end. In the synchronous version > however it is intended that all step control is performed by the > external program, using the previous solution values reported by qucs. I > have a prototype of this working using the Matlab/Octave solvers, which > determine the local gradient by fitting a 3-point spline to the previous > solution values. This is crude but seems to work ok for testing. Ah, I see. My main problem was the great code redundancy between e_trsolver.cpp and trsolver.cpp. Is it possible to refactor this? The e_trsolver class is a subclass of trsolver, but it removes functionality rather than adding to it. So shouldn't the hierarchy be the other way around then? Sören |
From: Richard C. <r.c...@ed...> - 2013-08-27 11:15:55
|
On 27/08/2013 11:43, Soeren D. Schulze wrote: > Am 27.08.2013 11:34, schrieb Richard Crozier: >> On 26/08/2013 19:29, Soeren D. Schulze wrote: >>> Hi again, >>> >>> * e_trsolver was temporarily disabled. What exactly does it do? >> >> e_trsolver is the nascent external interface (external trsolver) I am >> creating to the qucs transient solvers. It provides a means to interact >> with the solver between each time step, or even completely control the >> stepping from another program when qucs is compiled as a library (which >> has to be done manually at the moment). >> >> The external interface currently works in two ways, synchonous or >> asynchronous. In the asynchonous version it is intended that two time >> steps will be provided and the qucs solvers solve between these as they >> normally do, choosing their own minor time steps and simply reporting >> the values of the solution at the end. In the synchronous version >> however it is intended that all step control is performed by the >> external program, using the previous solution values reported by qucs. I >> have a prototype of this working using the Matlab/Octave solvers, which >> determine the local gradient by fitting a 3-point spline to the previous >> solution values. This is crude but seems to work ok for testing. > > Ah, I see. > > My main problem was the great code redundancy between e_trsolver.cpp and > trsolver.cpp. Is it possible to refactor this? > > The e_trsolver class is a subclass of trsolver, but it removes > functionality rather than adding to it. So shouldn't the hierarchy be > the other way around then? > > > Sören > Well, the intention of subclassing trsolver was that e_trsolver then had everything trsolver had. I wanted to reuse as much as possible of the trsolver code originally, just adding in a few more things, and overriding some of trsolvers methods. I also didn't want to modify trsolver heavily while I worked on this and figured out the best way to do it, and I also just wanted to see if it could work really. The architecture of this is still in a state of flux. I'm very busy right now, but soon I will have a bit more time to work on Qucs again, and will get this looking a bit nicer. In the meantime I'm open to suggestions/improvements. I'm not sure about trsolver inheriting from e_trsolver though, I'd rather modify trsolver where possible, as it seems more natural to me that trsolver is the parent. Richard -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. |
From: roucaries b. <rou...@gm...> - 2013-08-27 14:29:49
|
Could you test my improvement also ? i suppose my work clash with your work so I want to be merged ASAP in order to avoid conflict & rebase Bastien On Tue, Aug 27, 2013 at 1:15 PM, Richard Crozier <r.c...@ed...> wrote: > On 27/08/2013 11:43, Soeren D. Schulze wrote: >> Am 27.08.2013 11:34, schrieb Richard Crozier: >>> On 26/08/2013 19:29, Soeren D. Schulze wrote: >>>> Hi again, >>>> >>>> * e_trsolver was temporarily disabled. What exactly does it do? >>> >>> e_trsolver is the nascent external interface (external trsolver) I am >>> creating to the qucs transient solvers. It provides a means to interact >>> with the solver between each time step, or even completely control the >>> stepping from another program when qucs is compiled as a library (which >>> has to be done manually at the moment). >>> >>> The external interface currently works in two ways, synchonous or >>> asynchronous. In the asynchonous version it is intended that two time >>> steps will be provided and the qucs solvers solve between these as they >>> normally do, choosing their own minor time steps and simply reporting >>> the values of the solution at the end. In the synchronous version >>> however it is intended that all step control is performed by the >>> external program, using the previous solution values reported by qucs. I >>> have a prototype of this working using the Matlab/Octave solvers, which >>> determine the local gradient by fitting a 3-point spline to the previous >>> solution values. This is crude but seems to work ok for testing. >> >> Ah, I see. >> >> My main problem was the great code redundancy between e_trsolver.cpp and >> trsolver.cpp. Is it possible to refactor this? >> >> The e_trsolver class is a subclass of trsolver, but it removes >> functionality rather than adding to it. So shouldn't the hierarchy be >> the other way around then? >> >> >> Sören >> > > Well, the intention of subclassing trsolver was that e_trsolver then had > everything trsolver had. I wanted to reuse as much as possible of the > trsolver code originally, just adding in a few more things, and > overriding some of trsolvers methods. I also didn't want to modify > trsolver heavily while I worked on this and figured out the best way to > do it, and I also just wanted to see if it could work really. The > architecture of this is still in a state of flux. > > I'm very busy right now, but soon I will have a bit more time to work on > Qucs again, and will get this looking a bit nicer. In the meantime I'm > open to suggestions/improvements. I'm not sure about trsolver inheriting > from e_trsolver though, I'd rather modify trsolver where possible, as it > seems more natural to me that trsolver is the parent. > > Richard > > > > -- > The University of Edinburgh is a charitable body, registered in > Scotland, with registration number SC005336. > > > ------------------------------------------------------------------------------ > Introducing Performance Central, a new site from SourceForge and > AppDynamics. Performance Central is your source for news, insights, > analysis and resources for efficient Application Performance Management. > Visit us today! > http://pubads.g.doubleclick.net/gampad/clk?id=48897511&iu=/4140/ostg.clktrk > _______________________________________________ > Qucs-devel mailing list > Quc...@li... > https://lists.sourceforge.net/lists/listinfo/qucs-devel |
From: Soeren D. S. <soe...@gm...> - 2013-08-27 20:42:44
|
Am 27.08.2013 16:29, schrieb roucaries bastien: > Could you test my improvement also ? > > i suppose my work clash with your work so I want to be merged ASAP in > order to avoid conflict & rebase > > Bastien Was this in response to me? I just merged your branch into mine, and I don't see any issues. But it doesn't seem to improve anything, either, because I don't use tmatrix/tvector. I also uploaded a few speed-ups in my branch (mainly concerning matrix creation); now matrix computations are the bottleneck, so porting tmatrix/tvector to libeigen will probably cause some performance boost. Sören |
From: Kevin C. <cam...@gm...> - 2013-08-28 07:36:18
|
Just FYI I did an integration of Icarus Verilog with GNUcap and another with Spice3 for mixed signal with Icarus being loaded as a shared library - http://iverilog.wikia.com/wiki/GIT_Branch_Summary - embedded-vvp (the bulk of the code, C++) For GNUcap/Spice I just modified the PWL sources to call out to the dynamically loaded code for data points, and for Icarus I added VPI code to locate the (non-blocking) assignments and convert them to PWL data. I added truncation code in the PWL sources to catch the logic thresholds on A->D boundary. Control gets handed to the subordinate simulators between solving and final acceptance (as far as I can remember). Can probably get it working again if anyone is interested. Kev. On 08/27/2013 02:34 AM, Richard Crozier wrote: > On 26/08/2013 19:29, Soeren D. Schulze wrote: >> Hi again, >> >> >> * e_trsolver was temporarily disabled. What exactly does it do? >> >> > e_trsolver is the nascent external interface (external trsolver) I am > creating to the qucs transient solvers. It provides a means to interact > with the solver between each time step, or even completely control the > stepping from another program when qucs is compiled as a library (which > has to be done manually at the moment). > > The external interface currently works in two ways, synchonous or > asynchronous. In the asynchonous version it is intended that two time > steps will be provided and the qucs solvers solve between these as they > normally do, choosing their own minor time steps and simply reporting > the values of the solution at the end. In the synchronous version > however it is intended that all step control is performed by the > external program, using the previous solution values reported by qucs. I > have a prototype of this working using the Matlab/Octave solvers, which > determine the local gradient by fitting a 3-point spline to the previous > solution values. This is crude but seems to work ok for testing. > > Richard > > > |
From: Richard C. <r.c...@ed...> - 2013-08-28 08:03:31
|
On 28/08/2013 08:36, Kevin Cameron wrote: > > Just FYI I did an integration of Icarus Verilog with GNUcap and another with Spice3 for mixed signal with Icarus being loaded as a shared library - > > http://iverilog.wikia.com/wiki/GIT_Branch_Summary - embedded-vvp (the bulk of the code, C++) > > For GNUcap/Spice I just modified the PWL sources to call out to the dynamically loaded code for data points, and for Icarus I added VPI code to locate the (non-blocking) assignments and convert them to PWL data. I added truncation code in the PWL sources to catch the logic thresholds on A->D boundary. Control gets handed to the subordinate simulators between solving and final acceptance (as far as I can remember). > > Can probably get it working again if anyone is interested. > > Kev. > This probably is of interest, I've forwarded this to the gnucap list, I hope you don't mind. Richard -- The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336. |
From: Guilherme T. <gui...@gm...> - 2013-08-28 13:46:21
|
Hi Sören, On Mon, Aug 26, 2013 at 8:29 PM, Soeren D. Schulze <soe...@gm...>wrote: > > > <snip> > As a benchmark problem, I attached the circuit of the "Blameless > Amplifier" from "Douglas Self: Audio Power Amplifier Design Handbook". > It doesn't work in the old code primarily because no source stepping is > used in the Initial DC analysis. With the new code, Euler and > Trapezoidal work better, but the difference for Gear is small. > > The mailing list does not accept attachments. Can you make the benchmark available elsewhere? Guilherme Brondani Torri |
From: Soeren D. S. <soe...@gm...> - 2013-08-28 15:54:07
|
Am 28.08.2013 15:46, schrieb Guilherme Torri: > As a benchmark problem, I attached the circuit of the "Blameless > Amplifier" from "Douglas Self: Audio Power Amplifier Design Handbook". > It doesn't work in the old code primarily because no source stepping is > used in the Initial DC analysis. With the new code, Euler and > Trapezoidal work better, but the difference for Gear is small. > > The mailing list does not accept attachments. Can you make the benchmark > available elsewhere? Sorry, I didn't know that. Here it is: http://pastebin.com/Jv7erhMV It's from the book mentioned above (the circuit, not the .sch file), so I'm not sure if it's OK to include it in the repo. I created another circuit (a simple oscillator) which fails with the "master" version of Qucs but runs with my changes. This time, it's not because of the initial DC but actually because of the transient analysis: http://pastebin.com/9f3LDxeR This one is interesting because it only runs properly with Euler. It seems like a stability problem to me, so I'm curious if it would work with the Radau5 method (which has order 5 but the stability properties of the Euler method). Some additional notes: * It's currently impossible to specify initial conditions of capacitors/inductors manually (which is why it was hard to get the oscillator above started). I'm not entirely happy with the way it was implemented before – the initial DC solution was simply overridden. The right thing to do seems to be: model capacitors with a specified initial voltage as voltage sources in the initial DC simulation. But it's a bit non-trivial because this changes the size of the matrix. * If you enable STEPDEBUG, you will notice that the condition of the matrix becomes very bad in the second circuit if you increase the "off" resistance of the switch. This is because the resistance (not the conductance!) becomes an entry in the matrix. It's on the diagonal, so as long as we do LU decomposition, it's pretty harmless, but we should keep an eye on it. Fixing this is probably easy: just scaling the row should work. Sören |