From 73e259bb56cb120c9dbe204b334c26ae855d102c Mon Sep 17 00:00:00 2001 From: Holger Vogt Date: Sun, 26 Jul 2026 12:50:01 +0200 Subject: [PATCH] pa-134 VS --- src/frontend/Makefile.am | 2 + src/frontend/com_commands.h | 1 + src/frontend/com_hb.c | 77 +++++++++ src/frontend/com_hb.h | 7 + src/frontend/commands.c | 4 + src/include/ngspice/cktdefs.h | 1 + src/spicelib/analysis/dcpss.c | 300 ++++++++++++++++++++++++++++++++++ visualc/vngspice.vcxproj | 2 + 8 files changed, 394 insertions(+) create mode 100644 src/frontend/com_hb.c create mode 100644 src/frontend/com_hb.h diff --git a/src/frontend/Makefile.am b/src/frontend/Makefile.am index 9ba8b561c..8dd94a014 100644 --- a/src/frontend/Makefile.am +++ b/src/frontend/Makefile.am @@ -55,6 +55,8 @@ libfte_la_SOURCES = \ com_optimize.h \ com_qpss.c \ com_qpss.h \ + com_hb.c \ + com_hb.h \ com_checkpoint.c \ com_checkpoint.h \ com_option.c \ diff --git a/src/frontend/com_commands.h b/src/frontend/com_commands.h index 0c1b4e7c6..0fcbfb23b 100644 --- a/src/frontend/com_commands.h +++ b/src/frontend/com_commands.h @@ -8,6 +8,7 @@ void com_altermod(wordlist *wl); void com_alterparam(wordlist *wl); void com_optimize(wordlist *wl); /* Enhancement-130 */ void com_qpss(wordlist *wl); /* Enhancement-133 */ +void com_hb(wordlist *wl); /* Enhancement-134 */ void com_savestate(wordlist *wl); /* Enhancement-131 */ void com_loadstate(wordlist *wl); /* Enhancement-131 */ void com_meas(wordlist *wl); diff --git a/src/frontend/com_hb.c b/src/frontend/com_hb.c new file mode 100644 index 000000000..a539c3121 --- /dev/null +++ b/src/frontend/com_hb.c @@ -0,0 +1,77 @@ +/********** +Enhancement-134: Harmonic Balance -- `hb [points] [maxiter]`. + +Single-tone harmonic balance: find the periodic steady state in the FREQUENCY +domain by Newton, instead of integrating in time. Each node voltage is a truncated +Fourier series V(t)=sum_{k=-K..K} V_k e^{jk w0 t}; the KCL residual at each +node/harmonic is driven to zero with the E-121 conversion matrix as the Jacobian. +The heavy lifting is in HBanalyze() (spicelib/analysis/dcpss.c, which reuses the +conversion matrix + dense complex solver); this command parses the arguments, +makes sure the circuit is built, and runs it. +**********/ + +#include "ngspice/ngspice.h" +#include "ngspice/cpdefs.h" +#include "ngspice/cktdefs.h" +#include "ngspice/ftedefs.h" +#include "ngspice/fteext.h" +#include "ngspice/wordlist.h" +#include "ngspice/cpextern.h" + +#include "circuits.h" +#include "com_hb.h" + +static double hbnum(const char *w) +{ + char *s = (char *) w; + double v = 0.0; + if (ft_numparse(&s, FALSE, &v) < 0) + v = atof(w); + return v; +} + +void +com_hb(wordlist *wl) +{ + CKTcircuit *ckt; + double f0, tol = 1e-10; + int K, P = 0, maxiter = 50, verbose, err; + + if (!ft_curckt || !ft_curckt->ci_ckt) { + fprintf(cp_err, "Error: hb: there is no circuit loaded.\n"); + return; + } + ckt = ft_curckt->ci_ckt; + + if (!wl || !wl->wl_next) { + fprintf(cp_err, "Usage: hb [points] [maxiter]\n"); + return; + } + f0 = hbnum(wl->wl_word); + K = (int) hbnum(wl->wl_next->wl_word); + if (wl->wl_next->wl_next) { + P = (int) hbnum(wl->wl_next->wl_next->wl_word); + if (wl->wl_next->wl_next->wl_next) + maxiter = (int) hbnum(wl->wl_next->wl_next->wl_next->wl_word); + } + if (f0 <= 0.0 || K < 1) { + fprintf(cp_err, "Error: hb: need f0 > 0 and K >= 1.\n"); + return; + } + + /* make sure the circuit is built (matrix + states allocated) */ + if (ckt->CKTmatrix == NULL || SMPmatSize(ckt->CKTmatrix) <= 0) { + if ((err = CKTsetup(ckt)) != OK || (err = CKTtemp(ckt)) != OK) { + fprintf(cp_err, "Error: hb: circuit setup failed.\n"); + return; + } + } + + verbose = cp_getvar("hb_verbose", CP_BOOL, NULL, 0); + ft_curckt->ci_curTask = ft_curckt->ci_defTask; + ckt->CKTcurJob = ft_curckt->ci_defTask ? ft_curckt->ci_defTask->jobs : NULL; + + err = HBanalyze(ckt, f0, K, P, maxiter, tol, verbose ? 1 : 0); + if (err != OK) + fprintf(cp_err, "hb: harmonic balance did not complete (error %d).\n", err); +} diff --git a/src/frontend/com_hb.h b/src/frontend/com_hb.h new file mode 100644 index 000000000..0e75bb8b9 --- /dev/null +++ b/src/frontend/com_hb.h @@ -0,0 +1,7 @@ +#ifndef ngspice_COM_HB_H +#define ngspice_COM_HB_H + +/* Enhancement-134: harmonic balance. */ +void com_hb(wordlist *wl); + +#endif diff --git a/src/frontend/commands.c b/src/frontend/commands.c index 3673a379c..3fabde317 100644 --- a/src/frontend/commands.c +++ b/src/frontend/commands.c @@ -432,6 +432,10 @@ struct comm spcp_coms[] = { { 040, 040, 040, 040 }, E_DEFHMASK, 3, LOTS, NULL, "expr f1 f2 [periods] [maxorder] : two-tone quasi-periodic steady-state spectrum (intermodulation)." }, + { "hb", com_hb, TRUE, FALSE, /* Enhancement-134 */ + { 040, 040, 040, 040 }, E_DEFHMASK, 2, LOTS, + NULL, + "f0 K [points] [maxiter] : harmonic-balance steady-state spectrum (K harmonics)." }, { "savestate", com_savestate, FALSE, TRUE, /* Enhancement-131 */ { 1, 040000, 040000, 040000 }, E_DEFHMASK, 1, 1, NULL, diff --git a/src/include/ngspice/cktdefs.h b/src/include/ngspice/cktdefs.h index 687852e25..433a94275 100644 --- a/src/include/ngspice/cktdefs.h +++ b/src/include/ngspice/cktdefs.h @@ -460,6 +460,7 @@ extern int PSSaskQuest(CKTcircuit *, JOB *, int , IFvalue *); extern int PSSsetParm(CKTcircuit *, JOB *, int , IFvalue *); extern int PSSinit(CKTcircuit *, JOB *); extern int DCpss(CKTcircuit *, int); +extern int HBanalyze(CKTcircuit *, double f0, int K, int P, int maxiter, double tol, int verbose); /* E-134 */ #endif #ifdef RFSPICE diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c index 8aa98fd4d..a55f0f3d6 100644 --- a/src/spicelib/analysis/dcpss.c +++ b/src/spicelib/analysis/dcpss.c @@ -1227,6 +1227,306 @@ psp_sweep(CKTcircuit *ckt, PSSan *job) #endif +/* ===================== Enhancement-134: Harmonic Balance ===================== */ +/* + * Single-tone HB: solve the periodic steady state in the FREQUENCY domain. Each + * node voltage is a truncated Fourier series V(t) = sum_{k=-K..K} V_k e^{jk w0 t}; + * the KCL residual at every node/harmonic + * F_k = I_R,k(V) + [dq/dt]_k - Is_k = 0 + * is solved by Newton, with the (2K+1)N conversion matrix (E-121) as the Jacobian. + * - I_R(v(t_s)) : nonlinear RESISTIVE current, from a DC-mode device load at each of + * P time samples (residual current = G*v - rhs). + * - [dq/dt]_k : the REACTIVE current. dq/dt = C(v)*v' (chain rule), so its spectrum + * is the conversion matrix's reactive term (jm*w0*C_{k-m}) applied to V -- NONLINEAR + * charge is handled with NO per-device charge extraction, just the C(t) samples. + * - Is : the independent-source excitation spectrum (loaded at v=0, t=t_s). + * Reuses pac_build_matrix (Jacobian) and pss_csolve (dense complex Newton solve). + */ + +/* Sample the device residual + Jacobian at prescribed node voltages vsamp[s*N+(i-1)] + * (P samples). Fills hd with the G(t)/C(t) harmonics (h=0..2K) and returns the + * resistive-current harmonics in IRr/IRi (length (2K+1)*N). Returns 0 on success. */ +static int +hb_extract(CKTcircuit *ckt, const double *vsamp, int N, int P, int K, + struct pac_harm *hd, double *IRr, double *IRi) +{ + int H = 2 * K, i, r, c, e, h, nnz, s; + int *rr, *cc; + double *Gt, *Ct, *IRt, *cw, *sw, *Gmr, *Gmi, *Cmr, *Cmi, *bsave; + + memset(hd, 0, sizeof(*hd)); + if (P <= 0 || N <= 0 || K < 1) + return 1; + bsave = TMALLOC(double, N); + + spSetComplex(ckt->CKTmatrix->SPmatrix); + + /* establish structure at sample 0's bias */ + for (i = 1; i <= N; i++) + ckt->CKTrhsOld[i] = vsamp[i - 1]; + ckt->CKTrhsOld[0] = 0.0; + ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG; + CKTload(ckt); + ckt->CKTomega = 1.0; + ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEAC; + CKTacLoad(ckt); + + nnz = 0; + for (r = 1; r <= N; r++) + for (c = 1; c <= N; c++) + if (SMPfindElt(ckt->CKTmatrix, r, c, 0)) + nnz++; + if (nnz <= 0) + return 1; + rr = TMALLOC(int, nnz); + cc = TMALLOC(int, nnz); + e = 0; + for (r = 1; r <= N; r++) + for (c = 1; c <= N; c++) + if (SMPfindElt(ckt->CKTmatrix, r, c, 0)) { rr[e] = r; cc[e] = c; e++; } + + Gt = TMALLOC(double, (size_t)nnz * (size_t)P); + Ct = TMALLOC(double, (size_t)nnz * (size_t)P); + IRt = TMALLOC(double, (size_t)N * (size_t)P); + + for (s = 0; s < P; s++) { + double *Gv; + /* DC-mode load at v(t_s): matrix real = G, rhs = resistive companion */ + for (i = 1; i <= N; i++) + ckt->CKTrhsOld[i] = vsamp[(size_t)s * (size_t)N + (size_t)(i - 1)]; + ckt->CKTrhsOld[0] = 0.0; + for (i = 0; i <= N; i++) + ckt->CKTrhs[i] = 0.0; + ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEDCOP | MODEINITSMSIG; + CKTload(ckt); + /* companion source b = G*v - i(v) is in CKTrhs NOW; save it before acLoad + * clears it (so the resistive current is i(v) = G*v - b, the ACTUAL current, + * not the tangent G*v). */ + for (i = 1; i <= N; i++) + bsave[i - 1] = ckt->CKTrhs[i]; + ckt->CKTomega = 1.0; + ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODEAC; + CKTacLoad(ckt); /* clears + stamps G (real) and C (imag) cleanly */ + + /* resistive current I_R = G*v - b, using the clean G from acLoad */ + Gv = IRt + (size_t)s * (size_t)N; /* reuse row s as scratch, then subtract */ + for (i = 0; i < N; i++) + Gv[i] = 0.0; + for (e = 0; e < nnz; e++) { + double *el = (double *) SMPfindElt(ckt->CKTmatrix, rr[e], cc[e], 0); + double g = el ? el[0] : 0.0; + Gt[(size_t)e * (size_t)P + (size_t)s] = g; + Ct[(size_t)e * (size_t)P + (size_t)s] = el ? el[1] : 0.0; + Gv[rr[e] - 1] += g * vsamp[(size_t)s * (size_t)N + (size_t)(cc[e] - 1)]; + } + for (i = 1; i <= N; i++) + Gv[i - 1] -= bsave[i - 1]; /* I_R = G*v - b = i(v) */ + } + + /* DFT G(t), C(t) -> harmonics; I_R(t) -> IRr/IRi (harmonics -K..K packed 0..2K) */ + cw = TMALLOC(double, (size_t)(H + 1) * (size_t)P); + sw = TMALLOC(double, (size_t)(H + 1) * (size_t)P); + for (h = 0; h <= H; h++) + for (s = 0; s < P; s++) { + double ang = 2.0 * M_PI * (double)h * (double)s / (double)P; + cw[(size_t)h * (size_t)P + (size_t)s] = cos(ang); + sw[(size_t)h * (size_t)P + (size_t)s] = sin(ang); + } + Gmr = TMALLOC(double, (size_t)nnz * (size_t)(H + 1)); + Gmi = TMALLOC(double, (size_t)nnz * (size_t)(H + 1)); + Cmr = TMALLOC(double, (size_t)nnz * (size_t)(H + 1)); + Cmi = TMALLOC(double, (size_t)nnz * (size_t)(H + 1)); + for (e = 0; e < nnz; e++) + for (h = 0; h <= H; h++) { + double gr = 0, gi = 0, cr = 0, ci = 0; + for (s = 0; s < P; s++) { + double cs = cw[(size_t)h * (size_t)P + (size_t)s]; + double sn = sw[(size_t)h * (size_t)P + (size_t)s]; + double gv = Gt[(size_t)e * (size_t)P + (size_t)s]; + double cv = Ct[(size_t)e * (size_t)P + (size_t)s]; + gr += gv * cs; gi -= gv * sn; + cr += cv * cs; ci -= cv * sn; + } + Gmr[(size_t)e * (size_t)(H + 1) + (size_t)h] = gr / (double)P; + Gmi[(size_t)e * (size_t)(H + 1) + (size_t)h] = gi / (double)P; + Cmr[(size_t)e * (size_t)(H + 1) + (size_t)h] = cr / (double)P; + Cmi[(size_t)e * (size_t)(H + 1) + (size_t)h] = ci / (double)P; + } + /* I_R harmonics: full k=-K..K (row (k+K)*N + node) */ + for (r = 0; r < N; r++) + for (h = -K; h <= K; h++) { + double xr = 0, xi = 0; + int hh = h < 0 ? -h : h; + for (s = 0; s < P; s++) { + double cs = cw[(size_t)hh * (size_t)P + (size_t)s]; + double sn = sw[(size_t)hh * (size_t)P + (size_t)s]; + double x = IRt[(size_t)s * (size_t)N + (size_t)r]; + if (h >= 0) { xr += x * cs; xi -= x * sn; } + else { xr += x * cs; xi += x * sn; } /* conj for -h */ + } + IRr[(size_t)(h + K) * (size_t)N + (size_t)r] = xr / (double)P; + IRi[(size_t)(h + K) * (size_t)N + (size_t)r] = xi / (double)P; + } + + FREE(Gt); FREE(Ct); FREE(IRt); FREE(cw); FREE(sw); FREE(bsave); + hd->N = N; hd->M = K; hd->H = H; hd->nnz = nnz; hd->Ntot = (2*K + 1) * N; + hd->rr = rr; hd->cc = cc; + hd->Gmr = Gmr; hd->Gmi = Gmi; hd->Cmr = Cmr; hd->Cmi = Cmi; + hd->B0r = NULL; hd->B0i = NULL; hd->has_src = 0; + return 0; +} + +int +HBanalyze(CKTcircuit *ckt, double f0, int K, int Pin, int maxiter, double tol, int verbose) +{ + int N = SMPmatSize(ckt->CKTmatrix); + int P = Pin > 0 ? Pin : ((8 * K < 32) ? 32 : 8 * K); + int Ntot = (2 * K + 1) * N; + int iter, s, i, k, rc = 0; + double T = 1.0 / f0, w0 = 2.0 * M_PI * f0; + double *Vr, *Vi, *vsamp, *IRr, *IRi, *Isr, *Isi, *Fr, *Fi, *Jr, *Ji, *Kr, *Ki; + struct pac_harm hd; + + if (N <= 0 || K < 1) { fprintf(stderr, "HB: bad size.\n"); return E_PARMVAL; } + if ((2 * K + 1) * N > 900) { + fprintf(stderr, "HB: system %dx%d too large for the dense solver.\n", Ntot, Ntot); + return E_PARMVAL; + } + + Vr = TMALLOC(double, Ntot); Vi = TMALLOC(double, Ntot); + IRr = TMALLOC(double, Ntot); IRi = TMALLOC(double, Ntot); + Isr = TMALLOC(double, Ntot); Isi = TMALLOC(double, Ntot); + Fr = TMALLOC(double, Ntot); Fi = TMALLOC(double, Ntot); + Kr = TMALLOC(double, Ntot); Ki = TMALLOC(double, Ntot); + Jr = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); Ji = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); + vsamp = TMALLOC(double, (size_t)N * (size_t)P); + + /* --- source spectrum Is_k: load at v=0, sources evaluated at t_s --- */ + { + double *ist = TMALLOC(double, (size_t)N * (size_t)P); + for (s = 0; s < P; s++) { + for (i = 0; i <= N; i++) { ckt->CKTrhsOld[i] = 0.0; ckt->CKTrhs[i] = 0.0; } + ckt->CKTtime = (double)s * T / (double)P; + ckt->CKTmode = (ckt->CKTmode & MODEUIC) | MODETRAN | MODEINITTRAN; + CKTload(ckt); + for (i = 1; i <= N; i++) + ist[(size_t)s * (size_t)N + (size_t)(i - 1)] = ckt->CKTrhs[i]; + } + for (i = 0; i < N; i++) + for (k = -K; k <= K; k++) { + double xr = 0, xi = 0; int hh = k < 0 ? -k : k; + for (s = 0; s < P; s++) { + double ang = 2.0 * M_PI * hh * s / (double)P; + double x = ist[(size_t)s * (size_t)N + (size_t)i]; + xr += x * cos(ang); + xi += (k >= 0 ? -1.0 : 1.0) * x * sin(ang); + } + Isr[(size_t)(k + K) * (size_t)N + (size_t)i] = xr / P; + Isi[(size_t)(k + K) * (size_t)N + (size_t)i] = xi / P; + } + FREE(ist); + } + + /* --- Newton iterations --- */ + for (iter = 0; iter < maxiter; iter++) { + double fnorm = 0.0; + /* v(t_s) = Re sum_k V_k e^{j k w0 t_s} */ + for (s = 0; s < P; s++) + for (i = 0; i < N; i++) { + double v = 0.0; + for (k = -K; k <= K; k++) { + double ang = 2.0 * M_PI * k * s / (double)P; + v += Vr[(size_t)(k + K) * (size_t)N + (size_t)i] * cos(ang) + - Vi[(size_t)(k + K) * (size_t)N + (size_t)i] * sin(ang); + } + vsamp[(size_t)s * (size_t)N + (size_t)i] = v; + } + + if (hb_extract(ckt, vsamp, N, P, K, &hd, IRr, IRi)) { + fprintf(stderr, "HB: device extraction failed.\n"); rc = E_PARMVAL; break; + } + + /* build the full Jacobian J = G + jwC conversion matrix */ + pac_build_matrix(&hd, f0, 0.0, Jr, Ji); + + /* reactive current I_C = (J - Jg)*V where Jg is the resistive (G-only) + * conversion matrix -- i.e. the jwC part of J applied to V. */ + { + struct pac_harm hg = hd; + double *Jgr = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); + double *Jgi = TMALLOC(double, (size_t)Ntot * (size_t)Ntot); + hg.Cmr = TMALLOC(double, (size_t)hd.nnz * (size_t)(hd.H + 1)); /* zero C -> resistive only */ + hg.Cmi = TMALLOC(double, (size_t)hd.nnz * (size_t)(hd.H + 1)); + pac_build_matrix(&hg, f0, 0.0, Jgr, Jgi); + FREE(hg.Cmr); FREE(hg.Cmi); + /* I_C = (J - Jg) * V */ + for (i = 0; i < Ntot; i++) { + double cr = 0, ci = 0; + for (k = 0; k < Ntot; k++) { + double ar = Jr[(size_t)i * (size_t)Ntot + (size_t)k] - Jgr[(size_t)i * (size_t)Ntot + (size_t)k]; + double ai = Ji[(size_t)i * (size_t)Ntot + (size_t)k] - Jgi[(size_t)i * (size_t)Ntot + (size_t)k]; + cr += ar * Vr[k] - ai * Vi[k]; + ci += ar * Vi[k] + ai * Vr[k]; + } + Kr[i] = cr; Ki[i] = ci; /* Kr/Ki = reactive current I_C */ + } + FREE(Jgr); FREE(Jgi); + } + + /* residual F = I_R + I_C - Is */ + for (i = 0; i < Ntot; i++) { + Fr[i] = IRr[i] + Kr[i] - Isr[i]; + Fi[i] = IRi[i] + Ki[i] - Isi[i]; + fnorm += Fr[i] * Fr[i] + Fi[i] * Fi[i]; + } + fnorm = sqrt(fnorm); + if (verbose) + fprintf(stderr, "HB iter %2d: |F| = %.6e\n", iter, fnorm); + + /* Newton step: J * dV = -F */ + for (i = 0; i < Ntot; i++) { Fr[i] = -Fr[i]; Fi[i] = -Fi[i]; } + if (pss_csolve(Ntot, Jr, Ji, Fr, Fi)) { /* Fr/Fi <- dV */ + fprintf(stderr, "HB: singular Jacobian.\n"); rc = E_SINGULAR; + pac_free_harmonics(&hd); break; + } + for (i = 0; i < Ntot; i++) { Vr[i] += Fr[i]; Vi[i] += Fi[i]; } + pac_free_harmonics(&hd); + + if (fnorm < tol) { + fprintf(stdout, "HB: converged in %d iterations (|F| = %.3e).\n", iter + 1, fnorm); + break; + } + } + + /* --- output: labelled spectrum table, magnitude per node per harmonic --- */ + { + int numNames, error; + IFuid *nameList = NULL; + error = CKTnames(ckt, &numNames, &nameList); + fprintf(stdout, "\nHB: harmonic-balance spectrum (f0 = %g Hz, %d harmonics)\n" + " node harmonic frequency [Hz] |V| phase [deg]\n", + f0, K); + for (i = 0; i < N; i++) { + const char *nm = (!error && i < numNames) ? (const char *) nameList[i] : "?"; + for (k = 0; k <= K; k++) { + double sc = (k == 0) ? 1.0 : 2.0; /* single-sided amplitude */ + double vr = sc * Vr[(size_t)(k + K) * (size_t)N + (size_t)i]; + double vi = sc * Vi[(size_t)(k + K) * (size_t)N + (size_t)i]; + fprintf(stdout, " %-8s %6d %16.6e %14.6e %10.3f\n", + nm, k, k * f0, hypot(vr, vi), + (k == 0) ? 0.0 : atan2(vi, vr) * 180.0 / M_PI); + } + } + if (nameList) tfree(nameList); + } + (void) w0; + + FREE(Vr); FREE(Vi); FREE(IRr); FREE(IRi); FREE(Isr); FREE(Isi); + FREE(Fr); FREE(Fi); FREE(Jr); FREE(Ji); FREE(Kr); FREE(Ki); FREE(vsamp); + return rc; +} + + int DCpss(CKTcircuit *ckt, int restart) /* forced restart flag */ diff --git a/visualc/vngspice.vcxproj b/visualc/vngspice.vcxproj index 720fc0a73..b45c11426 100644 --- a/visualc/vngspice.vcxproj +++ b/visualc/vngspice.vcxproj @@ -886,6 +886,7 @@ + @@ -1503,6 +1504,7 @@ +