OSDI Extension: Delays and Discontinuities

[pepijndevos]

I've been wrestling with how to best add absdelay to OpenVAF, which makes a Delay Differential Equation.
After studying implementations in SciML and Diffrax I nudged Claude towards the following proposal, which also covers discontinuities which are currently silently ignored. Ignorant simulators can remain ignorant, while good simulators can do better.

Motivation

Verilog-A's absdelay and $discontinuity are currently poorly supported in the
OpenVAF/OSDI ecosystem. absdelay is lowered to a passthrough (the delayed value is
silently dropped), and $discontinuity is ignored except for $discontinuity(-1) inside
$limit contexts. This causes incorrect simulation results, as demonstrated by
OpenVAF#106 where a 20 GHz signal
through a 10ps delay showed attenuation and distortion instead of a clean time shift.

This proposal adds two small extensions to OSDI that let compiled models declare their
delay and discontinuity structure, and optionally let capable simulators provide
high-quality delay evaluation.

Design Principles

• Minimal API surface. Two new fields on OsdiDescriptor, two on OsdiSimInfo.
• Graceful degradation. Unaware simulators ignore the new fields entirely. Models
  include a built-in fallback for delay evaluation.
• Clean mapping to SciML. The descriptors map directly to DDE problem construction
  (constant_lags, dependent_lags, tstops, ContinuousCallback).

Extension 1: Discontinuity Descriptors

Appended to OsdiDescriptor:

uint32_t num_discontinuities;
OsdiDiscontinuity *discontinuities;

Each descriptor declares one source of discontinuities or delays in the model:

#define OSDI_DISCONT_FIXED_TIME  0  /* known time point, optionally periodic */
#define OSDI_DISCONT_CROSSING    1  /* state-dependent zero crossing        */
#define OSDI_DISCONT_CONST_DELAY 2  /* constant delay (absdelay with param) */
#define OSDI_DISCONT_VAR_DELAY   3  /* state-dependent delay                */

typedef struct OsdiDiscontinuity {
    uint32_t type;
    double time;     /* time point or delay value tau             */
    double period;   /* repetition period, 0 if not periodic      */
    double (*fn)(void *inst, void *model, double *state, double time);
    int32_t order;   /* discontinuity order, 0 if unused (delays) */
} OsdiDiscontinuity;

Field usage by type:

Type	time	period	fn	order
FIXED_TIME	time point	repeat period (0 = once)	NULL	0, 1, 2, ...
CROSSING	unused (0)	unused (0)	zero-crossing function	0, 1, 2, ...
CONST_DELAY	tau	unused (0)	NULL	unused (0)
VAR_DELAY	unused (0)	unused (0)	lag function	unused (0)

The function pointer signature is the same for crossings and variable delays:

double (*fn)(void *inst, void *model, double *state, double time);

For CROSSING: returns a value; the simulator root-finds for zero crossings.
For VAR_DELAY: returns the current delay tau; the simulator uses it for
discontinuity tracking (detecting when t - tau(state, t) crosses a known
past discontinuity).

Extension 2: Past State Query

Appended to OsdiSimInfo:

double (*query_past_state)(void *ctx, uint32_t state_idx, double time);
void *query_past_state_ctx;

The model calls this during eval() to get node voltages and branch currents at
past times. state_idx uses the same index space as the existing node mapping.

If query_past_state is NULL, the model falls back to an internal ring buffer.

OpenVAF Implementation

Lowering absdelay

Currently absdelay(expr, tau, max_delay) is lowered to just expr (line 719 of
hir_lower/src/expr.rs). The new lowering:

1. MIR analysis identifies the states feeding into expr (node voltages, branch
   currents). This dependency analysis already exists in OpenVAF for Jacobian computation.

2. Descriptor emission. For constant tau (a parameter), emit OSDI_DISCONT_CONST_DELAY
   with time = tau. For state-dependent tau, emit OSDI_DISCONT_VAR_DELAY with a
   compiled function that evaluates the delay expression.

3. Eval codegen. For each absdelay call, emit code that:

   • Computes t_delayed = abstime - tau
   • For each state feeding into expr: queries it at t_delayed via
     query_past_state if available, otherwise reads from an internal ring buffer
   • Evaluates expr using the delayed states

4. Internal ring buffer (fallback). Only allocated if query_past_state is NULL at
   setup_instance time. Sized from max_delay (the third argument Verilog-A requires).
   Stores (time, value) pairs for each state feeding into delayed expressions. Uses
   cubic Hermite interpolation. Written to at every eval() call. If the simulator
   provides query_past_state, no buffer is allocated and no memory is wasted.

Lowering $discontinuity

Currently ignored (except -1 in $limit). The new lowering:

• $discontinuity(n) at a statically-known time → emit OSDI_DISCONT_FIXED_TIME with
  the time and order = n.
• $discontinuity(n) inside a conditional on state → emit OSDI_DISCONT_CROSSING with a
  compiled function evaluating the branch condition, order = n. This lets the simulator
  root-find proactively rather than react after the fact.

All discontinuity reporting is purely declarative via descriptors. There is no runtime
flag — the existing EVAL_RET_FLAG_LIM remains as-is for its original purpose (Newton
convergence hinting for $limit), but is unrelated to this proposal.

PWL Sources (periodic breakpoints)

For a PWL source with finite repeat count: expand breakpoints at compile time into
individual OSDI_DISCONT_FIXED_TIME entries. No periodicity needed.

For truly periodic sources (indefinite repetition): emit one OSDI_DISCONT_FIXED_TIME
per breakpoint within one period, each with period set to the repetition period.

Simulator Integration

Ignorant Simulator (e.g. ngspice, no changes)

Does nothing. Sets neither field on OsdiSimInfo. The model uses its internal ring
buffer for delays. Results are adequate when the simulator timestep is small enough
relative to the signal bandwidth, which can be encouraged by the model writing
bound_step (already supported by OSDI).

The model can set bound_step = tau / N for some reasonable N when delays are present,
ensuring the ring buffer has enough samples even without simulator awareness.

Minimally Aware Simulator

Reads num_discontinuities and the descriptors. Uses the information:

• FIXED_TIME entries → add to tstops / forced stepping times
• CONST_DELAY entries → set max timestep to min(tau) / N
• Callback types → ignored

Still sets query_past_state = NULL. Model uses internal buffer. But timestep control
ensures the internal buffer has good resolution.

Fully Capable Simulator (e.g. Cadnip.jl)

Reads all descriptors and maps them to SciML constructs:

for d in descriptors
    if d.type == OSDI_DISCONT_FIXED_TIME
        if d.period > 0
            # periodic: add tstops at d.time + k*d.period for k = 0, 1, ...
            # up to method order for delay-induced, or indefinitely for sources
        else
            push!(tstops, d.time)
        end
        push!(d_discontinuities, Discontinuity(d.time, d.order))

    elseif d.type == OSDI_DISCONT_CROSSING
        push!(continuous_callbacks, ContinuousCallback(
            (u, t, integrator) -> d.fn(inst, model, u, t),
            integrator -> nothing  # affect: just step onto it
        ))

    elseif d.type == OSDI_DISCONT_CONST_DELAY
        push!(constant_lags, d.time)

    elseif d.type == OSDI_DISCONT_VAR_DELAY
        push!(dependent_lags, (u, p, t) -> d.fn(inst, model, u, t))
    end
end

Provides query_past_state using the solver's dense output interpolant:

function query_past_state(ctx, state_idx, t_delayed)
    sol = unsafe_pointer_to_objref(ctx)
    return sol(t_delayed; idxs=state_idx)
end

For delays, constructs a DDEProblem with MethodOfSteps and gets proper
discontinuity tracking, high-order interpolation of delayed states, and
unconstrained stepping with fixed-point iteration — all from SciML's existing
DDE infrastructure.

Summary

	Descriptors on OsdiDescriptor	Callbacks on OsdiSimInfo
New fields	num_discontinuities, *discontinuities	query_past_state, query_past_state_ctx
Types	4 (fixed time, crossing, const delay, var delay)	1 function pointer

Total: one struct definition, four type constants, four new fields across two
existing structs.

[javaNoviceProgrammer]

Can you please implement absdelay on a separate branch? I can give it a try and test it after it's implemented.

[javaNoviceProgrammer]

Hello,
I spent an hour with Claude Code and it was able to modify the source code of both openvaf and ngspice (osdi interface) for the correct implementation of absdelay. Here's the summary as a markdown file:

absdelay_implementation.md

Please review and implement.
@arpadbuermen

[arpadbuermen]

@javaNoviceProgrammer I don't see the sources, just a .md file so I have nothing to review. But judging from the md file there are two major shortcomings.

1. absdelay() must also contribute to the Jacobian. The proposal does nothing about it. The contribution is nonzero if the point in the past at which we are computing the delayed output is between the last stored timepoint and the timepoint at which we are solving the circuit. Not common, but happens when simulator increases the timestep.
2. You did not consider how to handle absdelay in AC analysis.

My suggestion: each absdelay(y,td) creates an internal node nint. In time domain analysis simulator provides the output value of absdelay(y,td) and d(absdelay(y,td))/dy, i.e. derivative wrt value of y at point where we are solving the circuit. This derivative is 0 if td>=currest timestep. Then simulator uses this information to fill the nint's row of Jacobian and the RHS. In expression where we have f(z,absdelay(y,td)) is replaced with f(z,nint). If you don't do this the simulator will encounter slow convergence or convergence failures when delay is <=current timestep. A similar bug can be found in spice ideal transmission line where this derivative is ignored.

Introduction of an internal node makes it possible to handle absdelay in the frequency domain. Because absdelay(y,td) is Y exp(-j omega td) by chain rule we have g exp(-j omega td) where g=dy/dx and x is an input variable on which y depends. Since it does not make sense to introduce a third Jacobian (for absdelay() induced contributions) it would make sense to add (to the descriptor) a table of all input variables (nodes or node differences) on which y depends and the corresponding offsets into the instance data structure where the values of g are stored. We would also need the value of td at the operating point. It should be stored in the instance data structure and an offset should be provided in the descriptor. From stored g values and stored delay the simulator could then compute the contributions to the AC matrix.

Doing it only for tran is a half baked solution that would have to be torn down when one implements the whole thing for the AC analysis. Doing it correctly is much more work than indicated by the attached md file.

This has been written at 2300 hrs and not by Claude so please excuse the minor errors. Also I think Claude should be used with care. I already had to fix 3 bugfixes created by Claude (they were all incorrect). Up to now cleaning up after Claude turned out to be a total waste of time, unfortunately. It is good for reverse engineering and documentation, though.

An afterthought. Introducing 2 nodes per absdelay(y, td) would simplify things even more. One node would hold the the value of y, the other one the value of z=absdelay(y, td). Now the derivatives (g) simply become parts of the resistive Jacobian in row corresponding to y. They are generated by the existing derivatives machinery of OpenVAF. The AC analysis would only have to fill the extra equation (row corresponding to z) of the form

Y exp(-j omega td) - Z = 0

In time domain this equation would be

absdelay(y, td) - z = 0

Its Jacobian contribution

+-d(absdelay(y,td))/dy in the y's columns and -+1 in the z's columns. Note that y and z are input variables so they can be a node difference or a single ended node.

In this solution OSDI exposes the two nodes corresponding to y, the two nodes corresponding to z, and the row corresponding to z. It also exposes the offsets into the instance structure where the current value of td (td can be an expression) and the value of tdmax (I think there exists a variant of absdelay with tdmax, must read the LRM though) are held.

I would go with this 2 node approach.

To summarize. I am not going to implement things the way they are suggested in the .md file.

[javaNoviceProgrammer]

@javaNoviceProgrammer I don't see the sources, just a .md file so I have nothing to review. But judging from the md file there are two major shortcomings.

1. absdelay() must also contribute to the Jacobian. The proposal does nothing about it. The contribution is nonzero if the point in the past at which we are computing the delayed output is between the last stored timepoint and the timepoint at which we are solving the circuit. Not common, but happens when simulator increases the timestep.
2. You did not consider how to handle absdelay in AC analysis.

My suggestion: each absdelay(y,td) creates an internal node nint. In time domain analysis simulator provides the output value of absdelay(y,td) and d(absdelay(y,td))/dy, i.e. derivative wrt value of y at point where we are solving the circuit. This derivative is 0 if td>=currest timestep. Then simulator uses this information to fill the nint's row of Jacobian and the RHS. In expression where we have f(z,absdelay(y,td)) is replaced with f(z,nint). If you don't do this the simulator will encounter slow convergence or convergence failures when delay is <=current timestep. A similar bug can be found in spice ideal transmission line where this derivative is ignored.

Introduction of an internal node makes it possible to handle absdelay in the frequency domain. Because absdelay(y,td) is Y exp(-j omega td) by chain rule we have g exp(-j omega td) where g=dy/dx and x is an input variable on which y depends. Since it does not make sense to introduce a third Jacobian (for absdelay() induced contributions) it would make sense to add (to the descriptor) a table of all input variables (nodes or node differences) on which y depends and the corresponding offsets into the instance data structure where the values of g are stored. We would also need the value of td at the operating point. It should be stored in the instance data structure and an offset should be provided in the descriptor. From stored g values and stored delay the simulator could then compute the contributions to the AC matrix.

Doing it only for tran is a half baked solution that would have to be torn down when one implements the whole thing for the AC analysis. Doing it correctly is much more work than indicated by the attached md file.

This has been written at 2300 hrs and not by Claude so please excuse the minor errors. Also I think Claude should be used with care. I already had to fix 3 bugfixes created by Claude (they were all incorrect). Up to now cleaning up after Claude turned out to be a total waste of time, unfortunately. It is good for reverse engineering and documentation, though.

An afterthought. Introducing 2 nodes per absdelay(y, td) would simplify things even more. One node would hold the the value of y, the other one the value of z=absdelay(y, td). Now the derivatives (g) simply become parts of the resistive Jacobian in row corresponding to y. They are generated by the existing derivatives machinery of OpenVAF. The AC analysis would only have to fill the extra equation (row corresponding to z) of the form

Y exp(-j omega td) - Z = 0

In time domain this equation would be

absdelay(y, td) - z = 0

Its Jacobian contribution

+-d(absdelay(y,td))/dy in the y's columns and -+1 in the z's columns. Note that y and z are input variables so they can be a node difference or a single ended node.

In this solution OSDI exposes the two nodes corresponding to y, the two nodes corresponding to z, and the row corresponding to z. It also exposes the offsets into the instance structure where the current value of td (td can be an expression) and the value of tdmax (I think there exists a variant of absdelay with tdmax, must read the LRM though) are held.

I would go with this 2 node approach.

To summarize. I am not going to implement things the way they are suggested in the .md file.

@arpadbuermen : ok, thanks for the insightful comment. I will continue working with Claude to improve and check all possible simulations. I had checked DC and transient, but not AC. So, will focus on that.
Do you think the suggested implementation on the openvaf/Rust side is good enough?

As a side note, about a year ago I discovered that the osdi devices (N-devices) had a shortcoming in the netlist such that they could not handle more than 15 nodes in netlist (I think it was ngspice 43). I reported that to the ngspice team and they fixed it. In ngspice 46, the number of nodes of an N-device is fully dynamic and can be 50 or more. See link below:
https://sourceforge.net/p/ngspice/discussion/133842/thread/4145b7074d/?limit=25#fa7c

[arpadbuermen]

@javaNoviceProgrammer Sorry it was long. But it was a chain of thought written down at 11pm. :)

Your proposed approach will lead to convergence problems and make it very hard to correctly implement ac or any other frequency domain analysis (I dare not to say impossible, bc there is always somebody who proves you wrong by compiling Doom to run on an air conditioner). If your goal is transient analysis only, you should still implement the absdelay Jacobian contribution.

I still think you should switch to the 2 nodes per absdelay approach. It is simpler to implement and allows for ac analysis support. As soon as you will try to add the Jacobian contribution with your approach it will become complicated.

The mentioned shortcoming is a Ngspice shortcoming. It has nothing to do with OpenVAF. Good to hear they fixed it.

In my opinion Claude would do a much better job reverse engineering and documenting OpenVAF sources. This is something that is needed badly. Claude should be used to write code in a complex system like OpenVAF only if you know how to write the same code by hand.

[pepijndevos]

I have found SciML an invaluable resource of thorough and correct diffeq implementations, and delay equations seem no exception: https://docs.sciml.ai/DiffEqDocs/stable/tutorials/dde_example/

If you want to use Claude I would suggest to let it study their implementation and OpenVAF to plan a correct approach that considers the jumps in derivatives that the method of steps handles.

IIRC I made a start doing this and the correct implementation requires an API change. Let me dig it up... oh wait it's actually this issue lol

[javaNoviceProgrammer]

@arpadbuermen : your suggestion on two-node implementation was a good idea. My verison1 implementation was based on signal history storage which Claude eventually couldn't reconcile with AC sim. So, I switched to your suggestion and it seems to work correctly for all sims (DC, AC, Transient). One extra challenge was to modify both SPARSE and KLU solvers, which took extra time.

Here is my current git repo where I have the modified source code and documentation (Enhancement-1.md) and the binarires for apple-silicon. I plan to fully implement github CI for all-platform builds.

https://github.com/javaNoviceProgrammer/Ngspice_OpenVAF_Enhancements

https://github.com/javaNoviceProgrammer/Ngspice_OpenVAF_Enhancements/blob/main/Enhancement-1.md

@pepijndevos : I have run some test runs and benchmarks for DC, AC, Transient sims in this folder. It would be great if you can run more test cases and let me know if there is any bug:

https://github.com/javaNoviceProgrammer/Ngspice_OpenVAF_Enhancements/tree/main/absdelay_examples

My next step to implement "implicit contribution statements" which openvaf doesn't seem to currently support.

[arpadbuermen]

@javaNoviceProgrammer I have a few suggestions on how to organize the OSDI API for absdelay(). A few guidelines

• Avoid adding functions to the API unless necessary. Function calls waste time.
• A model is just a model. Keep analysis-specific computation in the simulator. Therefore value history and interpolation should be managed by the simulator.

My suggestion for the API is:

1. For each absdelay() the instance structure has one real slot for the current value of delay and one optional real slot for the current value of maxdelay. The descriptor holds the offsets of these slots in a table, one row for each absdelay(). If maxdelay is not present its offset is INT32_MAX. delay/absdelay are computed by eval().

2. For absdelay(V(a,b),...) or for absdelay(V(a), ...) the input node pair is (a, b) or (a, GND). The node pair local indices for each absdelay() are stored in a table in the descriptor. One row holds the local node indices of the pair corresponding to one absdelay(). Note that GND is indicated by UINT32_MAX. The simulator can compute the absdelay() input from the circuit's solution based on this information, instance's node bindings, and the current solution vector. For absdelay(some_expression, ...) an internal node x is created. The equation corresponding to that node is
   expression - V(x) = 0.
   OpenVAF is responsible for computing and loading the Jacobian and the residual of this equation. Since it contains nothing exotic (like absdelay which has been replaced by its output node value) this can be done by the existing mechanisms of OpenVAF. The input node pair reported in the descriptor is for this case (x, GND).

3. Each absdelay() creates an output node y. The output node's local index is in a table in the descriptor, one row per absdelay().

You can merge tables from 1) and 2) into a single table, one row per absdelay(), 2 columns for input, one column for output node.

Simulator's job

The equation in the time domain corresponding to node y is
delayed(V(x)) - V(y) = 0

The simulator is responsible for computing and loading the residual and the Jacobian for this equation. The history of x is managed by the simulator. The simulator can compute the new entry for the history from the input nodes table in the descriptor, the instance node bindings, and the current solution. The simulator is responsible for computing the delayed value of x based on the stored history and the computed delay/maxdelay obtained from the instance structure's slots. The Jacobian of this equation (transient analysis) is

row/column   x node 1           x node 2          y node
y node       d delayed(x)/dx   -d delayed(x)/dx   -1

d delayed(x)/dx is zero if td >= h where h is the current timestep (from last stored solution to the solution being computed now). Otherwise it depends on the interpolation algorithm and td.

The naive approach sets d delayed(x)/dx = 0 always.

In AC analysis the equation corresponding to node y is

exp(-j omega td) V(x) - V(y) = 0

Complex entries in the matrix row y are

row/column  x node 1           x node 2             y node
y node      exp(-j omega td)   -exp(-j omega td)    -1

This approach introduces strict separation between the model and the analyses. Simulator only gets information on absdelay input/output nodes, the value of delay, and the value of maxdelay. Everything else is computed in the simulator. This way also advanced analyses (PSS, HB) can be implemented without having to make any changes in OpenVAF. I know that for tran/ac the job is simple. But for anything advanced (RF analyses) this strict separation is a necessity.

If you have comments, let me know so we can polish this specification. I would appreciate if you could implement things this way. Once reviewed I will gladly merge the absdelay() support with OpenVAF-reloaded.

[arpadbuermen]

@javaNoviceProgrammer Just read your Enhancement-1.md file. The implementation is mostly in line with what I wrote above. If I understand correctly, no functions are added to the API.

I would like to expose max_delay (just like delay) in the instance data structure with an offset given in the descriptor (or UINT32_MAX when no maxdelay is given).

I would merge the extra data structures into the descriptor at its end. This way simulators ignorant of absdelay() will still work with models that do not use absdelay(). OSDI version can remain 0.4. We add

• uint32_t absdelay_count .. scalar integer count of absdelays in the module
• const OsdiAbsDelayInfo *absdelay_info .. array of entries for this module, one per absdelay
  to the descriptor and we are done.
  If we keep adding "partial descriptors" API will become a mess.

Also I would like to keep the special case when x can be expressed as a difference of two unknowns or equals an unknown. In this case we can drop the input internal node. The price is that x is now defined as the difference between two unknowns, i.e. nodes, and the absdelay descriptor holds two node indices for delay's input instead of one.

I am a bit baffled... why did you have to modify SPARSE/KLU?

[javaNoviceProgrammer]

@javaNoviceProgrammer Just read your Enhancement-1.md file. The implementation is mostly in line with what I wrote above. If I understand correctly, no functions are added to the API.

I would like to expose max_delay (just like delay) in the instance data structure with an offset given in the descriptor (or UINT32_MAX when no maxdelay is given).

I would merge the extra data structures into the descriptor at its end. This way simulators ignorant of absdelay() will still work with models that do not use absdelay(). OSDI version can remain 0.4. We add

• uint32_t absdelay_count .. scalar integer count of absdelays in the module
• const OsdiAbsDelayInfo *absdelay_info .. array of entries for this module, one per absdelay
  to the descriptor and we are done.
  If we keep adding "partial descriptors" API will become a mess.

Also I would like to keep the special case when x can be expressed as a difference of two unknowns or equals an unknown. In this case we can drop the input internal node. The price is that x is now defined as the difference between two unknowns, i.e. nodes, and the absdelay descriptor holds two node indices for delay's input instead of one.

I am a bit baffled... why did you have to modify SPARSE/KLU?

@arpadbuermen : ok, my friend. I leave it up to you how you'd like to implement/support absdelay() function in OpenVAF. I have a test folder in my github fork that you may use to double check your implementation for DC, AC, and Transient simulations.

https://github.com/javaNoviceProgrammer/Ngspice_OpenVAF_Enhancements/tree/main/absdelay_examples

I added GitHub CI to my fork repo and was able to set things up for binary builds of ngspice-46 and openvaf for "linux-intel", "linux-arm", "macos-apple-silicon" and "windows-intel" based on LLVM 18.

https://github.com/javaNoviceProgrammer/Ngspice_OpenVAF_Enhancements/tree/main/bin