Resistive residual with limiting

[arpadbuermen]

I have a problem with resistive residual under limiting. Suppose I have the following model (dumbest possible diode)

`include "constants.vams"
`include "disciplines.vams"

module simdio(A,C);
    // Simplified SPICE diode

    // Terminals and internal nodes
    inout A, C;
    electrical A,C;

    // Branches
    branch (A,C) br_a_c;

    // Model parameters
    (*desc = "Saturation current", units = "A"*) 
      parameter real Is = 1e-14 from [0:inf];
    (*desc = "Emission coefficient", units = ""*) 
      parameter real N = 1 from [0:inf];

    // Opvars
    (*desc= "Diode ohmic current", units ="A" *) real If;
    (*desc= "Differential conductance", units ="A/V" *) real gd;

    // Variables
    real Vfu, Vf, VT;
    real Tdev;

    analog begin
        // Device temperature and parameter measurement temperature (in K)
        Tdev = $temperature;

        // Thermal voltage
        VT = `P_K*Tdev/`P_Q;

        // Branch voltages
        Vfu = V(br_a_c);
        $display("Vf         %.5e", Vfu);
        Vf = $limit(V(br_a_c), "pnjlim", $vt, 0.7);
        $display("Vf limited %.5e", Vf);

        // Ohmic current of diode branch
        If = Is * (exp(Vf/(N*VT))-1);

        gd = ddx(If, V(A));

        $display("If limited %.5e", If);
        $display("gd limited %.5e", gd);

        // Diode branch current
        I(br_a_c) <+ If;
    end
endmodule

I build a circuit with a voltage source (2V) powering a diode. Spice3-like netlist

V1 (1 0) dc=2
D1 (1 0) dumbdiode

The old solution at iteration 1 is set to v(1)=1, i(V1)=0

I ran eval() with the following flags CALC_RESIST_RESIDUAL | CALC_RESIST_JACOBIAN | CALC_OP | ENABLE_LIM | CALC_RESIST_LIM_RHS | ANALYSIS_DC | ANALYSIS_STATIC | INIT_LIM (in first iteration that I described, in later iterations this flag is off)

At first iteration I get OSDI(debug) d1: Vf 1.00000e+00 OSDI(debug) d1: Vf limited 7.00000e-01 OSDI(debug) d1: If limited 8.22997e+00 OSDI(debug) d1: gd limited 3.49642e+02

Obviously limiting reduced the diode voltage from 1.0 to 0.7.

This is the resistive residual loaded by load_residual_resist() 1 : 8.22997 v1:flow(br) : 1 This looks fine. The only thing that contributes a residual to node 1 is the dode. The voltage source current is 0 so its contribution to node 1 residual should be 0. The residual at node 1 is equal to the diode current at limited Vf.

Now this is the limited residual loaded by load_limit_rhs_resist() 1 : -349.642 v1:flow(br) : 0

According to the diode example in OSDI manual I expected the limited residual to be gd(Vflimited-Vf) which should result in -0.3349.6=-104.9. Is something wrong with load_limit_rhs_resist() or am I getting something wrong?

Secondly, if I dump the resistive residual and the limited resistive residual from the diode instance via the offsets specified in the device descriptor (resist_residual_off, resist_limit_rhs_off) I get the following residual contributed by the diode anode d1, A residual=8.22997 The limited residual for anode is not available (resist_limit_rhs_off==UINT32_MAX). I expected resist_limit_rhs_off to be different from UINT32_MAX and at least provide the value loaded by load_limit_rhs_resist().

[pascalkuthe]

thanks for reporting this! You are right that the expected value would be gd*(Vflimited-Vf). I looked trough the code and the implementation does seem correct to me and I did successfully test a similar example with ngspice when I added this feature.

I will probably need to add some more robust testing of the OSDI output files to figure this case out (I can only really test them right now by testing ngspice/melange). That will probably require some extra work to setup. I am currently focused on finishing the noise integration so it may take a bit of time until I can look at something larger like that.

[arpadbuermen]

I will prepare an example with NGSpice patches that demostrates the problem as soon as I have time.

I have another comment regarding residual and limited residual.

load_residual_resist() loads the resistive residual that has optionally been limited. In the Netwon-Raphson iteration the sign of this contribution changes as it goes on the RHS of the linear system, provided that we want to solve for xdelta which is used then in computing the new candidate solution as xnew = xold + xdelta

On the other hand load_limit_rhs_resist() (which computes the linearized component of the residual when limiting is applied) produces the contribution to the RHS (provided that we are solving for xdelta). The sign of the load_residual_resist contribution and load_limit_rhs_resist contribution differs. Therefore we must either 1) call load_residual_resist(x) for all instances, then change the sign of x, and finally, call load_limit_rhs_resist(x) 2) call load_residual_resist(x) followed by load_limit_rhs_resist(y) for each instance, and then compute the residual as x-y (or RHS as y-x)

Neither of the two options is optimal 1) is bad for cache locality, while 2) requires an extra vector

It would make sense if load_residual_resist and load_limit_rhs_resist would agree in the sign of their contribution, or at least provide a function load_limit_residual_resist() that would load the linearized residual contributions in the same way as load_residual_resist() does.

[pascalkuthe]

thanks for the reports here. This has been fixed as far as I am willing to in #92.

The original bug report about noise was totally correct and a large oversight thanks for that ! I now have testcases to ensure these functienos work correctly.

Regarding the second point I agree that the current function was incorrect (it was actually a bug/oversight). The sign of load_limit_rhs was supposed to be negative.

While it is true that the newton algorithm is often written in mathematical literator as:

$$ J_f \Delta x = -f(x) \quad x^\prime = x + \Delta x $$

This formalism was always strange to me. I think a better way to implement the same equation is:

$$ J_f \Delta x = f(x) \quad x^\prime = x - \Delta x $$

By changing the definition/sign of delta x you avoid negating the RHS. While it ulitmetly doesn't matter too much where you perform the negation I think its clearer to not have the residual negated (and easier to integrate with models).

This is the formalism openvaf/osdi is build around.

In accordance with that formalism the limit_rhs needs to be negative (it wasn't before that was a bug). I am not willing to deviate from that formalism as I see Openvaf/OSDI closer to the model (which does not have a negative sign) than the simulator. I only really added the spice functions because reconstructing the spicelike RHS from a sparse matrix is quite involved/inefficient.

If you want to use the opposite formalism I would recommend just accumulating all contributions from all OSDI models in a single vector and then subtract that from your actual residual. Compared to model code/sparse solver a single such subtraction shouldn't matter much (and would vectorized nicely/have good cache locality).