feat(cllc): P4 — SPICE rewrite (active SR, snubbers, GEAR, IC, Evab, sense, analytical fallback)

Per CLLC_REWRITE_PLAN.md §6:
- Active SR bridge Sa..Sd replacing Ds1..Ds4 (synchronous gating tied
  to pwm1/pwm2) — eliminates B-source ideal-diode chatter that broke
  ngspice convergence in forward mode.
- Per-switch RC snubbers (Rsn_*=1k, Csn_*=1nF) across each switch.
- SW1 model VT=2.5 VH=0.8 RON=0.01 ROFF=1Meg (Guide §5).
- Evab VCVS probe `V(node_a)-V(node_b)` for converter-port stream.
- METHOD=GEAR + TRTOL=7 solver options.
- IC pre-charge: v(vout_p)=Vo, v(pri_l1_in)=0, v(sec_c2_in)=0.
- Vout_sense zero-volt source on the output rail.
- Time step period/200 (period/500 above 1.5·fr).
- Divide-by-zero guard: Rload = max(Vo/max(Io, 1e-9), 1e-3).
- Throws on switchingFrequency <= 0.
- Analytical fallback in simulate_and_extract_operating_points when
  ngspice unavailable or sim fails — tags OP name with [analytical]
  (DAB.cpp:1249–1299 pattern).

iPri NRMSE on PtP fixtures (gate 30%):
  Telecom-500W   4.42%
  Telecom-250W   3.56%
  Infineon-11kW 15.45%   (was 23.2% pre-SR rewrite)

Tests: +Test_Cllc_SPICE_Netlist, +Test_Cllc_DivideByZero_Guard,
updated Test_CllcConverter_NetlistGeneration (Sa..Sd checks).
26/26 [cllc-topology] tests passing.

Author: AlfVII

src/converter_models/Cllc.cpp

@@ -4,6 +4,7 @@
 #include "support/Utils.h"
 #include <cfloat>
 #include <cmath>
+#include <iostream>
 #include <complex>
 #include <limits>
 #include "support/Exceptions.h"
@@ -1179,6 +1180,10 @@ double get_value_or(T&& val, double default_val) {
         double inputVoltage = inputVoltages[inputVoltageIndex];
         auto opPoint = get_operating_points()[operatingPointIndex];
         double switchingFrequency = get_value_or(opPoint.get_switching_frequency(), 0.0);
+        if (switchingFrequency <= 0) {
+            throw std::invalid_argument(
+                "CLLC SPICE: operating point has invalid switching frequency");
+        }
         double outputVoltage = opPoint.get_output_voltages()[0];
         double outputCurrent = opPoint.get_output_currents()[0];
         double n = turnsRatio;
@@ -1188,6 +1193,7 @@ double get_value_or(T&& val, double default_val) {
         double L2 = params.secondaryResonantInductance;
         double C2 = params.secondaryResonantCapacitance;
         double Lm = params.magnetizingInductance;
+        double fr = params.resonantFrequency;
 
         double period = 1.0 / switchingFrequency;
         double halfPeriod = period / 2.0;
@@ -1206,115 +1212,152 @@ double get_value_or(T&& val, double default_val) {
         int totalPeriods = steadyPeriods + periodsToExtract;
         double simTime = totalPeriods * period;
         double startTime = steadyPeriods * period;
-        double stepTime = period / 200.0;
+        // Time step: period/200 baseline, period/500 when fs > 1.5·fr
+        // (super-resonant ringing demands tighter resolution).
+        // CLLC_REWRITE_PLAN §6.2.
+        double stepDivisor = (fr > 0 && switchingFrequency > 1.5 * fr) ? 500.0 : 200.0;
+        double stepTime = period / stepDivisor;
 
-        // Load resistance
-        double Rload = outputVoltage / outputCurrent;
+        // Load resistance — divide-by-zero guard per plan §6.2 (R_load floor).
+        double Rload = std::max(outputVoltage / std::max(outputCurrent, 1e-9), 1e-3);
 
         // Secondary inductance for coupled inductor model
         double Lsec = Lm / (n * n);  // Lm referred to secondary for coupling
 
         std::ostringstream circuit;
+        circuit << std::scientific;
 
         circuit << "* CLLC Bidirectional Resonant Converter - Generated by OpenMagnetics\n";
         circuit << "* Vin=" << inputVoltage << "V, Vout=" << outputVoltage
-                << "V, f=" << (switchingFrequency/1e3) << "kHz\n";
+                << "V, Pout=" << (outputVoltage*outputCurrent) << "W, fs="
+                << (switchingFrequency/1e3) << "kHz, fr=" << (fr/1e3) << "kHz\n";
         circuit << "* n=" << n << ", L1=" << (L1*1e6) << "uH, C1=" << (C1*1e9) << "nF\n";
         circuit << "* Lm=" << (Lm*1e6) << "uH, L2=" << (L2*1e6) << "uH, C2=" << (C2*1e9) << "nF\n\n";
 
         // DC Input
         circuit << "* DC Input\n";
         circuit << "Vin vin_p 0 " << inputVoltage << "\n\n";
 
-        // Switch models
-        circuit << "* Switch and diode models\n";
-        circuit << ".model SW1 SW VT=2.5 VH=0.5\n";
-        circuit << ".model DIDEAL D(IS=1e-14 RS=0.01 CJO=1e-12)\n\n";
+        // Switch model: VT=2.5 VH=0.8 RON=0.01 ROFF=1Meg per plan §6.2
+        circuit << "* Switch model (CLLC_REWRITE_PLAN §6.2)\n";
+        circuit << ".model SW1 SW(VT=2.5 VH=0.8 RON=0.01 ROFF=1Meg)\n\n";
 
         // PWM control signals for full bridge
         // Pair 1: S1 (high-side leg A), S4 (low-side leg B) - first half cycle
         // Pair 2: S2 (low-side leg A), S3 (high-side leg B) - second half cycle
         circuit << "* PWM control signals (complementary pairs with dead time)\n";
-        circuit << "Vpwm1 pwm1 0 PULSE(0 5 0 10n 10n " << std::scientific << tOn
-                << " " << period << std::fixed << ")\n";
-        circuit << "Vpwm2 pwm2 0 PULSE(0 5 " << std::scientific << halfPeriod
-                << " 10n 10n " << tOn << " " << period << std::fixed << ")\n\n";
-
-        // Primary full bridge
-        circuit << "* Primary Full Bridge\n";
-        circuit << "* Leg A: S1 (high-side), S2 (low-side)\n";
+        circuit << "Vpwm1 pwm1 0 PULSE(0 5 0 10n 10n " << tOn
+                << " " << period << ")\n";
+        circuit << "Vpwm2 pwm2 0 PULSE(0 5 " << halfPeriod
+                << " 10n 10n " << tOn << " " << period << ")\n\n";
+
+        // Primary full bridge — 4 active switches.
+        circuit << "* Primary Full Bridge (4 active switches)\n";
         circuit << "S1 vin_p node_a pwm1 0 SW1\n";
-        circuit << "S2 node_a 0 pwm2 0 SW1\n";
-        circuit << "* Leg B: S3 (high-side), S4 (low-side)\n";
+        circuit << "S2 node_a 0     pwm2 0 SW1\n";
         circuit << "S3 vin_p node_b pwm2 0 SW1\n";
-        circuit << "S4 node_b 0 pwm1 0 SW1\n\n";
+        circuit << "S4 node_b 0     pwm1 0 SW1\n";
+        // Per-switch snubbers (1k + 1nF) — 4 primary, 4 secondary = 8 total
+        circuit << "* Primary snubbers (1k + 1nF across each switch)\n";
+        circuit << "Rsn_S1 vin_p ns1 1k\n   Csn_S1 ns1 node_a 1n\n";
+        circuit << "Rsn_S2 node_a ns2 1k\n  Csn_S2 ns2 0 1n\n";
+        circuit << "Rsn_S3 vin_p ns3 1k\n   Csn_S3 ns3 node_b 1n\n";
+        circuit << "Rsn_S4 node_b ns4 1k\n  Csn_S4 ns4 0 1n\n\n";
+
+        // Bridge midpoint differential probe.
+        circuit << "* Bridge differential voltage probe\n";
+        circuit << "Evab vab 0 VALUE={V(node_a) - V(node_b)}\n\n";
 
         // Primary current sense
         circuit << "* Primary current sense\n";
         circuit << "Vpri_sense node_a pri_c1_in 0\n\n";
 
         // Primary resonant tank: C1 in series with L1
-        circuit << "* Primary Resonant Tank (C1 series with L1)\n";
-        circuit << "C_res1 pri_c1_in pri_l1_in " << std::scientific << C1 << std::fixed << "\n";
-        circuit << "L_res1 pri_l1_in pri_trafo_in " << std::scientific << L1 << std::fixed << "\n\n";
-
-        // Transformer: coupled inductors (Lpri = Lm, Lsec = Lm/n²)
-        // The magnetizing inductance IS the primary inductor of the transformer
-        circuit << "* Transformer (coupled inductors)\n";
-        circuit << "Lpri pri_trafo_in node_b " << std::scientific << Lm << std::fixed << "\n";
-        circuit << "Lsec sec_trafo_p sec_trafo_n " << std::scientific << Lsec << std::fixed << "\n";
+        circuit << "* Primary Resonant Tank (Cr1 series with Lr1)\n";
+        circuit << "C_res1 pri_c1_in pri_l1_in " << C1 << "\n";
+        circuit << "L_res1 pri_l1_in pri_trafo_in " << L1 << "\n\n";
+
+        // Transformer: coupled inductors (Lpri = Lm, Lsec = Lm/n²). K = 0.9999
+        // (never 1.0) so SPICE doesn't choke on a singular coupling matrix.
+        circuit << "* Transformer (coupled inductors, K=0.9999 — never 1.0)\n";
+        circuit << "Lpri pri_trafo_in node_b " << Lm << "\n";
+        circuit << "Lsec sec_trafo_p sec_trafo_n " << Lsec << "\n";
         circuit << "Kpri_sec Lpri Lsec 0.9999\n\n";
 
         // Secondary resonant tank: L2 in series with C2
-        circuit << "* Secondary Resonant Tank (L2 series with C2)\n";
+        circuit << "* Secondary Resonant Tank (Lr2 series with Cr2) + sense\n";
         circuit << "Vsec_sense sec_trafo_p sec_l2_in 0\n";
-        circuit << "L_res2 sec_l2_in sec_c2_in " << std::scientific << L2 << std::fixed << "\n";
-        circuit << "C_res2 sec_c2_in node_c " << std::scientific << C2 << std::fixed << "\n\n";
+        circuit << "L_res2 sec_l2_in sec_c2_in " << L2 << "\n";
+        circuit << "C_res2 sec_c2_in node_c " << C2 << "\n\n";
 
-        // Secondary node_d = sec_trafo_n
-        circuit << "* Secondary bridge reference\n";
+        // node_d = sec_trafo_n; high-Z DC path for solver convergence
+        circuit << "* Secondary bridge reference node\n";
         circuit << "Vd_ref sec_trafo_n node_d 0\n";
-        // Add high-value resistor to provide DC path for convergence
         circuit << "Rdc_sec sec_trafo_n 0 1G\n\n";
 
-        // Secondary full bridge rectifier (diodes for simplicity)
-        // In forward mode: Ds1,Ds4 conduct first half; Ds2,Ds3 conduct second half
-        circuit << "* Secondary Full Bridge Rectifier (diodes)\n";
-        circuit << "Ds1 node_c vout_p DIDEAL\n";
-        circuit << "Ds2 vout_n node_c DIDEAL\n";
-        circuit << "Ds3 node_d vout_p DIDEAL\n";
-        circuit << "Ds4 vout_n node_d DIDEAL\n\n";
-
-        // Snubber resistors for convergence
-        circuit << "* Snubber resistors for convergence\n";
-        circuit << "Rsnub1 node_c vout_p 1MEG\n";
-        circuit << "Rsnub2 vout_n node_c 1MEG\n";
-        circuit << "Rsnub3 node_d vout_p 1MEG\n";
-        circuit << "Rsnub4 vout_n node_d 1MEG\n\n";
+        // ---- Active synchronous-rectifier full bridge on secondary ----
+        // Per plan §6.1: in forward mode the SR is gated synchronously with
+        // the primary PWM (pwm_s1 = pwm_p1, pwm_s2 = pwm_p2). This avoids
+        // the chatter that an ideal-diode B-source emulator suffers at the
+        // V(D,S) ≈ 0 transition (timestep collapse). Reverse mode (P8)
+        // will drive the secondary as inverter and re-purpose the primary
+        // as the synchronous rectifier with its own gating.
+        //
+        // Conduction map (forward, full-bridge):
+        //   pwm1 high  → S1 (high-A) + S4 (low-B) ON  → node_a > node_b
+        //              → secondary node_c > node_d
+        //              → SR Sa (node_c → vout_p) and Sd (vout_n → node_d) ON
+        //   pwm2 high  → opposite polarity  → SR Sb (vout_n → node_c) and
+        //                                       Sc (node_d → vout_p) ON
+        circuit << "* Secondary Active Synchronous Rectifier (4 SR switches)\n";
+        circuit << "* pwm1 -> Sa, Sd    pwm2 -> Sb, Sc   (forward mode SR)\n";
+        circuit << "Sa node_c vout_p pwm1 0 SW1\n";
+        circuit << "Sb vout_n node_c pwm2 0 SW1\n";
+        circuit << "Sc node_d vout_p pwm2 0 SW1\n";
+        circuit << "Sd vout_n node_d pwm1 0 SW1\n";
+        // Per-switch snubbers (1k + 1nF) on the SR side
+        circuit << "* Secondary snubbers (1k + 1nF across each SR switch)\n";
+        circuit << "Rsn_Sa node_c nsa 1k\n  Csn_Sa nsa vout_p 1n\n";
+        circuit << "Rsn_Sb vout_n nsb 1k\n  Csn_Sb nsb node_c 1n\n";
+        circuit << "Rsn_Sc node_d nsc 1k\n  Csn_Sc nsc vout_p 1n\n";
+        circuit << "Rsn_Sd vout_n nsd 1k\n  Csn_Sd nsd node_d 1n\n\n";
 
         // Ground reference for secondary
         circuit << "* Secondary ground reference\n";
         circuit << "Vgnd_sec vout_n 0 0\n\n";
 
-        // Output filter and load
-        circuit << "* Output filter and load\n";
+        // Output filter, sense source, and load
+        circuit << "* Output filter, current sense, and load\n";
         circuit << "Cout vout_p vout_n 100u IC=" << outputVoltage << "\n";
-        circuit << "Rload vout_p vout_n " << Rload << "\n\n";
-
-        // Transient analysis
-        circuit << "* Transient Analysis\n";
-        circuit << ".tran " << std::scientific << stepTime << " " << simTime
-                << " " << startTime << std::fixed << " UIC\n\n";
-
-        // Save signals
-        circuit << "* Output signals\n";
-        circuit << ".save v(vin_p) v(pri_trafo_in) v(node_b) i(Vin) i(Vpri_sense)"
-                << " v(sec_trafo_p) v(sec_trafo_n) i(Vsec_sense)"
-                << " v(vout_p) v(vout_n)\n\n";
-
-        // Options
-        circuit << ".options RELTOL=0.003 ABSTOL=1e-8 VNTOL=1e-5 TRTOL=10 ITL1=500 ITL4=100\n";
-        circuit << ".ic v(vout_p)=" << outputVoltage << "\n\n";
+        circuit << "Vout_sense vout_p vout_load 0\n";
+        circuit << "Rload vout_load vout_n " << std::fixed << Rload
+                << std::scientific << "\n\n";
+
+        // Transient analysis with UIC honoring the .ic statements below.
+        circuit << "* Transient Analysis (UIC honors .ic pre-charge)\n";
+        circuit << ".tran " << stepTime << " " << simTime
+                << " " << startTime << " UIC\n\n";
+
+        // Save signals — node-level and current branches that the
+        // extractor maps into Primary/Secondary windings and ports.
+        circuit << "* Save signals\n";
+        circuit << ".save v(vab) v(vin_p) v(pri_trafo_in) v(node_a) v(node_b)\n";
+        circuit << "+ v(node_c) v(node_d) v(sec_trafo_p) v(sec_trafo_n)\n";
+        circuit << "+ v(vout_p) v(vout_n) v(pri_c1_in) v(sec_c2_in)\n";
+        circuit << "+ i(Vin) i(Vpri_sense) i(Vsec_sense) i(Vout_sense)\n\n";
+
+        // Solver options — verbatim from plan §6.2
+        circuit << "* Solver options (CLLC_REWRITE_PLAN §6.2)\n";
+        circuit << ".options RELTOL=0.01 ABSTOL=1e-7 VNTOL=1e-4 ITL1=500 ITL4=500\n";
+        circuit << ".options METHOD=GEAR TRTOL=7\n\n";
+
+        // Initial conditions: pre-charge output cap to Vo, leave resonant
+        // caps at zero. SPICE node names for cap voltages: the node BETWEEN
+        // C1 (pri_c1_in -- pri_l1_in) and BETWEEN C2 (sec_c2_in -- node_c).
+        circuit << "* Initial conditions (UIC pre-charge)\n";
+        circuit << ".ic v(vout_p)=" << std::fixed << outputVoltage
+                << std::scientific << "\n";
+        circuit << ".ic v(pri_l1_in)=0 v(sec_c2_in)=0\n\n";
 
         circuit << ".end\n";
 
@@ -1333,7 +1376,9 @@ double get_value_or(T&& val, double default_val) {
 
         NgspiceRunner runner;
         if (!runner.is_available()) {
-            throw std::runtime_error("ngspice is not available for simulation");
+            std::cerr << "[CLLC] ngspice not available — falling back to analytical "
+                         "operating-point generation" << std::endl;
+            return process_operating_points(turnsRatios, magnetizingInductance);
         }
 
         std::vector<double> inputVoltages;
@@ -1370,7 +1415,19 @@ double get_value_or(T&& val, double default_val) {
                 auto simResult = runner.run_simulation(netlist, config);
 
                 if (!simResult.success) {
-                    throw std::runtime_error("CLLC Simulation failed: " + simResult.errorMessage);
+                    std::cerr << "[CLLC] ngspice run failed: " << simResult.errorMessage
+                              << " — falling back to analytical for this op." << std::endl;
+                    auto analytical = process_operating_point_for_input_voltage(
+                        inputVoltages[inputVoltageIndex], cllcOpPoint, n,
+                        magnetizingInductance, params);
+                    std::string aname = inputVoltagesNames[inputVoltageIndex] + " input volt.";
+                    if (get_operating_points().size() > 1) {
+                        aname += " op. point " + std::to_string(opIndex);
+                    }
+                    aname += " [analytical]";
+                    analytical.set_name(aname);
+                    operatingPoints.push_back(analytical);
+                    continue;
                 }
 
                 // Map waveform names to winding excitations