feat(core): mutual inductance K, ic= ICs, .tran uic + Chua & Lin showcase

examples/showcase/main.tsx

@@ -205,6 +205,38 @@ D1 sw out DMOD
 C1 out 0 100u
 Rload out 0 10
 .tran 50n 500u`,
+  },
+  {
+    id: 'chua-lin',
+    name: 'Chua & Lin 8-7',
+    desc: 'Coupled inductors + UIC',
+    icon: '\u29c9',
+    group: 'Benchmarks',
+    tag: '.tran',
+    signals: ['x', 'y', 'out'],
+    tranNetlist: `
+* Chua & Lin "Computer-Aided Analysis" 8-7 (page 343).
+* Three magnetically-coupled inductors (K1/K2/K3) with two capacitor and
+* one inductor initial conditions. A classic SPICE benchmark \u2014 every
+* simulator gives a slightly different answer because it forces the
+* engine to handle three things at once: the K-element coupling matrix,
+* the .tran "uic" flag, and a stiff multi-time-scale transient.
+V1 inp 0 DC 1
+R10 inp x 0.5
+C2 x 0 2
+L8 x c 2H ic=2
+C12 x y 5 ic=2
+L9 y c 2H
+C3 y 0 4 ic=5
+R11 y out 0.25
+R5 out 0 2
+I1 0 out DC 1
+L6 c n1 4H
+R4 n1 0 1
+K1 L6 L8 -0.3535534
+K2 L6 L9 -0.3535534
+K3 L8 L9 -0.5
+.tran 0.2 200 0 0.2 uic`,
   },
   {
     integrationMethod: 'gear2',

packages/core/src/analysis/transient-driver.ts

@@ -104,6 +104,13 @@ interface InternalTransientConfig {
   maxTimestep: number;
   /** Optional pre-computed DC solution. When provided, skips internal DC op point. */
   initialSolution?: Float64Array;
+  /**
+   * The seed is a UIC initial-condition solution rather than a DC operating
+   * point. The companion model has no valid history at t=0, so the first
+   * step must bootstrap with Backward Euler and a cut dt — same pattern the
+   * driver uses after crossing a waveform breakpoint.
+   */
+  seedIsUIC?: boolean;
 }
 
 class TransientSimImpl implements TransientSim {
@@ -136,7 +143,15 @@ class TransientSimImpl implements TransientSim {
     if (config.initialSolution) {
       // Caller already computed DC — skip internal DC and seed directly.
       this.assembler.solution.set(config.initialSolution);
-      this.stampPrevB();
+      if (config.seedIsUIC) {
+        // Force BE bootstrap on the first step: leave prevB undefined and
+        // trigger the post-breakpoint dt cut. The seeded UIC solution does
+        // not satisfy any prior companion equation, so any trap/gear-2
+        // history term we built here would be mathematically inconsistent.
+        this.justCrossedBreakpoint = true;
+      } else {
+        this.stampPrevB();
+      }
     } else {
       this.initDC();
     }
@@ -369,13 +384,15 @@ export function createDriverFromCompiled(
     timestep: number;
     maxTimestep: number;
     initialSolution?: Float64Array;
+    seedIsUIC?: boolean;
   },
 ): TransientSim & { peekInitialStep(): TransientStep } {
   const impl = new TransientSimImpl(compiled, options, {
     stopTime: config.stopTime,
     timestep: config.timestep,
     maxTimestep: config.maxTimestep,
     initialSolution: config.initialSolution,
+    seedIsUIC: config.seedIsUIC,
   });
   return impl;
 }

packages/core/src/analysis/transient.ts

@@ -12,13 +12,15 @@ export function solveTransient(
   analysis: TransientAnalysis,
   options: ResolvedOptions,
   initialSolution?: Float64Array,
+  seedIsUIC?: boolean,
 ): TransientResult {
   const { nodeNames, branchNames } = compiled;
   const driver = createDriverFromCompiled(compiled, options, {
     stopTime: analysis.stopTime,
     timestep: analysis.timestep,
     maxTimestep: analysis.maxTimestep ?? Math.min(analysis.timestep, analysis.stopTime / 50),
     initialSolution,
+    seedIsUIC,
   });
 
   const timePoints: number[] = [0];

packages/core/src/analysis/uic.ts

@@ -0,0 +1,73 @@
+import type { Circuit, CompiledCircuit } from '../circuit.js';
+import type { ResolvedOptions } from '../types.js';
+import { solveDCOperatingPoint } from './dc.js';
+import { Inductor } from '../devices/inductor.js';
+
+/**
+ * Compute the initial solution vector for a `.tran ... uic` run.
+ *
+ * Substitutes capacitors with `ic=` for DC voltage sources, inductors with
+ * `ic=` for DC current sources, and drops mutual-inductance K-elements
+ * (which contribute nothing at DC). Solves the resulting linear DC system,
+ * then maps node voltages and branch currents back to the *original*
+ * compiled circuit's solution layout — which uses different branch indices
+ * because inductors-with-ic become branchless current sources, and
+ * capacitors-with-ic gain new branches.
+ *
+ * Storage elements without `ic` keep their default DC behaviour: caps act
+ * as opens, inductors as shorts. This matches the ngspice convention where
+ * the unspecified element settles to the value implied by the surrounding
+ * topology rather than being forced to zero.
+ *
+ * @returns Initial solution sized for the original `compiled` system.
+ */
+export function computeUICInitialSolution(
+  circuit: Circuit,
+  compiled: CompiledCircuit,
+  options: ResolvedOptions,
+): Float64Array {
+  const icCompiled = circuit.compileForUIC();
+  const { assembler: icAsm } = solveDCOperatingPoint(icCompiled, options);
+
+  const seed = new Float64Array(compiled.nodeCount + compiled.branchCount);
+
+  // Map node voltages by name. The IC system has the same node set as the
+  // original (substitutions only swap device types, not topology).
+  for (const [name, origIdx] of compiled.nodeIndexMap) {
+    if (origIdx < 0) continue;
+    const icIdx = icCompiled.nodeIndexMap.get(name);
+    if (icIdx !== undefined && icIdx >= 0) {
+      seed[origIdx] = icAsm.solution[icIdx];
+    }
+  }
+
+  // Map branch currents. Original branches come from V sources, inductors
+  // (all of them), and controlled sources (E/H). After substitution:
+  //   - V sources keep their branches (currents transfer directly).
+  //   - Inductors WITH ic became I sources (no branch in IC system); use the
+  //     declared ic as the seed current, with the same dot/sign convention
+  //     as the inductor's own KCL stamp (branch current enters n+).
+  //   - Inductors WITHOUT ic stayed as inductors (DC short, branch present).
+  //   - E/H controlled sources keep their branches.
+  for (let i = 0; i < compiled.branchNames.length; i++) {
+    const name = compiled.branchNames[i];
+    const origBranchAbs = compiled.nodeCount + i;
+
+    // Look up the original device to know what kind of element this branch
+    // belongs to.
+    const dev = compiled.devices.find(d => d.name === name);
+
+    if (dev instanceof Inductor && dev.ic !== undefined) {
+      seed[origBranchAbs] = dev.ic;
+      continue;
+    }
+
+    // Otherwise, find the same name in the IC system's branch layout.
+    const icBranchIdx = icCompiled.branchNames.indexOf(name);
+    if (icBranchIdx >= 0) {
+      seed[origBranchAbs] = icAsm.solution[icCompiled.nodeCount + icBranchIdx];
+    }
+  }
+
+  return seed;
+}

packages/core/src/circuit.ts

@@ -7,6 +7,7 @@ import { VoltageSource } from './devices/voltage-source.js';
 import { CurrentSource } from './devices/current-source.js';
 import { Capacitor } from './devices/capacitor.js';
 import { Inductor } from './devices/inductor.js';
+import { MutualInductor } from './devices/mutual-inductor.js';
 import { Diode } from './devices/diode.js';
 import { BJT } from './devices/bjt.js';
 import { MOSFET } from './devices/mosfet.js';
@@ -60,6 +61,11 @@ interface DeviceDescriptor {
   modelName?: string;
   params?: Record<string, number>;
   controlSource?: string;
+  /** Initial condition (volts for caps, amps for inductors). */
+  ic?: number;
+  /** For 'K' descriptors: names of the two inductors being coupled. */
+  coupledA?: string;
+  coupledB?: string;
 }
 
 /**
@@ -130,10 +136,10 @@ export class Circuit {
    * @param nodeNeg - Negative terminal node
    * @param capacitance - Capacitance value in farads
    */
-  addCapacitor(name: string, nodePos: string, nodeNeg: string, capacitance: number): void {
+  addCapacitor(name: string, nodePos: string, nodeNeg: string, capacitance: number, ic?: number): void {
     this.nodeSet.add(nodePos);
     this.nodeSet.add(nodeNeg);
-    this.descriptors.push({ type: 'C', name, nodes: [nodePos, nodeNeg], value: capacitance });
+    this.descriptors.push({ type: 'C', name, nodes: [nodePos, nodeNeg], value: capacitance, ic });
   }
 
   /**
@@ -144,10 +150,30 @@ export class Circuit {
    * @param nodeNeg - Negative terminal node
    * @param inductance - Inductance value in henries
    */
-  addInductor(name: string, nodePos: string, nodeNeg: string, inductance: number): void {
+  addInductor(name: string, nodePos: string, nodeNeg: string, inductance: number, ic?: number): void {
     this.nodeSet.add(nodePos);
     this.nodeSet.add(nodeNeg);
-    this.descriptors.push({ type: 'L', name, nodes: [nodePos, nodeNeg], value: inductance });
+    this.descriptors.push({ type: 'L', name, nodes: [nodePos, nodeNeg], value: inductance, ic });
+  }
+
+  /**
+   * Couple two inductors via a coupling coefficient k (K-element).
+   *
+   * The mutual inductance is M = k * sqrt(L_a * L_b), with |k| ≤ 1 for a
+   * physically realizable coupling. Both inductors must already exist via
+   * {@link addInductor}; their dot convention is implicit in their declared
+   * node order.
+   *
+   * @param name - K-element name (e.g., `'K1'`)
+   * @param indA - First inductor name (e.g., `'L1'`)
+   * @param indB - Second inductor name (e.g., `'L2'`)
+   * @param k - Coupling coefficient (positive or negative)
+   */
+  addInductorCoupling(name: string, indA: string, indB: string, k: number): void {
+    this.descriptors.push({
+      type: 'K', name, nodes: [],
+      coupledA: indA, coupledB: indB, value: k,
+    });
   }
 
   /**
@@ -361,7 +387,7 @@ export class Circuit {
    */
   addAnalysis(type: 'op'): void;
   addAnalysis(type: 'dc', params: { source: string; start: number; stop: number; step: number }): void;
-  addAnalysis(type: 'tran', params: { timestep: number; stopTime: number; startTime?: number; maxTimestep?: number }): void;
+  addAnalysis(type: 'tran', params: { timestep: number; stopTime: number; startTime?: number; maxTimestep?: number; uic?: boolean }): void;
   addAnalysis(type: 'ac', params: { variation: 'dec' | 'oct' | 'lin'; points: number; startFreq: number; stopFreq: number }): void;
   addAnalysis(type: string, params?: Record<string, unknown>): void {
     switch (type) {
@@ -378,13 +404,14 @@ export class Circuit {
         });
         break;
       case 'tran': {
-        const tranCmd: { type: 'tran'; timestep: number; stopTime: number; startTime?: number; maxTimestep?: number } = {
+        const tranCmd: { type: 'tran'; timestep: number; stopTime: number; startTime?: number; maxTimestep?: number; uic?: boolean } = {
           type: 'tran',
           timestep: params!.timestep as number,
           stopTime: params!.stopTime as number,
         };
         if (params?.startTime !== undefined) tranCmd.startTime = params.startTime as number;
         if (params?.maxTimestep !== undefined) tranCmd.maxTimestep = params.maxTimestep as number;
+        if (params?.uic) tranCmd.uic = true;
         this._analyses.push(tranCmd);
         break;
       }
@@ -450,9 +477,47 @@ export class Circuit {
    * @throws Error if a referenced subcircuit or control source is undefined
    * @throws {@link CycleError} if subcircuit instances form a circular dependency
    */
+  /**
+   * Compile a UIC seed circuit: capacitors with `ic=` are replaced by DC
+   * voltage sources of that value, and inductors with `ic=` are replaced
+   * by DC current sources. Storage elements without `ic` keep their default
+   * DC behaviour (caps open, inductors short). K-elements are dropped (they
+   * have no DC contribution). Used to compute the t=0 seed for `.tran ... uic`.
+   */
+  compileForUIC(): CompiledCircuit {
+    return this.compileWith(d => this.toUICDescriptor(d));
+  }
+
+  private toUICDescriptor(desc: DeviceDescriptor): DeviceDescriptor | null {
+    if (desc.type === 'C' && desc.ic !== undefined) {
+      return {
+        type: 'V', name: desc.name, nodes: desc.nodes,
+        waveform: { type: 'dc', value: desc.ic },
+      };
+    }
+    if (desc.type === 'L' && desc.ic !== undefined) {
+      // Inductor convention: positive branch current flows from n+ to n-,
+      // i.e., LEAVES n+. SPICE current source convention: positive value
+      // INJECTS into n+. Flip the node order so an I source of +ic models
+      // an inductor with branch current +ic.
+      return {
+        type: 'I', name: desc.name, nodes: [desc.nodes[1], desc.nodes[0]],
+        waveform: { type: 'dc', value: desc.ic },
+      };
+    }
+    if (desc.type === 'K') return null; // no DC effect
+    return desc;
+  }
+
   compile(): CompiledCircuit {
+    return this.compileWith(d => d);
+  }
+
+  private compileWith(transform: (d: DeviceDescriptor) => DeviceDescriptor | null): CompiledCircuit {
     // Pre-expand subcircuit instances into flat device descriptors
-    const expandedDescriptors = this.expandAllSubcircuits();
+    const expandedDescriptors = this.expandAllSubcircuits()
+      .map(transform)
+      .filter((d): d is DeviceDescriptor => d !== null);
 
     // Collect all nodes from expanded descriptors
     for (const desc of expandedDescriptors) {
@@ -503,12 +568,24 @@ export class Circuit {
           devices.push(new CurrentSource(desc.name, nodeIndices, resolveWaveform(desc.waveform)));
           break;
         case 'C':
-          devices.push(new Capacitor(desc.name, nodeIndices, desc.value!));
+          devices.push(new Capacitor(desc.name, nodeIndices, desc.value!, desc.ic));
           break;
         case 'L': {
           const bi = branchIndex++;
           branchNames.push(desc.name);
-          devices.push(new Inductor(desc.name, nodeIndices, bi, desc.value!));
+          devices.push(new Inductor(desc.name, nodeIndices, bi, desc.value!, desc.ic));
+          break;
+        }
+        case 'K': {
+          // Resolve the two coupled inductor names to their device instances.
+          const a = deviceMap.get(desc.coupledA!);
+          const b = deviceMap.get(desc.coupledB!);
+          if (!(a instanceof Inductor) || !(b instanceof Inductor)) {
+            throw new Error(
+              `K-element '${desc.name}' references unknown or non-inductor device(s): ${desc.coupledA}, ${desc.coupledB}`,
+            );
+          }
+          devices.push(new MutualInductor(desc.name, a, b, desc.value!));
           break;
         }
         case 'D': {

packages/core/src/devices/capacitor.ts

@@ -8,6 +8,7 @@ export class Capacitor implements DeviceModel {
     readonly name: string,
     readonly nodes: number[],
     public capacitance: number,
+    public ic?: number,
   ) {}
 
   setParameter(value: number): void {

packages/core/src/devices/index.ts

@@ -2,6 +2,7 @@ export type { DeviceModel, StampContext } from './device.js';
 export { Resistor } from './resistor.js';
 export { Capacitor } from './capacitor.js';
 export { Inductor } from './inductor.js';
+export { MutualInductor } from './mutual-inductor.js';
 export { VoltageSource } from './voltage-source.js';
 export { CurrentSource } from './current-source.js';
 export { Diode } from './diode.js';

packages/core/src/devices/inductor.ts

@@ -9,6 +9,7 @@ export class Inductor implements DeviceModel {
     readonly nodes: number[],
     readonly branchIndex: number,
     public inductance: number,
+    public ic?: number,
   ) {
     this.branches = [branchIndex];
   }