spice21/
analysis.rs

1//!
2//! # Spice21 Analyses
3//!
4use num::{Complex, Float, Zero};
5use serde::{Deserialize, Serialize};
6use std::collections::HashMap;
7use std::ops::Index;
8
9use crate::circuit::{Ckt, NodeRef};
10use crate::comps::{Component, ComponentSolver};
11use crate::defs;
12use crate::sparse21::{Eindex, Matrix};
13use crate::{sperror, SpNum, SpResult};
14
15///
16/// # Matrix Stamps
17///
18/// `Stamps` are the interface between Components and Solvers.
19/// Each Component returns `Stamps` from each call to `load`,
20/// conveying its Matrix-contributions in `Stamps.g`
21/// and its RHS contributions in `Stamps.b`.
22#[derive(Debug)]
23pub(crate) struct Stamps<NumT> {
24    pub(crate) g: Vec<(Option<Eindex>, NumT)>,
25    pub(crate) b: Vec<(Option<VarIndex>, NumT)>,
26}
27impl<NumT: SpNum> Stamps<NumT> {
28    pub fn new() -> Stamps<NumT> {
29        Stamps { g: vec![], b: vec![] }
30    }
31}
32
33#[derive(Debug, Clone, Copy, Deserialize, Serialize)]
34pub(crate) enum VarKind {
35    V = 0,
36    I,
37    Q,
38}
39
40#[derive(Clone, Copy, Debug, Deserialize, Serialize)]
41pub struct VarIndex(pub usize);
42
43#[derive(Debug, Deserialize, Serialize)]
44pub(crate) struct Variables<NumT> {
45    kinds: Vec<VarKind>,
46    values: Vec<NumT>,
47    names: Vec<String>,
48}
49impl<NumT: SpNum> Variables<NumT> {
50    pub fn new() -> Self {
51        Variables {
52            kinds: vec![],
53            values: vec![],
54            names: vec![],
55        }
56    }
57    /// Convert Variables<OtherT> to Variables<NumT>
58    /// Keeps all `kinds`, while resetting all values to zero.
59    fn from<OtherT>(other: Variables<OtherT>) -> Self {
60        Variables {
61            kinds: other.kinds,
62            names: other.names,
63            values: vec![NumT::zero(); other.values.len()],
64        }
65    }
66    /// Add a new Variable with attributes `name` and `kind`.
67    pub fn add(&mut self, name: String, kind: VarKind) -> VarIndex {
68        // FIXME: check if present
69        self.kinds.push(kind);
70        self.names.push(name);
71        self.values.push(NumT::zero());
72        return VarIndex(self.kinds.len() - 1);
73    }
74    /// Add a new Voltage Variable with name `name`.
75    pub fn addv(&mut self, name: String) -> VarIndex {
76        self.add(name, VarKind::V)
77    }
78    /// Add a new Current Variable with name `name`.
79    pub fn addi(&mut self, name: String) -> VarIndex {
80        self.add(name, VarKind::I)
81    }
82    /// Find a variable named `name`. Returns `VarIndex` if found, `None` if not present.
83    pub fn find<S: Into<String>>(&self, name: S) -> Option<VarIndex> {
84        let n = name.into().clone();
85        match self.names.iter().position(|x| *x == n) {
86            Some(i) => Some(VarIndex(i)),
87            None => None,
88        }
89    }
90    /// Retrieve the Variable corresponding to Node `node`,
91    /// creating it if necessary.
92    pub fn find_or_create(&mut self, node: NodeRef) -> Option<VarIndex> {
93        match node {
94            NodeRef::Gnd => None,
95            NodeRef::Name(name) => {
96                // FIXME: shouldn't have to clone all the names here
97                match self.names.iter().cloned().position(|x| x == name) {
98                    Some(i) => Some(VarIndex(i)),
99                    None => Some(self.add(name.clone(), VarKind::V)),
100                }
101            }
102            NodeRef::Num(num) => {
103                let name = num.to_string();
104                // FIXME: shouldn't have to clone all the names here
105                match self.names.iter().cloned().position(|x| x == name) {
106                    Some(i) => Some(VarIndex(i)),
107                    None => Some(self.add(name.clone(), VarKind::V)),
108                }
109            }
110        }
111    }
112    /// Retrieve a Variable value.
113    /// `None` represents "ground" and always has value zero.
114    pub fn get(&self, i: Option<VarIndex>) -> NumT {
115        match i {
116            None => NumT::zero(),
117            Some(ii) => self.values[ii.0],
118        }
119    }
120    pub fn len(&self) -> usize {
121        self.kinds.len()
122    }
123}
124
125/// Solver Iteration Struct
126/// Largely for debug of convergence and progress
127#[allow(dead_code)] // Used for debug
128struct Iteration<NumT: SpNum> {
129    n: usize,
130    x: Vec<NumT>,
131    dx: Vec<NumT>,
132    vtol: Vec<bool>,
133    itol: Vec<bool>,
134}
135
136/// Newton-Style Iterative Solver
137/// Owns each of its circuit's ComponentSolvers,
138/// its SparseMatrix, and Variables.
139pub(crate) struct Solver<'a, NumT: SpNum> {
140    pub(crate) comps: Vec<ComponentSolver<'a>>,
141    pub(crate) vars: Variables<NumT>,
142    pub(crate) mat: Matrix<NumT>,
143    pub(crate) rhs: Vec<NumT>,
144    pub(crate) history: Vec<Vec<NumT>>,
145    pub(crate) defs: defs::Defs,
146    pub(crate) opts: Options,
147}
148
149/// Real-valued Solver specifics
150/// FIXME: nearly all of this *should* eventually be share-able with the Complex Solver
151impl Solver<'_, f64> {
152    /// Collect and incorporate updates from all components
153    fn update(&mut self, an: &AnalysisInfo) {
154        for comp in self.comps.iter_mut() {
155            let updates = comp.load(&self.vars, an, &self.opts);
156            // Make updates for G and b
157            for upd in updates.g.iter() {
158                if let (Some(ei), val) = *upd {
159                    self.mat.update(ei, val);
160                }
161            }
162            for upd in updates.b.iter() {
163                if let (Some(ei), val) = *upd {
164                    self.rhs[ei.0] += val;
165                }
166            }
167        }
168    }
169    fn solve(&mut self, an: &AnalysisInfo) -> SpResult<Vec<f64>> {
170        self.history = vec![]; // Reset our guess-history
171        let mut dx = vec![0.0; self.vars.len()];
172
173        for _k in 0..100 {
174            // FIXME: number of iterations
175            // Make a copy of state for tracking
176            self.history.push(self.vars.values.clone());
177            // Reset our matrix and RHS vector
178            self.mat.reset();
179            self.rhs = vec![0.0; self.vars.len()];
180
181            // Load up component updates
182            self.update(an);
183
184            // Calculate the residual error
185            let res: Vec<f64> = self.mat.res(&self.vars.values, &self.rhs)?;
186
187            // Check convergence
188            if self.converged(&dx, &res) {
189                // Converged. Commit component states
190                for c in self.comps.iter_mut() {
191                    c.commit();
192                }
193                return Ok(self.vars.values.clone()); // FIXME: stop cloning
194            }
195            // Haven't Converged. Solve for our update.
196            dx = self.mat.solve(res)?;
197            let max_step = 1000e-3;
198            let max_abs = dx.iter().fold(0.0, |s, v| if v.abs() > s { v.abs() } else { s });
199            if max_abs > max_step {
200                for r in 0..dx.len() {
201                    dx[r] = dx[r] * max_step / max_abs;
202                }
203            }
204            // And update our guess
205            for r in 0..self.vars.len() {
206                self.vars.values[r] += dx[r];
207            }
208        }
209        return Err(sperror("Convergence Failed"));
210    }
211}
212
213/// Complex-Valued Solver Specifics
214/// FIXME: nearly all of this *should* eventually be share-able with the Real Solver
215impl<'a> Solver<'a, Complex<f64>> {
216    /// Create a Complex solver from a real-valued one.
217    /// Commonly deployed when moving from DCOP to AC analysis.
218    fn from(re: Solver<'a, f64>) -> Self {
219        let mut op = Solver::<'a, Complex<f64>> {
220            comps: re.comps,
221            vars: Variables::<Complex<f64>>::from(re.vars),
222            mat: Matrix::new(),
223            rhs: vec![],
224            history: vec![],
225            defs: re.defs,
226            opts: re.opts,
227        };
228
229        // Create matrix elements, over-writing each Component's pointers
230        for comp in op.comps.iter_mut() {
231            comp.create_matrix_elems(&mut op.mat);
232        }
233        return op;
234    }
235
236    /// Collect and incorporate updates from all components
237    fn update(&mut self, an: &AnalysisInfo) {
238        for comp in self.comps.iter_mut() {
239            let updates = comp.load_ac(&self.vars, an, &self.opts);
240            // Make updates for G and b
241            for upd in updates.g.iter() {
242                if let (Some(ei), val) = *upd {
243                    self.mat.update(ei, val);
244                }
245            }
246            for upd in updates.b.iter() {
247                if let (Some(ei), val) = *upd {
248                    self.rhs[ei.0] += val;
249                }
250            }
251        }
252    }
253    fn solve(&mut self, an: &AnalysisInfo) -> SpResult<Vec<Complex<f64>>> {
254        self.history = vec![]; // Reset our guess-history
255        let mut dx = vec![Complex::zero(); self.vars.len()];
256        let mut iters: Vec<Iteration<Complex<f64>>> = vec![];
257
258        for _k in 0..20 {
259            // FIXME: number of iterations
260            // Make a copy of state for tracking
261            self.history.push(self.vars.values.clone());
262            // Reset our matrix and RHS vector
263            self.mat.reset();
264            self.rhs = vec![Complex::zero(); self.vars.len()];
265
266            // Load up component updates
267            self.update(an);
268
269            // Calculate the residual error
270            let res: Vec<Complex<f64>> = self.mat.res(&self.vars.values, &self.rhs)?;
271            let vtol: Vec<bool> = dx.iter().map(|&v| v.norm() < 1e-3).collect();
272            let itol: Vec<bool> = res.iter().map(|&v| v.norm() < 1e-9).collect();
273            // Check convergence
274            if vtol.iter().all(|v| *v) && itol.iter().all(|v| *v) {
275                // Commit component results
276                for c in self.comps.iter_mut() {
277                    c.commit();
278                }
279                return Ok(self.vars.values.clone());
280            }
281            // Solve for our update
282            dx = self.mat.solve(res)?;
283            let max_step = 1.0;
284            let max_abs = dx.iter().fold(0.0, |s, v| if v.norm() > s { v.norm() } else { s });
285            if max_abs > max_step {
286                for r in 0..dx.len() {
287                    dx[r] = dx[r] * max_step / max_abs;
288                }
289            }
290            // And update our guess
291            for r in 0..self.vars.len() {
292                self.vars.values[r] += dx[r];
293            }
294            iters.push(Iteration {
295                n: _k,
296                x: self.vars.values.clone(),
297                dx: dx.clone(),
298                vtol: vtol,
299                itol: itol,
300            });
301        }
302        return Err(sperror("Convergence Failed"));
303    }
304}
305
306impl<'a, NumT: SpNum> Solver<'a, NumT> {
307    /// Create a new Solver, translate `Ckt` Components into its `ComponentSolvers`.
308    pub(crate) fn new(ckt: Ckt, opts: Options) -> Solver<'a, NumT> {
309        // Elaborate the circuit
310        use crate::elab::{elaborate, Elaborator};
311        let e = elaborate(ckt, opts);
312        let Elaborator {
313            defs, mut comps, vars, opts, ..
314        } = e;
315        // Create our matrix and its elements
316        let mut mat = Matrix::new();
317        for comp in comps.iter_mut() {
318            comp.create_matrix_elems(&mut mat);
319        }
320        // And return a Solver with the combination
321        Solver {
322            comps,
323            vars,
324            mat,
325            rhs: Vec::new(),
326            history: Vec::new(),
327            defs,
328            opts,
329        }
330    }
331    fn converged(&self, dx: &Vec<NumT>, res: &Vec<NumT>) -> bool {
332        // Inter-step Newton convergence
333        for e in dx.iter() {
334            if e.absv() > self.opts.reltol {
335                return false;
336            }
337        }
338        // KCL convergence
339        for e in res.iter() {
340            if e.absv() > self.opts.iabstol {
341                return false;
342            }
343        }
344        return true;
345    }
346}
347
348/// Operating Point Result
349#[derive(Debug)]
350pub struct OpResult {
351    pub names: Vec<String>,
352    pub values: Vec<f64>,
353    pub map: HashMap<String, f64>,
354}
355impl OpResult {
356    /// Create an OpResult from a (typically final) set of `Variables`.
357    fn from(vars: Variables<f64>) -> Self {
358        let mut map: HashMap<String, f64> = HashMap::new();
359        for i in 0..vars.names.len() {
360            map.insert(vars.names[i].clone(), vars.values[i]);
361        }
362        let Variables { names, values, .. } = vars;
363        OpResult { names, values, map }
364    }
365    /// Get the value of signal `signame`, or an `SpError` if not present
366    pub(crate) fn get<S: Into<String>>(&self, signame: S) -> SpResult<f64> {
367        match self.map.get(&signame.into()) {
368            Some(v) => Ok(v.clone()),
369            None => Err(sperror("Signal Not Found")),
370        }
371    }
372}
373/// Maintain much (most?) of our original vector-result-format
374/// via enabling integer indexing
375impl Index<usize> for OpResult {
376    type Output = f64;
377    fn index(&self, t: usize) -> &f64 {
378        &self.values[t]
379    }
380}
381
382/// Dc Operating Point Analysis
383pub fn dcop(ckt: Ckt, opts: Option<Options>) -> SpResult<OpResult> {
384    let o = if let Some(o) = opts { o } else { Options::default() };
385    let mut s = Solver::<f64>::new(ckt, o);
386    let _r = s.solve(&AnalysisInfo::OP)?;
387    return Ok(OpResult::from(s.vars));
388}
389pub(crate) enum AnalysisInfo<'a> {
390    OP,
391    TRAN(&'a TranOptions, &'a TranState),
392    AC(&'a AcOptions, &'a AcState),
393}
394
395pub(crate) enum NumericalIntegration {
396    BE,
397    TRAP,
398}
399
400impl Default for NumericalIntegration {
401    fn default() -> NumericalIntegration {
402        NumericalIntegration::BE
403    }
404}
405
406/// # TranState
407///
408/// Internal state of transient analysis
409/// Often passed to Components for timesteps etc
410#[derive(Default)]
411pub(crate) struct TranState {
412    pub(crate) t: f64,
413    pub(crate) dt: f64,
414    pub(crate) vic: Vec<usize>,
415    pub(crate) ric: Vec<usize>,
416    pub(crate) ni: NumericalIntegration,
417}
418impl TranState {
419    /// Numerical Integration
420    pub fn integrate(&self, dq: f64, dq_dv: f64, vguess: f64, _ip: f64) -> (f64, f64, f64) {
421        let dt = self.dt;
422        match self.ni {
423            NumericalIntegration::BE => {
424                let g = dq_dv / dt;
425                let i = dq / dt;
426                let rhs = i - g * vguess;
427                (g, i, rhs)
428            }
429            NumericalIntegration::TRAP => {
430                panic!("Trapezoidal Integration is Disabled!");
431                // let g = 2.0 * dq_dv / dt;
432                // let i = 2.0 * dq / dt - ip;
433                // let rhs = i - g * vguess;
434                // (g, i, rhs)
435            }
436        }
437    }
438    pub fn integq(&self, dq: f64, dq_dv: f64, vguess: f64, _ip: f64) -> ChargeInteg {
439        let (g, i, rhs) = self.integrate(dq, dq_dv, vguess, _ip);
440        ChargeInteg { g, i, rhs }
441    }
442}
443/// Result of numerical integration for a charge-element
444#[derive(Default, Clone)]
445pub(crate) struct ChargeInteg {
446    pub(crate) g: f64,
447    pub(crate) i: f64,
448    pub(crate) rhs: f64,
449}
450/// Transient Analysis Options
451#[derive(Debug)]
452pub struct TranOptions {
453    pub tstep: f64,
454    pub tstop: f64,
455    pub ic: Vec<(NodeRef, f64)>,
456}
457impl TranOptions {
458    pub fn decode(bytes_: &[u8]) -> SpResult<Self> {
459        // FIXME: remove
460        use prost::Message;
461        use std::io::Cursor;
462
463        // Decode the protobuf version
464        let proto = proto::TranOptions::decode(&mut Cursor::new(bytes_))?;
465        // And convert into the real thing
466        Ok(Self::from(proto))
467    }
468}
469impl From<proto::TranOptions> for TranOptions {
470    fn from(i: proto::TranOptions) -> Self {
471        use super::circuit::n;
472        let mut ic: Vec<(NodeRef, f64)> = vec![];
473        for (name, val) in &i.ic {
474            ic.push((n(name), val.clone()));
475        }
476        Self {
477            tstep: i.tstep,
478            tstop: i.tstop,
479            ic,
480        }
481    }
482}
483impl Default for TranOptions {
484    fn default() -> TranOptions {
485        Self::from(proto::TranOptions::default())
486    }
487}
488
489pub(crate) struct Tran<'a> {
490    solver: Solver<'a, f64>,
491    state: TranState,
492    pub(crate) opts: TranOptions,
493}
494
495impl<'a> Tran<'a> {
496    pub fn new(ckt: Ckt, opts: Options, args: TranOptions) -> Tran<'a> {
497        let solver = Solver::new(ckt, opts);
498        let ics = args.ic.clone();
499        let mut t = Tran {
500            solver,
501            opts: args,
502            state: TranState::default(),
503        };
504        for (node, val) in &ics {
505            t.ic(node.clone(), *val);
506        }
507        t
508    }
509    /// Create and set an initial condition on Node `n`, value `val`.
510    pub fn ic(&mut self, n: NodeRef, val: f64) {
511        use crate::comps::{Resistor, Vsrc};
512
513        // Create two new variables: the forcing voltage, and current in its source
514        let fnode = self.solver.vars.add(format!(".{}.vic", n.to_string()), VarKind::V);
515        let ivar = self.solver.vars.add(format!(".{}.iic", n.to_string()), VarKind::I);
516
517        let mut r = Resistor::new(1.0, Some(fnode), self.solver.vars.find_or_create(n)); // FIXME: rforce value
518        r.create_matrix_elems(&mut self.solver.mat);
519        self.solver.comps.push(r.into());
520        self.state.ric.push(self.solver.comps.len() - 1);
521        let mut v = Vsrc::new(val, 0.0, Some(fnode), None, ivar);
522        v.create_matrix_elems(&mut self.solver.mat);
523        self.solver.comps.push(v.into());
524        self.state.vic.push(self.solver.comps.len() - 1);
525    }
526    pub fn solve(&mut self) -> SpResult<TranResult> {
527        // Initialize results
528        let mut results = TranResult::new();
529        results.signals(&self.solver.vars);
530
531        // Solve for our initial condition
532        let tsoln = self.solver.solve(&AnalysisInfo::OP);
533        let tdata = match tsoln {
534            Ok(x) => x,
535            Err(e) => {
536                println!("Failed to find initial solution");
537                return Err(e);
538            }
539        };
540        results.push(self.state.t, &tdata);
541
542        // Update initial-condition sources and resistances
543        // FIXME: whether to change the voltages
544        for c in self.state.vic.iter() {
545            self.solver.comps[*c].update(0.0);
546        }
547        for c in self.state.ric.iter() {
548            self.solver.comps[*c].update(1e-9);
549        }
550
551        let mut tpoint: usize = 0;
552        let max_tpoints: usize = 1e9 as usize;
553        self.state.t = self.opts.tstep;
554        self.state.dt = self.opts.tstep;
555        while self.state.t < self.opts.tstop && tpoint < max_tpoints {
556            let aninfo = AnalysisInfo::TRAN(&self.opts, &self.state);
557            let tsoln = self.solver.solve(&aninfo);
558            let tdata = match tsoln {
559                Ok(x) => x,
560                Err(e) => {
561                    println!("Failed at t={}", self.state.t);
562                    return Err(e);
563                }
564            };
565            results.push(self.state.t, &tdata);
566
567            // self.state.ni = NumericalIntegration::TRAP; // FIXME!
568            tpoint += 1;
569            self.state.t += self.opts.tstep;
570        }
571        results.end();
572        Ok(results)
573    }
574}
575/// # TranResult
576/// In-Memory Store for transient data
577#[derive(Default, Serialize, Deserialize)]
578pub struct TranResult {
579    pub signals: Vec<String>,
580    pub time: Vec<f64>,
581    pub data: Vec<Vec<f64>>,
582    pub map: HashMap<String, Vec<f64>>,
583}
584impl TranResult {
585    pub fn new() -> Self {
586        Self {
587            signals: vec![],
588            time: vec![],
589            data: vec![],
590            map: HashMap::new(),
591        }
592    }
593    fn signals(&mut self, vars: &Variables<f64>) {
594        for name in vars.names.iter() {
595            self.signals.push(name.to_string());
596        }
597    }
598    fn push(&mut self, t: f64, vals: &Vec<f64>) {
599        self.time.push(t);
600        self.data.push(vals.clone());
601        // FIXME: filter out un-saved and internal variables
602    }
603    /// Simulation complete, re-org data into hash-map of signals
604    fn end(&mut self) {
605        self.map.insert("time".to_string(), self.time.clone());
606        for i in 0..self.signals.len() {
607            let mut vals: Vec<f64> = vec![];
608            for v in 0..self.time.len() {
609                vals.push(self.data[v][i]);
610            }
611            self.map.insert(self.signals[i].clone(), vals);
612        }
613    }
614    pub fn len(&self) -> usize {
615        self.time.len()
616    }
617    /// Retrieve values of signal `name`
618    pub fn get(&self, name: &str) -> SpResult<&Vec<f64>> {
619        match self.map.get(name) {
620            Some(v) => Ok(v),
621            None => Err(sperror(format!("Signal Not Found: {}", name))),
622        }
623    }
624}
625/// Maintain much (most?) of our original vector-result-format
626/// via enabling integer indexing
627impl Index<usize> for TranResult {
628    type Output = Vec<f64>;
629    fn index(&self, t: usize) -> &Vec<f64> {
630        &self.data[t]
631    }
632}
633
634/// Transient Analysis
635pub fn tran(ckt: Ckt, opts: Option<Options>, args: Option<TranOptions>) -> SpResult<TranResult> {
636    let o = if let Some(val) = opts { val } else { Options::default() };
637    let a = if let Some(val) = args { val } else { TranOptions::default() };
638    return Tran::new(ckt, o, a).solve();
639}
640
641/// Simulation Options
642pub struct Options {
643    pub temp: f64,
644    pub tnom: f64,
645    pub gmin: f64,
646    pub iabstol: f64,
647    pub reltol: f64,
648    pub chgtol: f64,
649    pub volt_tol: f64,
650    pub trtol: usize,
651    pub tran_max_iter: usize,
652    pub dc_max_iter: usize,
653    pub dc_trcv_max_iter: usize,
654    pub integrate_method: usize,
655    pub order: usize,
656    pub max_order: usize,
657    pub pivot_abs_tol: f64,
658    pub pivot_rel_tol: f64,
659    pub src_factor: f64,
660    pub diag_gmin: f64,
661}
662
663use crate::proto;
664
665impl From<proto::SimOptions> for Options {
666    fn from(i: proto::SimOptions) -> Self {
667        Self {
668            temp: if let Some(val) = i.temp { val } else { 300.15 },
669            tnom: if let Some(val) = i.tnom { val } else { 300.15 },
670            gmin: if let Some(val) = i.gmin { val } else { 1e-12 },
671            iabstol: if let Some(val) = i.iabstol { val } else { 1e-12 },
672            reltol: if let Some(val) = i.reltol { val } else { 1e-3 },
673            chgtol: 1e-14,
674            volt_tol: 1e-6,
675            trtol: 7,
676            tran_max_iter: 10,
677            dc_max_iter: 100,
678            dc_trcv_max_iter: 50,
679            integrate_method: 0,
680            order: 1,
681            max_order: 2,
682            pivot_abs_tol: 1e-13,
683            pivot_rel_tol: 1e-3,
684            src_factor: 1.0,
685            diag_gmin: 0.0,
686        }
687    }
688}
689impl Default for Options {
690    fn default() -> Self {
691        Self::from(proto::SimOptions::default())
692    }
693}
694
695#[derive(Default)]
696pub(crate) struct AcState {
697    pub omega: f64,
698}
699
700/// AC Analysis Options
701pub struct AcOptions {
702    pub fstart: usize,
703    pub fstop: usize,
704    pub npts: usize, // Total, not "per decade"
705}
706impl From<proto::AcOptions> for AcOptions {
707    fn from(i: proto::AcOptions) -> Self {
708        Self {
709            fstart: i.fstart as usize,
710            fstop: i.fstop as usize,
711            npts: i.npts as usize,
712        }
713    }
714}
715impl Default for AcOptions {
716    fn default() -> Self {
717        Self::from(proto::AcOptions::default())
718    }
719}
720
721/// AcResult
722/// In-Memory Store for Complex-Valued AC Data
723#[derive(Default, Serialize, Deserialize)]
724pub struct AcResult {
725    pub signals: Vec<String>,
726    pub freq: Vec<f64>,
727    pub data: Vec<Vec<Complex<f64>>>,
728    pub map: HashMap<String, Vec<Complex<f64>>>,
729}
730impl AcResult {
731    fn new() -> Self {
732        Self::default()
733    }
734    fn signals<T>(&mut self, vars: &Variables<T>) {
735        for name in vars.names.iter() {
736            self.signals.push(name.to_string());
737        }
738    }
739    fn push(&mut self, f: f64, vals: &Vec<Complex<f64>>) {
740        self.freq.push(f);
741        self.data.push(vals.clone());
742        // FIXME: filter out un-saved and internal variables
743    }
744    /// Simulation complete, re-org data into hash-map of signals
745    fn end(&mut self) {
746        // self.map.insert("freq".to_string(), self.freq.clone()); // FIXME: will require multi datatypes per result
747        for i in 0..self.signals.len() {
748            let mut vals: Vec<Complex<f64>> = vec![];
749            for v in 0..self.freq.len() {
750                vals.push(self.data[v][i]);
751            }
752            self.map.insert(self.signals[i].clone(), vals);
753        }
754    }
755    pub fn len(&self) -> usize {
756        self.freq.len()
757    }
758}
759
760/// AC Analysis
761pub fn ac(ckt: Ckt, opts: Option<Options>, args: Option<AcOptions>) -> SpResult<AcResult> {
762    /// FIXME: result saving is in flux, and essentially on three tracks:
763    /// * The in-memory format used by unit-tests returns vectors of complex numbers
764    /// * The first on-disk format, streaming JSON, falls down for nested data. It has complex numbers flattened, along with frequency.
765    /// * The AcResult struct holds all relevant data (in memory), and can serialize it to JSON (or any other serde format).
766    ///
767    use serde::ser::{SerializeSeq, Serializer};
768    use std::fs::File;
769
770    let opts = if let Some(val) = opts { val } else { Options::default() };
771    let args = if let Some(val) = args { val } else { AcOptions::default() };
772
773    // Initial DCOP solver and solution
774    let mut solver = Solver::<f64>::new(ckt, opts);
775    let _dc_soln = solver.solve(&AnalysisInfo::OP)?;
776
777    // Convert to an AC solver
778    let mut solver = Solver::<Complex<f64>>::from(solver);
779    let mut state = AcState::default();
780    let mut soln = vec![];
781
782    // Set up streaming writer
783    let rf = File::create("stream.ac.json").unwrap(); // FIXME: name
784    let mut ser = serde_json::Serializer::new(rf);
785    let mut seq = ser.serialize_seq(None).unwrap();
786
787    // Initialize results
788    let mut results = AcResult::new();
789    results.signals(&solver.vars);
790
791    // Set up frequency sweep
792    let mut f = args.fstart as f64;
793    let fstop = args.fstop as f64;
794    let fstep = (10.0).powf(f64::log10(fstop / f) / args.npts as f64);
795
796    // Main Frequency Loop
797    while f <= fstop {
798        use std::f64::consts::PI;
799        state.omega = 2.0 * PI * f;
800        let an = AnalysisInfo::AC(&args, &state);
801        let fsoln = solver.solve(&an)?;
802
803        // Push to our in-mem data
804        results.push(f, &fsoln);
805        // AND push to the flattened, streaming data
806        let mut flat: Vec<f64> = vec![f];
807        for pt in fsoln.iter() {
808            flat.push(pt.re);
809            flat.push(pt.im);
810        }
811        seq.serialize_element(&flat).unwrap();
812        // AND push to our simple vector-data
813        soln.push(fsoln);
814        // Last-iteration handling
815        if f == fstop {
816            break;
817        }
818        f = f64::min(f * fstep, fstop);
819    }
820    // Close up streaming results
821    SerializeSeq::end(seq).unwrap();
822
823    // Write the full-batch JSON result
824    use std::io::prelude::*;
825    let mut rfj = File::create("ac.json").unwrap(); // FIXME: name
826    let s = serde_json::to_string(&results).unwrap();
827    rfj.write_all(s.as_bytes()).unwrap();
828
829    // And return our results
830    results.end();
831    return Ok(results);
832}
spice21/analysis.rs

spice21/
analysis.rs
