diff --git a/ChangeLog b/ChangeLog index 1fd8831d4..042c55878 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,24 @@ +2011-08-04 + + * src/frontend/commands.c , + * src/frontend/runcoms.c , + * src/frontend/runcoms.h , + * src/frontend/shyu.c , + * src/frontend/spiceif.c , + * src/include/Makefile.am , + * src/include/cktdefs.h , + * src/spicelib/analysis/Makefile.am , + * src/spicelib/analysis/analysis.c , + * src/spicelib/parser/inp2dot.c , + * (new) src/include/pssdefs.h , + + * (new) src/spicelib/analysis/dcpss.c , + * (new) src/spicelib/analysis/pssaskq.c , + * (new) src/spicelib/analysis/pssinit.c , + * (new) src/spicelib/analysis/psssetp.c : + + paolo, please fix this <------------- + 2011-08-04 Robert Larice * src/spicelib/analysis/cktsetbk.c : bug fix, broken insert at the front of CKTbreaks[] diff --git a/src/frontend/commands.c b/src/frontend/commands.c index a92f2dd65..089b7812d 100644 --- a/src/frontend/commands.c +++ b/src/frontend/commands.c @@ -277,6 +277,12 @@ struct comm spcp_coms[] = { { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, NULL, "[.tran line args] : Do a transient analysis." } , +/* SP: Steady State Analysis */ + { "pss", com_pss, TRUE, TRUE, + { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, + NULL, + "[.pss line args] : Do a periodic state analysis." } , +/* SP */ { "ac", com_ac, TRUE, TRUE, { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, NULL, @@ -683,6 +689,12 @@ struct comm nutcp_coms[] = { { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, NULL, "[.tran line args] : Do a transient analysis." } , +/* SP: Steady State Analysis */ + { "pss", NULL, TRUE, TRUE, + { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, + NULL, + "[.pss line args] : Do a periodic steady state analysis." } , +/* SP */ { "ac", NULL, TRUE, TRUE, { 0, 0, 0, 0 }, E_DEFHMASK, 0, LOTS, NULL, diff --git a/src/frontend/runcoms.c b/src/frontend/runcoms.c index 74225db1a..b5e94ee17 100644 --- a/src/frontend/runcoms.c +++ b/src/frontend/runcoms.c @@ -169,6 +169,15 @@ com_noise(wordlist *wl) return; } +/* SP: Steady State Analysis */ +void +com_pss(wordlist *wl) +{ + dosim("pss", wl); + return; +} +/* SP */ + static int dosim( char *what, /* in: command (pz,op,dc,ac,tf,tran,sens,disto,noise,run) */ diff --git a/src/frontend/runcoms.h b/src/frontend/runcoms.h index 5dab2771b..4d6f3715f 100644 --- a/src/frontend/runcoms.h +++ b/src/frontend/runcoms.h @@ -14,6 +14,9 @@ void com_dc(wordlist *wl); void com_ac(wordlist *wl); void com_tf(wordlist *wl); void com_tran(wordlist *wl); +/* SP: Stady State Analysis */ +void com_pss(wordlist *wl); +/* SP */ void com_sens(wordlist *wl); void com_disto(wordlist *wl); void com_noise(wordlist *wl); diff --git a/src/frontend/shyu.c b/src/frontend/shyu.c index ba7e3c0ab..99ad703f8 100644 --- a/src/frontend/shyu.c +++ b/src/frontend/shyu.c @@ -32,6 +32,7 @@ if_sens_run(CKTcircuit *ckt, wordlist *args, INPtables *tab) JOB *opJob; JOB *dcJob; JOB *tranJob; + JOB *pssJob; card *current; IFvalue ptemp; IFvalue *parm; @@ -285,6 +286,51 @@ uic: } } + /* *********************** */ + /* PSS - Spertica - 100910 */ + /* *********************** */ + if(strcmp(token ,"pss")==0){ + which = -1; + for(j=0;jnumAnalyses;j++) { + if(strcmp(ft_sim->analyses[j]->name,"PSS")==0) { + which = j; + break; + } + } + if(which != -1) { + err = (*(ft_sim->newAnalysis))(ft_curckt->ci_ckt,which,"pssan", + &(tranJob),ft_curckt->ci_specTask); + if(err) { + ft_sperror(err,"createPSS"); + return(0); + } + } + else{ + current->error = INPerrCat(current->error,INPmkTemp( + "periodic steady state analysis unsupported\n")); + } + + parm=INPgetValue(ckt,&line,IF_REAL,tab);/* Guessed Frequency */ + error = INPapName(ckt,which,pssJob,"fguess",parm); + if(error) current->error = INPerrCat(current->error, + INPerror(error)); + + parm = INPgetValue(ckt,&line,IF_REAL,tab); /* Stabilization time */ + error = INPapName(ckt,which,pssJob,"stabtime",parm); + if(error) current->error = INPerrCat(current->error, + INPerror(error)); + + parm = INPgetValue(ckt,&line,IF_INTEGER,tab); /* PSS points */ + error = INPapName(ckt,which,pssJob,"points",parm); + if(error) current->error = INPerrCat(current->error, + INPerror(error)); + + parm = INPgetValue(ckt,&line,IF_INTEGER,tab); /* PSS points */ + error = INPapName(ckt,which,pssJob,"harmonics",parm); + if(error) current->error = INPerrCat(current->error, + INPerror(error)); + } + next: while(*line) { /* read the entire line */ if(flag){ diff --git a/src/frontend/spiceif.c b/src/frontend/spiceif.c index 01a9183b2..717648835 100644 --- a/src/frontend/spiceif.c +++ b/src/frontend/spiceif.c @@ -215,7 +215,10 @@ if_run(CKTcircuit *ckt, char *what, wordlist *args, INPtables *tab) || eq(what, "adjsen") || eq(what, "sens") || eq(what,"tf") - || eq(what, "noise")) + || eq(what, "noise") +/* SP: Steady State Analysis */ + || eq(what, "pss")) +/* SP */ { s = wl_flatten(args); /* va: tfree char's tmalloc'ed in wl_flatten */ (void) sprintf(buf, ".%s", s); @@ -329,6 +332,9 @@ ci_specTask will point to it CDHW*/ ||(eq(what, "adjsen")) ||(eq(what, "sens")) ||(eq(what,"tf")) +/* SP: Steady State Analysis */ + ||(eq(what, "pss")) +/* SP */ ||(eq(what, "run")) ) { /*CDHW Run the analysis pointed to by ci_curTask CDHW*/ diff --git a/src/include/Makefile.am b/src/include/Makefile.am index 838e26c99..4c0f0d404 100644 --- a/src/include/Makefile.am +++ b/src/include/Makefile.am @@ -93,6 +93,7 @@ noinst_HEADERS = \ plot.h \ pnode.h \ profile.h \ + pssdefs.h \ pzdefs.h \ sen2defs.h \ sensdefs.h \ diff --git a/src/include/cktdefs.h b/src/include/cktdefs.h index 4b14864e2..19bf463fb 100644 --- a/src/include/cktdefs.h +++ b/src/include/cktdefs.h @@ -265,6 +265,15 @@ struct CKTcircuit { Enh_Ckt_Data_t *enh; /* data used by general enhancements */ #endif /* gtri - evt - wbk - 5/20/91 - add event-driven and enhancements data */ +/* SP: Steady State Analysis - 100609 */ + double CKTstabTime; /* PSS stab time */ + double CKTguessedFreq; /* PSS guessed frequency */ + int CKTharms; /* PSS harmonics */ + long int CKTpsspoints; /* PSS number of samples */ + char *CKToscNode; /* PSS oscnode */ + double CKTsteady_coeff; + int CKTsc_iter; +/* SP: 100609 */ }; @@ -379,6 +388,12 @@ extern int TFsetParm(CKTcircuit *, JOB *, int , IFvalue *); extern int TRANaskQuest(CKTcircuit *, JOB *, int , IFvalue *); extern int TRANsetParm(CKTcircuit *, JOB *, int , IFvalue *); extern int TRANinit(CKTcircuit *, JOB *); +/* SP: Steady State Analysis */ +extern int PSSaskQuest( CKTcircuit *, JOB *, int , IFvalue *); +extern int PSSsetParm( CKTcircuit *, JOB *, int , IFvalue *); +extern int PSSinit(CKTcircuit *, JOB *); +extern int DCpss( CKTcircuit *, int ); +/* SP */ extern int NaskQuest(CKTcircuit *, JOB *, int, IFvalue *); extern int NsetParm(CKTcircuit *, JOB *, int, IFvalue *); extern int NIacIter(CKTcircuit *); diff --git a/src/include/pssdefs.h b/src/include/pssdefs.h new file mode 100755 index 000000000..3e1f6cc1c --- /dev/null +++ b/src/include/pssdefs.h @@ -0,0 +1,38 @@ +/********** +Author: 2010-05 Stefano Perticaroli ``spertica'' +**********/ + +#ifndef PSS +#define PSS + +#include "jobdefs.h" +#include "tskdefs.h" + /* + * PSSdefs.h - defs for pss analyses + */ + +typedef struct { + int JOBtype; + JOB *JOBnextJob; + char *JOBname; + double PSSguessedFreq; + CKTnode *PSSoscNode; + double PSSstabTime; + long PSSmode; + long int PSSpoints; + int PSSharms; + void *PSSplot_td; + void *PSSplot_fd; + int sc_iter; + double steady_coeff; +} PSSan; + +#define GUESSED_FREQ 1 +#define STAB_TIME 2 +#define OSC_NODE 3 +#define PSS_POINTS 4 +#define PSS_HARMS 5 +#define PSS_UIC 6 +#define SC_ITER 7 +#define STEADY_COEFF 8 +#endif /*PSS*/ diff --git a/src/spicelib/analysis/Makefile.am b/src/spicelib/analysis/Makefile.am index 4966f667f..43ad2e7ef 100644 --- a/src/spicelib/analysis/Makefile.am +++ b/src/spicelib/analysis/Makefile.am @@ -72,6 +72,7 @@ libckt_la_SOURCES = \ dcoaskq.c \ dcop.c \ dcosetp.c \ + dcpss.c \ dctaskq.c \ dctran.c \ dctrcurv.c \ @@ -85,6 +86,9 @@ libckt_la_SOURCES = \ ninteg.c \ noisean.c \ nsetparm.c \ + pssaskq.c \ + pssinit.c \ + psssetp.c \ pzan.c \ pzaskq.c \ pzsetp.c \ diff --git a/src/spicelib/analysis/analysis.c b/src/spicelib/analysis/analysis.c index c73786552..5fc78fd62 100644 --- a/src/spicelib/analysis/analysis.c +++ b/src/spicelib/analysis/analysis.c @@ -17,7 +17,9 @@ extern SPICEanalysis TFinfo; extern SPICEanalysis DISTOinfo; extern SPICEanalysis NOISEinfo; extern SPICEanalysis SENSinfo; - +/* SP: Periodic Steady State*/ +extern SPICEanalysis PSSinfo; +/* SP */ SPICEanalysis *analInfo[] = { &OPTinfo, @@ -30,6 +32,7 @@ SPICEanalysis *analInfo[] = { &DISTOinfo, &NOISEinfo, &SENSinfo, + &PSSinfo, }; diff --git a/src/spicelib/analysis/dcpss.c b/src/spicelib/analysis/dcpss.c new file mode 100755 index 000000000..c01ec9934 --- /dev/null +++ b/src/spicelib/analysis/dcpss.c @@ -0,0 +1,1455 @@ +/********** + Author: 2010-05 Stefano Perticaroli ``spertica'' + **********/ + +/* subroutine to do DC PSS analysis */ + +#include "ngspice.h" +#include "cktdefs.h" +#include "pssdefs.h" +#include "sperror.h" + +/* for FFT */ +#include "ftedefs.h" +#include "dvec.h" +#include "sim.h" +#include "fteparse.h" +#include "const.h" +#include "../../frontend/fourier.h" +#include "../../frontend/variable.h" + +#include "config.h" +#include "cktaccept.h" +#include "trandefs.h" +#include "fteext.h" +#include "missing_math.h" + +#ifdef XSPICE +/* gtri - add - wbk - Add headers */ +#include "miftypes.h" + +#include "evt.h" +#include "mif.h" +#include "evtproto.h" +#include "ipctiein.h" +/* gtri - end - wbk - Add headers */ +#endif + +#ifdef CLUSTER +#include "cluster.h" +#endif + +#ifdef HAS_WINDOWS /* hvogt 10.03.99, nach W. Mues */ +void SetAnalyse( char * Analyse, int Percent); +#endif + + +int +CKTfour(int, int, double *, double *, double *, double, double *, double *, double *, double *,double *); + +int +DCpss(CKTcircuit *ckt, int restart) +{ + int oscnNode; + int i; + double olddelta; + double delta; + double new; + double *temp; + double startdTime; + double startsTime; + double startlTime; + double startcTime; + double startkTime; + double startTime; + int startIters; + int converged; + int firsttime; + int error; + #ifdef WANT_SENSE2 + #ifdef SENSDEBUG + FILE *outsen; + #endif /* SENSDEBUG */ + #endif + int save_order; + long save_mode; + IFuid timeUid; + IFuid freqUid; + IFuid *nameList; + int numNames; + double maxstepsize=0.0; + + int ltra_num; + CKTnode *node; + #ifdef PARALLEL_ARCH + long type = MT_TRANAN, length = 1; + #endif /* PARALLEL_ARCH */ + #ifdef XSPICE + /* gtri - add - wbk - 12/19/90 - Add IPC stuff */ + Ipc_Boolean_t ipc_firsttime = IPC_TRUE; + Ipc_Boolean_t ipc_secondtime = IPC_FALSE; + Ipc_Boolean_t ipc_delta_cut = IPC_FALSE; + double ipc_last_time = 0.0; + double ipc_last_delta = 0.0; + /* gtri - end - wbk - 12/19/90 - Add IPC stuff */ + #endif + #ifdef CLUSTER + int redostep; + #endif + /* new variables - to be reorganized */ + double time_temp, gf_history[1024], rr_history[1024]; + int msize, cont_1=1, in_stabilization=1, in_pss=0,shooting_cycle_counter=0, k; + long int nextTime_count=0, ntc_start_sh=0; + double *RHS_copy_se, err, err_0=1.0e30, time_temp_0, delta_t; + double time_err_min_1=0, time_err_min_0=0, err_min_0=1.0e30, err_min_1, delta_0, delta_1; + double flag_tu_2, times_fft[8192], err_1=0, err_max, time_err_max; + int flag_tu_1=0, pss_cycle_counter=1, pss_points_cycle=0, i4, i5, k1,rest; + int count_1, count_2, count_3, count_4, count_5, count_6, count_7, dynamic_test=0; + double ntc_mv, ntc_vec[4], ntc_old, gf_last_0=1e+30, gf_last_1=313; + double err_last=0, f_proj, thd; + double *psstimes, *pssvalues, *pssValues, tv_01, tv_03, tv_04, + *pssfreqs, *pssmags, *pssphases, *pssnmags, *pssnphases, *pssResults, + *RHS_max, *RHS_min, err_conv_ref, *S_old, *S_diff; + + printf("Periodic Steady State analysis started.\n"); + + PSSan *job = (PSSan*) (ckt->CKTcurJob); + oscnNode = job->PSSoscNode->number; + printf("PSS guessed frequency %g.\n", ckt->CKTguessedFreq); + printf("PSS points %ld.\n", ckt->CKTpsspoints); + printf("PSS harmonics number %d.\n", ckt->CKTharms); + printf("PSS steady coefficient %g.\n", ckt->CKTsteady_coeff); + + /* set first delta time step and circuit time */ + delta=ckt->CKTstep; + ckt->CKTtime=ckt->CKTinitTime; + ckt->CKTfinalTime=ckt->CKTstabTime; + /* printf("initial delta: %g\n", ckt->CKTdelta); */ + + /* 100906 - Paste from dctran.c */ + if(restart || ckt->CKTtime == 0) { + delta=MIN(1/(ckt->CKTguessedFreq)/100,ckt->CKTstep); + + #ifdef STEPDEBUG + printf("delta = %g finalTime/200: %g CKTstep: %g\n",delta,ckt->CKTfinalTime/200,ckt->CKTstep); + #endif + /* begin LTRA code addition */ + if (ckt->CKTtimePoints != NULL) + FREE(ckt->CKTtimePoints); + + if (ckt->CKTstep >= ckt->CKTmaxStep) + maxstepsize = ckt->CKTstep; + else + maxstepsize = ckt->CKTmaxStep; + + ckt->CKTsizeIncr = 10; + ckt->CKTtimeIndex = -1; /* before the DC soln has been stored */ + ckt->CKTtimeListSize = (int)(1 / (ckt->CKTguessedFreq) / maxstepsize + 0.5); + ltra_num = CKTtypelook("LTRA"); + if (ltra_num >= 0 && ckt->CKThead[ltra_num] != NULL) + ckt->CKTtimePoints = NEWN(double, ckt->CKTtimeListSize); + /* end LTRA code addition */ + + if(ckt->CKTbreaks) FREE(ckt->CKTbreaks); + ckt->CKTbreaks = TMALLOC(double, 2); + if(ckt->CKTbreaks == NULL) return(E_NOMEM); + *(ckt->CKTbreaks)=0; + *(ckt->CKTbreaks+1)=ckt->CKTfinalTime; + ckt->CKTbreakSize=2; + + #ifdef XSPICE + /* gtri - begin - wbk - 12/19/90 - Modify setting of CKTminBreak */ + /* if(ckt->CKTminBreak==0) ckt->CKTminBreak=ckt->CKTmaxStep*5e-5; */ + /* Set to 10 times delmin for ATESSE 1 compatibity */ + if(ckt->CKTminBreak==0) ckt->CKTminBreak = 10.0 * ckt->CKTdelmin; + /* gtri - end - wbk - 12/19/90 - Modify setting of CKTminBreak */ + #else + if(ckt->CKTminBreak==0) ckt->CKTminBreak=ckt->CKTmaxStep*5e-5; + #endif + + #ifdef XSPICE + /* gtri - add - wbk - 12/19/90 - Add IPC stuff and set anal_init and anal_type */ + + /* Tell the beginPlot routine what mode we're in */ + g_ipc.anal_type = IPC_ANAL_TRAN; + + /* Tell the code models what mode we're in */ + g_mif_info.circuit.anal_type = MIF_DC; + + g_mif_info.circuit.anal_init = MIF_TRUE; + + /* gtri - end - wbk */ + #endif + error = CKTnames(ckt,&numNames,&nameList); + if(error) return(error); + SPfrontEnd->IFnewUid (ckt, &timeUid, NULL, + "time", UID_OTHER, NULL); + error = SPfrontEnd->OUTpBeginPlot (ckt, + ckt->CKTcurJob, + "Time Domain Periodic Steady State",timeUid,IF_REAL,numNames,nameList, + IF_REAL,&(((PSSan*)ckt->CKTcurJob)->PSSplot_td)); + tfree(nameList); + if(error) return(error); + + ckt->CKTtime = 0; + ckt->CKTdelta = 0; + ckt->CKTbreak=1; + firsttime = 1; + save_mode = (ckt->CKTmode&MODEUIC)|MODETRANOP | MODEINITJCT; + save_order = ckt->CKTorder; + #ifdef XSPICE + /* gtri - begin - wbk - set a breakpoint at end of supply ramping time */ + /* must do this after CKTtime set to 0 above */ + if(ckt->enh->ramp.ramptime > 0.0) + CKTsetBreak(ckt, ckt->enh->ramp.ramptime); + /* gtri - end - wbk - set a breakpoint at end of supply ramping time */ + + /* gtri - begin - wbk - Call EVTop if event-driven instances exist */ + if(ckt->evt->counts.num_insts != 0){ + /* use new DCOP algorithm */ + converged = EVTop(ckt, + (ckt->CKTmode & MODEUIC) | MODETRANOP | MODEINITJCT, + (ckt->CKTmode & MODEUIC)|MODETRANOP| MODEINITFLOAT, + ckt->CKTdcMaxIter, + MIF_TRUE); + EVTdump(ckt, IPC_ANAL_DCOP, 0.0); + + EVTop_save(ckt, MIF_FALSE, 0.0); + + /* gtri - end - wbk - Call EVTop if event-driven instances exist */ + } else + #endif + converged = CKTop(ckt, + (ckt->CKTmode & MODEUIC)|MODETRANOP| MODEINITJCT, + (ckt->CKTmode & MODEUIC)|MODETRANOP| MODEINITFLOAT, + ckt->CKTdcMaxIter); + #ifdef STEPDEBUG + if(converged != 0) { + fprintf(stdout,"\nTransient solution failed -\n"); + CKTncDump(ckt); + fprintf(stdout,"\n"); + fflush(stdout); + } else if (!ft_noacctprint && !ft_noinitprint) { + fprintf(stdout,"\nInitial Transient Solution\n"); + fprintf(stdout,"--------------------------\n\n"); + fprintf(stdout,"%-30s %15s\n", "Node", "Voltage"); + fprintf(stdout,"%-30s %15s\n", "----", "-------"); + for(node=ckt->CKTnodes->next;node;node=node->next) { + if (strstr(node->name, "#branch") || !strstr(node->name, "#")) + fprintf(stdout,"%-30s %15g\n", node->name, + *(ckt->CKTrhsOld+node->number)); + } + fprintf(stdout,"\n"); + fflush(stdout); + } + #endif + if(converged != 0) return(converged); + #ifdef XSPICE + /* gtri - add - wbk - 12/19/90 - Add IPC stuff */ + + /* Send the operating point results for Mspice compatibility */ + if(g_ipc.enabled) { + ipc_send_dcop_prefix(); + CKTdump(ckt,(double)0,(((PSSan*)ckt->CKTcurJob)->PSSplot_td)); + ipc_send_dcop_suffix(); + } + + /* gtri - end - wbk */ + + + /* gtri - add - wbk - 12/19/90 - set anal_init and anal_type */ + + g_mif_info.circuit.anal_init = MIF_TRUE; + + /* Tell the code models what mode we're in */ + g_mif_info.circuit.anal_type = MIF_TRAN; + + /* gtri - end - wbk */ + + /* gtri - begin - wbk - Add Breakpoint stuff */ + + /* Initialize the temporary breakpoint variables to infinity */ + g_mif_info.breakpoint.current = 1.0e30; + g_mif_info.breakpoint.last = 1.0e30; + + /* gtri - end - wbk - Add Breakpoint stuff */ + #endif + ckt->CKTstat->STATtimePts ++; + ckt->CKTorder=1; + for(i=0;i<7;i++) { + ckt->CKTdeltaOld[i]=ckt->CKTmaxStep; + } + ckt->CKTdelta = delta; + #ifdef STEPDEBUG + (void)printf("delta initialized to %g\n",ckt->CKTdelta); + #endif + ckt->CKTsaveDelta = ckt->CKTfinalTime/50; + + #ifdef WANT_SENSE2 + if(ckt->CKTsenInfo && (ckt->CKTsenInfo->SENmode & TRANSEN)){ + #ifdef SENSDEBUG + printf("\nTransient Sensitivity Results\n\n"); + CKTsenPrint(ckt); + #endif /* SENSDEBUG */ + save = ckt->CKTsenInfo->SENmode; + ckt->CKTsenInfo->SENmode = TRANSEN; + save1 = ckt->CKTmode; + save2 = ckt->CKTorder; + ckt->CKTmode = save_mode; + ckt->CKTorder = save_order; + if(error = CKTsenDCtran(ckt)) return(error); + ckt->CKTmode = save1; + ckt->CKTorder = save2; + } + #endif + /*print*/ + ckt->CKTmode = (ckt->CKTmode&MODEUIC)|MODETRAN | MODEINITTRAN; + /* modeinittran set here */ + ckt->CKTag[0]=ckt->CKTag[1]=0; + bcopy(ckt->CKTstate0, ckt->CKTstate1, + (size_t) ckt->CKTnumStates * sizeof(double)); + + #ifdef WANT_SENSE2 + if(ckt->CKTsenInfo && (ckt->CKTsenInfo->SENmode & TRANSEN)){ + size = SMPmatSize(ckt->CKTmatrix); + for(i = 1; i<=size ; i++) + *(ckt->CKTrhsOp + i) = *(ckt->CKTrhsOld + i); + } + #endif + + startTime = SPfrontEnd->IFseconds(); + startIters = ckt->CKTstat->STATnumIter; + startdTime = ckt->CKTstat->STATdecompTime; + startsTime = ckt->CKTstat->STATsolveTime; + startlTime = ckt->CKTstat->STATloadTime; + startcTime = ckt->CKTstat->STATcombineTime; + startkTime = ckt->CKTstat->STATsyncTime; + #ifdef CLUSTER + CLUsetup(ckt); + #endif + } else { + /*saj As traninit resets CKTmode */ + ckt->CKTmode = (ckt->CKTmode&MODEUIC)|MODETRAN | MODEINITPRED; + /* saj */ + startTime = SPfrontEnd->IFseconds(); + startIters = ckt->CKTstat->STATnumIter; + startdTime = ckt->CKTstat->STATdecompTime; + startsTime = ckt->CKTstat->STATsolveTime; + startlTime = ckt->CKTstat->STATloadTime; + startcTime = ckt->CKTstat->STATcombineTime; + startkTime = ckt->CKTstat->STATsyncTime; + if(ckt->CKTminBreak==0) ckt->CKTminBreak=ckt->CKTmaxStep*5e-5; + firsttime=0; + /* To get rawfile working saj*/ + error = SPfrontEnd->OUTpBeginPlot + (NULL, NULL, NULL, NULL, 0, 666, NULL, 666, &(((PSSan*)ckt->CKTcurJob)->PSSplot_td)); + if(error) { + fprintf(stderr, "Couldn't relink rawfile\n"); + return error; + } + /*end saj*/ + goto resume; + } + + /* 650 */ + nextTime: + + nextTime_count=nextTime_count+1; + /* Does not start from initial time avoiding IC issues */ + if ( nextTime_count>=10 && nextTime_count<8202 && in_stabilization) { + times_fft[nextTime_count-10]=ckt->CKTtime; + } + + /* begin LTRA code addition */ + if (ckt->CKTtimePoints) { + ckt->CKTtimeIndex++; + if (ckt->CKTtimeIndex >= ckt->CKTtimeListSize) { + /* need more space */ + int need; + if (in_stabilization) { + need = (int)(0.5 + (ckt->CKTstabTime - ckt->CKTtime) / maxstepsize); /* FIXME, ceil ? */ + } else { + need = (int)(0.5 + (time_temp + (1 / ckt->CKTguessedFreq) - ckt->CKTtime) / maxstepsize); + } + if (need < ckt->CKTsizeIncr) + need = ckt->CKTsizeIncr; + ckt->CKTtimeListSize += need; + ckt->CKTtimePoints = TREALLOC(double, ckt->CKTtimePoints, ckt->CKTtimeListSize); + ckt->CKTsizeIncr = (int) ceil(1.4 * ckt->CKTsizeIncr); + } + *(ckt->CKTtimePoints + ckt->CKTtimeIndex) = ckt->CKTtime; + } + /* end LTRA code addition */ + error = CKTaccept(ckt); + /* check if current breakpoint is outdated; if so, clear */ + if (ckt->CKTtime > *(ckt->CKTbreaks)) CKTclrBreak(ckt); + + /* + * Breakpoint handling scheme: + * When a timepoint t is accepted (by CKTaccept), clear all previous + * breakpoints, because they will never be needed again. + * + * t may itself be a breakpoint, or indistinguishably close. DON'T + * clear t itself; recognise it as a breakpoint and act accordingly + * + * if t is not a breakpoint, limit the timestep so that the next + * breakpoint is not crossed + */ + + #ifdef STEPDEBUG + printf("Delta %g accepted at time %g (finaltime: %g)\n",ckt->CKTdelta,ckt->CKTtime,ckt->CKTfinalTime); + fflush(stdout); + #endif /* STEPDEBUG */ + ckt->CKTstat->STATaccepted ++; + ckt->CKTbreak=0; + /* XXX Error will cause single process to bail. */ + if(error) { + ckt->CKTcurrentAnalysis = DOING_TRAN; + ckt->CKTstat->STATtranTime += SPfrontEnd->IFseconds() - startTime; + ckt->CKTstat->STATtranIter += ckt->CKTstat->STATnumIter - startIters; + ckt->CKTstat->STATtranDecompTime += ckt->CKTstat->STATdecompTime - + startdTime; + ckt->CKTstat->STATtranSolveTime += ckt->CKTstat->STATsolveTime - + startsTime; + ckt->CKTstat->STATtranLoadTime += ckt->CKTstat->STATloadTime - + startlTime; + ckt->CKTstat->STATtranCombTime += ckt->CKTstat->STATcombineTime - + startcTime; + ckt->CKTstat->STATtranSyncTime += ckt->CKTstat->STATsyncTime - + startkTime; + return(error); + } + #ifdef XSPICE + /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */ + + if(g_ipc.enabled) { + + if ( in_pss && pss_cycle_counter==1 ) { + /* Send event-driven results */ + EVTdump(ckt, IPC_ANAL_TRAN, 0.0); + + /* Then follow with analog results... */ + + /* Test to see if delta was cut by a breakpoint, */ + /* a non-convergence, or a too large truncation error */ + if(ipc_firsttime) + ipc_delta_cut = IPC_FALSE; + else if(ckt->CKTtime < (ipc_last_time + (0.999 * ipc_last_delta))) + ipc_delta_cut = IPC_TRUE; + else + ipc_delta_cut = IPC_FALSE; + + /* Record the data required to check for delta cuts */ + ipc_last_time = ckt->CKTtime; + ipc_last_delta = MIN(ckt->CKTdelta, ckt->CKTmaxStep); + + /* Send results data if time since last dump is greater */ + /* than 'mintime', or if first or second timepoints, */ + /* or if delta was cut */ + if( (ckt->CKTtime >= (g_ipc.mintime + g_ipc.last_time)) || + ipc_firsttime || ipc_secondtime || ipc_delta_cut ) { + ipc_send_data_prefix(ckt->CKTtime); + CKTdump(ckt,ckt->CKTtime, + (((PSSan*)ckt->CKTcurJob)->PSSplot_td)); + ipc_send_data_suffix(); + + if(ipc_firsttime) { + ipc_firsttime = IPC_FALSE; + ipc_secondtime = IPC_TRUE; + } + else if(ipc_secondtime) + ipc_secondtime = IPC_FALSE; + + g_ipc.last_time = ckt->CKTtime; + } + } + } + else + /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */ + #endif + #ifdef CLUSTER + if ( in_pss && pss_cycle_counter==1 ) { + CLUoutput(ckt); + } + #endif + if ( in_pss && pss_cycle_counter==1 ) { + if(ckt->CKTtime >= ckt->CKTinitTime) CKTdump( ckt, ckt->CKTtime, ( ((PSSan*)ckt->CKTcurJob)->PSSplot_td) ); + *(psstimes+pss_points_cycle)=ckt->CKTtime; + for(count_1=1; count_1CKTrhsOld+count_1); + } + pss_points_cycle++; + } + #ifdef XSPICE + /* gtri - begin - wbk - Update event queues/data for accepted timepoint */ + /* Note: this must be done AFTER sending results to SI so it can't */ + /* go next to CKTaccept() above */ + if(ckt->evt->counts.num_insts > 0) + EVTaccept(ckt, ckt->CKTtime); + /* gtri - end - wbk - Update event queues/data for accepted timepoint */ + #endif + ckt->CKTstat->STAToldIter = ckt->CKTstat->STATnumIter; + + /* ***********************************/ + /* ******* SHOOTING CODE BLOCK *******/ + /* ***********************************/ + if (in_stabilization) { + /* test if stabTime has been reached */ + if ( AlmostEqualUlps( ckt->CKTtime , ckt->CKTstabTime, 100 ) ) + { + /* use the 'dirty hack' to get near the fundamental frequency */ + if (nextTime_count<8192) { + for(count_2=1; count_20) { + delta_t += (times_fft[count_2]-times_fft[count_2-1])/nextTime_count; + } + } + } else { + for(count_2=1; count_2<8192; count_2++) { + if (times_fft[count_2]>0) { + delta_t += (times_fft[count_2]-times_fft[count_2-1])/8192; + } + } + } + if (ckt->CKTguessedFreq>1.1/delta_t/10 || ckt->CKTguessedFreq<0.9/delta_t/10) { + ckt->CKTguessedFreq=1.0/delta_t/10; /*** FREQUENCY INITIAL GUESS ***/ + printf("Frequency initial guess changed to %g from stabilization transient analysis.\n",ckt->CKTguessedFreq); + } + cont_1=1; + time_temp=ckt->CKTtime; + ckt->CKTfinalTime=time_temp+2/(ckt->CKTguessedFreq); + /* set the first requested breakpoint */ + CKTsetBreak(ckt, time_temp+1/(ckt->CKTguessedFreq)); + printf("Exiting from stabilization.\n"); + printf("Time of first shooting evaluation will be %1.10g.\n",time_temp+1/(ckt->CKTguessedFreq)); + /* next time is no more in stab - unset the flag */ + in_stabilization=0; + /* get matrix size and allocate memory */ + msize = SMPmatSize(ckt->CKTmatrix); + RHS_copy_se = TMALLOC(double, msize); + RHS_max = TMALLOC(double, msize); + RHS_min = TMALLOC(double, msize); + S_old = TMALLOC(double, msize); + S_diff = TMALLOC(double, msize); + /* print RHS on exiting from stab */ + printf("RHS on exiting from stabilization: "); + for(count_3 = 1; count_3 <= msize; count_3++) { + *(RHS_copy_se+count_3-1)= *(ckt->CKTrhsOld+count_3); + printf("%-15g ", *(RHS_copy_se+count_3-1)); + } + printf("\n"); + } + + /* ELSE not in stabilization but in shooting */ + } else if ( !in_pss ) { + + /* error norms of RHS solution on every accepted nextTime if out of stabilization */ + err=0; + for(count_4 = 1; count_4 <= msize; count_4++) + { + double diff = *(ckt->CKTrhsOld+count_4) - *(RHS_copy_se+count_4-1); + err += diff * diff; + /* save max and min per node or branch on every estimated period */ + if (*(RHS_max+count_4-1) < *(ckt->CKTrhsOld+count_4)) + *(RHS_max+count_4-1) = *(ckt->CKTrhsOld+count_4); + if (*(RHS_min+count_4-1) > *(ckt->CKTrhsOld+count_4)) + *(RHS_min+count_4-1) = *(ckt->CKTrhsOld+count_4); + } + err=sqrt(err); + /*** frequency projection ***/ + f_proj=(err-err_last)*(ckt->CKTguessedFreq); + + err_last=err; + + /* Start frequency estimation */ + if (errCKTtime>=time_temp+0.5/ckt->CKTguessedFreq) /* far enough from time temp... */ + { + if (errCKTtime; /* store minimum of RHS vector error time */ + delta_1=delta_0; + delta_0=ckt->CKTdelta; + } + } + err_0=err; + + if (err>err_1 && ckt->CKTtime>=time_temp+0.1/ckt->CKTguessedFreq) /* far enough from time temp... */ + { + if (err>err_max) { + err_max=err; /* store maximum of RHS vector error */ + time_err_max=ckt->CKTtime; /* store maximum of RHS vector error time */ + } + } + err_1=err; + + /* if evolution is near shooting */ + if ( AlmostEqualUlps( ckt->CKTtime, time_temp+1/ckt->CKTguessedFreq, 10 ) || (ckt->CKTtime > time_temp+1/ckt->CKTguessedFreq) ) + { + if (shooting_cycle_counter == 0) { + /* If first time in shooting we warn about that ! */ + ntc_start_sh=nextTime_count; + printf("In shooting...\n"); + } + + /* Take mean value of number of next time steps in one shooting evaluation - 4 frame window */ + ntc_vec[3]=ntc_vec[2]; + ntc_vec[2]=ntc_vec[1]; + ntc_vec[1]=ntc_vec[0]; + ntc_vec[0]=nextTime_count-ntc_old; + ntc_mv=(ntc_vec[0]+ntc_vec[1]+ntc_vec[2]+ntc_vec[3])*0.25; + ntc_old=nextTime_count; + printf("\n----------------\n"); + printf("Shooting cycle iteration number: %3d ||", shooting_cycle_counter); + printf("NTC_MV: %g || f_proj: %g\n", ntc_mv, f_proj); /* for debugging purpose */ + printf("Print of dynamically consistent nodes voltages or branches currents:\n"); + + for(i=1, node=ckt->CKTnodes->next; node; i++, node=node->next) { + + if (!strstr(node->name, "#")) { + + + + if (fabs(*(RHS_max+i-1))>fabs(*(RHS_min+i-1))) { + tv_01=fabs(*(RHS_max+i-1)); + } else { + tv_01=fabs(*(RHS_min+i-1)); + } + err_conv_ref += ( ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )*1e-3+1e-6)*7*ckt->CKTsteady_coeff; + if ( fabs( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) ) > 10*1e-6) + { + *(S_diff+i-1)= ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )/tv_01 - *(S_old+i-1); + *(S_old+i-1)= ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )/tv_01; + if(fabs(*(S_old+i-1))>0.1) printf("Node voltage %15s: RHS diff. %-15g || Conv. ref. %-15g || RHS max. %-15g || RHS min. %-15g || Snode %-15g\n",node->name, + (*(ckt->CKTrhsOld+i)-*(RHS_copy_se+i-1)), + ( ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )*1e-3+1e-6)*7*ckt->CKTsteady_coeff, + *(RHS_max+i-1), + *(RHS_min+i-1), + *(S_diff+i-1) + //( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) ) / ( (*(RHS_max+i-1)) + (*(RHS_min+i-1)) ) / 0.5 + ); + dynamic_test++; /* test on voltage dynamic consistence */ + } else + { + *(S_old+i-1)=0; + *(S_diff+i-1)=0; + } + } else { + + + + if (fabs(*(RHS_max+i-1))>fabs(*(RHS_min+i-1))) { + tv_01=fabs(*(RHS_max+i-1)); + } else { + tv_01=fabs(*(RHS_min+i-1)); + } + err_conv_ref += ( ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )*1e-3+1e-9)*7*ckt->CKTsteady_coeff; + if ( fabs( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) ) > 10*1e-9) + { + *(S_diff+i-1)= ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )/tv_01 - *(S_old+i-1); + *(S_old+i-1)= ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )/tv_01; + if(fabs(*(S_old+i-1))>0.1) printf("Branch current %15s: RHS diff. %-15g || Conv. ref. %-15g || RHS max. %-15g || RHS min. %-15g || Sbranch %-15g \n",node->name, + (*(ckt->CKTrhsOld+i)-*(RHS_copy_se+i-1)), + ( ( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) )*1e-3+1e-9)*7*ckt->CKTsteady_coeff, + *(RHS_max+i-1), + *(RHS_min+i-1), + *(S_diff+i-1) + //( (*(RHS_max+i-1)) - (*(RHS_min+i-1)) ) / ( (*(RHS_max+i-1)) + (*(RHS_min+i-1)) ) / 0.5 + ); + dynamic_test++; /* test on current dynamic consistence */ + } else + { + *(S_old+i-1)=0; + *(S_diff+i-1)=0; + } + } + } + if (dynamic_test==0) { + /* Test for dynamic existence */ + printf("No detectable dynamic on voltages nodes or currents branches. PSS analysis aborted.\n"); + FREE(RHS_copy_se); + FREE(RHS_max); + FREE(RHS_min); + FREE(S_old); + FREE(S_diff); + return(OK); + } + if ((time_err_min_0-time_temp)<0) { + /* Something has gone wrong... */ + printf("Cannot find a minimum for error vector in estimated period. Try to adjust tstab! PSS analysis aborted.\n"); + FREE(RHS_copy_se); + FREE(RHS_max); + FREE(RHS_min); + FREE(S_old); + FREE(S_diff); + return(OK); + } + printf("Global Convergence Error reference: %g.\n", err_conv_ref/dynamic_test); + /*** FREQUENCY ESTIMATION UPDATE ***/ + if ( err_min_0==err || err_min_0==1e+30 ) { + ckt->CKTguessedFreq=(ckt->CKTguessedFreq)+f_proj; + #ifdef STEPDEBUG + printf("Frequency DOWN: est per %g, err min %g, err min 1 %g, err max %g, ntc %ld, err %g, err_last %g\n", time_err_min_0-time_temp,err_min_0,err_min_1,err_max,nextTime_count,err,err_last); + #endif + gf_last_1=gf_last_0; + gf_last_0=ckt->CKTguessedFreq; + } else { + ckt->CKTguessedFreq=1/(time_err_min_0-time_temp); + #ifdef STEPDEBUG + printf("Frequency UP: est per %g, err min %g, err min 1 %g, err max %g, ntc %ld, err %g, err_last %g\n", time_err_min_0-time_temp,err_min_0,err_min_1,err_max,nextTime_count,err,err_last); + #endif + gf_last_1=gf_last_0; + gf_last_0=ckt->CKTguessedFreq; + } + /* Store auxiliary variable of time_temp */ + time_temp_0=time_temp; + /* Next evaluation of shooting will be updated time (time_temp) summed to updated guessed period */ + time_temp=ckt->CKTtime; + /* IMPORTANT! Final time must be updated! Otherwise delta time can be wrongly calculated */ + ckt->CKTfinalTime=time_temp+2/ckt->CKTguessedFreq; + /* Set next the breakpoint */ + CKTsetBreak(ckt, time_temp+1/(ckt->CKTguessedFreq)); + /* Store error history */ + rr_history[shooting_cycle_counter]=err; + gf_history[shooting_cycle_counter]=ckt->CKTguessedFreq; + shooting_cycle_counter++; + printf("Updated guessed frequency: %1.10lg .\n",ckt->CKTguessedFreq); + printf("Next shooting evaluation time is %1.10g and current time is %1.10g.\n",time_temp+1/(ckt->CKTguessedFreq),ckt->CKTtime); + /* shooting exit condition */ + if ( shooting_cycle_counter>ckt->CKTsc_iter || (rr_history[shooting_cycle_counter-1] FE: %15.10g || RR: %15.10g\n",count_5+1,1.0/delta_t/10,rr_history[count_5+1]); /* the very lucky case */ + } else { + printf("%-3d -> FE: %15.10g || RR: %15.10g\n",count_5+1,gf_history[count_5],rr_history[count_5+1]); + } + #endif + /* reuse variables */ + if (rr_history[count_5]CKTsc_iter) { + ckt->CKTguessedFreq=gf_history[shooting_cycle_counter-1]; + printf("\nConvergence reached. Final circuit time is %1.10g s and predicted fundamental frequency is %g Hz.\n",ckt->CKTtime,ckt->CKTguessedFreq); + in_pss=1; /* PERIODIC STEADY STATE NOT REACHED however set the flag */ + } else { + ckt->CKTguessedFreq=gf_history[k-1]; + printf("\nConvergence not reached. However the most near convergence iteration has predicted (iteration %d) a fundamental frequency of %g Hz.\n",k,ckt->CKTguessedFreq); + in_pss=1; /* PERIODIC STEADY STATE REACHED set the flag */ + } + /* Allocates memory for nodes data in PSS */ + psstimes = TMALLOC(double, ckt->CKTpsspoints); + pssvalues = TMALLOC(double, msize*ckt->CKTpsspoints); + pssValues = TMALLOC(double, ckt->CKTpsspoints); + pssfreqs = TMALLOC(double, ckt->CKTharms); + pssmags = TMALLOC(double, ckt->CKTharms); + pssphases = TMALLOC(double, ckt->CKTharms); + pssnmags = TMALLOC(double, ckt->CKTharms); + pssnphases = TMALLOC(double, ckt->CKTharms); + pssResults = TMALLOC(double, msize*ckt->CKTharms); + } + /* restore maximum and minimum error for next search */ + err_min_0=1.0e+30; + err_max=-1.0e+30; + err_0=1.e30; + err_1=-1.0e+30; + tv_03=err_conv_ref; + err_conv_ref=0; + tv_04=dynamic_test; + dynamic_test=0; + rest=shooting_cycle_counter%2; + /* Reset actual RHS reference for next shooting evaluation */ + for(count_6 = 1; count_6 <= msize; count_6++) { + *(RHS_copy_se+count_6-1)= *(ckt->CKTrhsOld+count_6); + } + #ifdef STEPDEBUG + printf("RHS on new shooting cycle: "); + for(count_3 = 1; count_3 <= msize; count_3++) { + printf("%-15g ", *(RHS_copy_se+count_3-1)); + } + printf("\n"); + #endif + if (in_pss!=1) { + for(count_7 = 1; count_7 <= msize; count_7++) + { + /* reset max and min per node or branch on every shooting cycle */ + *(RHS_max+count_7-1)=-1.0e+30; + *(RHS_min+count_7-1)=1.0e+30; + } + } + printf("----------------\n\n"); + } + } else { + /* return on the converged shooting condition */ + if ( AlmostEqualUlps( ckt->CKTtime, time_temp+1/ckt->CKTguessedFreq, 10 ) || (ckt->CKTtime > time_temp+1/ckt->CKTguessedFreq) ) + { + /* restore time */ + ckt->CKTtime=time_temp; + /* Set the breakpoint */ + CKTsetBreak(ckt, time_temp+1/(ckt->CKTguessedFreq)); + pss_cycle_counter++; + if (pss_cycle_counter>1) { + /* End plot in time domain */ + SPfrontEnd->OUTendPlot (((PSSan*)ckt->CKTcurJob)->PSSplot_td); + /* The following line must be placed just before a new OUTpBeginPlot is called */ + error = CKTnames(ckt,&numNames,&nameList); + if (error) return (error); + SPfrontEnd->IFnewUid (ckt, &freqUid, NULL, + "frequency", UID_OTHER, NULL); + error = SPfrontEnd->OUTpBeginPlot (ckt, + ckt->CKTcurJob, + "Frequency Domain Periodic Steady State",freqUid,IF_REAL,numNames,nameList, + IF_REAL,&(((PSSan*)ckt->CKTcurJob)->PSSplot_fd)); + tfree(nameList); + /* ************************* */ + /* Fourier transform on data */ + /* ************************* */ + for(i4 = 1; i4 <= msize; i4++) { + for(k1=0; k1CKTpsspoints; k1++){ + *(pssValues+k1)=*(pssvalues+k1*msize+i4-1); + } + CKTfour(ckt->CKTpsspoints, ckt->CKTharms, &thd, psstimes, pssValues, + ckt->CKTguessedFreq, pssfreqs, pssmags, pssphases, pssnmags, + pssnphases); + for(k1 = 1; k1 <= ckt->CKTharms; k1++) { + *(pssResults+k1-1+(i4-1)*msize)=*(pssmags+k1-1); + } + } + for(k1 = 1; k1 <= ckt->CKTharms; k1++) { + for(i4 = 1; i4 <= msize; i4++) { + *(ckt->CKTrhsOld+i4)=*(pssResults+(k1-1)+(i4-1)*msize); + } + CKTdump( ckt, *(pssfreqs+(k1-1)), ( ((PSSan*)ckt->CKTcurJob)->PSSplot_fd) ); + } + /* End plot in freq domain */ + SPfrontEnd->OUTendPlot (((PSSan*)ckt->CKTcurJob)->PSSplot_fd); + FREE(RHS_copy_se); + FREE(RHS_max); + FREE(RHS_min); + FREE(S_old); + FREE(S_diff); + FREE(pssfreqs); + FREE(psstimes); + FREE(pssValues); + FREE(pssResults); + FREE(pssmags); + FREE(pssphases); + FREE(pssnmags); + FREE(pssnphases); + return(OK); + } + } + } + /* ********************************** */ + /* **** END SHOOTING CODE BLOCK ***** */ + /* ********************************** */ + + if( SPfrontEnd->IFpauseTest() ) { + /* user requested pause... */ + ckt->CKTcurrentAnalysis = DOING_TRAN; + ckt->CKTstat->STATtranTime += SPfrontEnd->IFseconds() - startTime; + ckt->CKTstat->STATtranIter += ckt->CKTstat->STATnumIter - startIters; + ckt->CKTstat->STATtranDecompTime += ckt->CKTstat->STATdecompTime - + startdTime; + ckt->CKTstat->STATtranSolveTime += ckt->CKTstat->STATsolveTime - + startsTime; + ckt->CKTstat->STATtranLoadTime += ckt->CKTstat->STATloadTime - + startlTime; + ckt->CKTstat->STATtranCombTime += ckt->CKTstat->STATcombineTime - + startcTime; + ckt->CKTstat->STATtranSyncTime += ckt->CKTstat->STATsyncTime - + startkTime; + return(E_PAUSE); + } + + /* RESUME */ + resume: + if( (ckt->CKTdelta <= ckt->CKTfinalTime/50) && + (ckt->CKTdelta <= ckt->CKTmaxStep)) { + ; + } else { + #ifdef STEPDEBUG + if(ckt->CKTfinalTime/50CKTmaxStep) { + (void)printf("limited by Tstop/50\n"); + } else { + (void)printf("limited by Tmax == %g\n",ckt->CKTmaxStep); + } + #endif + } + #ifdef HAS_WINDOWS + if (ckt->CKTtime == 0.) + SetAnalyse( "tran init", 0); + else + SetAnalyse( "tran", (int)((ckt->CKTtime * 1000.) / ckt->CKTfinalTime)); + #endif + ckt->CKTdelta = + MIN(ckt->CKTdelta,ckt->CKTmaxStep); + #ifdef XSPICE + /* gtri - begin - wbk - Cut integration order if first timepoint after breakpoint */ + //if(ckt->CKTtime == g_mif_info.breakpoint.last) + if ( AlmostEqualUlps( ckt->CKTtime, g_mif_info.breakpoint.last, 100 ) ) + ckt->CKTorder = 1; + /* gtri - end - wbk - Cut integration order if first timepoint after breakpoint */ + + #endif + + + /* are we at a breakpoint, or indistinguishably close? */ + //if ((ckt->CKTtime == *(ckt->CKTbreaks)) || (*(ckt->CKTbreaks) - + if ( (*(ckt->CKTbreaks)-(ckt->CKTtime) <= ckt->CKTdelmin) ) { + /*if ( AlmostEqualUlps( ckt->CKTtime, *(ckt->CKTbreaks), 100 ) || (*(ckt->CKTbreaks) - + * (ckt->CKTtime) <= ckt->CKTdelmin)) {*/ + /* first timepoint after a breakpoint - cut integration order */ + /* and limit timestep to .1 times minimum of time to next breakpoint, + * and previous timestep + */ + ckt->CKTorder = 1; + if( (ckt->CKTdelta >.1* ckt->CKTsaveDelta) || + (ckt->CKTdelta > .1*(*(ckt->CKTbreaks+1)-*(ckt->CKTbreaks))) ) { + if(ckt->CKTsaveDelta < (*(ckt->CKTbreaks+1)-*(ckt->CKTbreaks))) { + #ifdef STEPDEBUG + (void)printf("limited by pre-breakpoint delta (saveDelta: %1.10g, nxt_breakpt: %1.10g, curr_breakpt: %1.10g and CKTtime: %1.10g\n", + ckt->CKTsaveDelta, *(ckt->CKTbreaks+1), *(ckt->CKTbreaks), ckt->CKTtime); + #endif + } else { + #ifdef STEPDEBUG + (void)printf("limited by next breakpoint\n"); + (void)printf("(saveDelta: %1.10g, Delta: %1.10g, CKTtime: %1.10g and delmin: %1.10g\n",ckt->CKTsaveDelta,ckt->CKTdelta,ckt->CKTtime,ckt->CKTdelmin); + #endif + } + } + + if ( (*(ckt->CKTbreaks+1)-*(ckt->CKTbreaks)) == 0 ) { + ckt->CKTdelta = ckt->CKTdelmin; + } else { + ckt->CKTdelta = MIN(ckt->CKTdelta, .1 * MIN(ckt->CKTsaveDelta, + *(ckt->CKTbreaks+1)-*(ckt->CKTbreaks))); + } + + if(firsttime) { + ckt->CKTdelta /= 10; + #ifdef STEPDEBUG + (void)printf("delta cut for initial timepoint\n"); + #endif + } + + #ifdef XSPICE + } + + /* gtri - begin - wbk - Add Breakpoint stuff */ + + if(ckt->CKTtime + ckt->CKTdelta >= g_mif_info.breakpoint.current) { + /* If next time > temporary breakpoint, force it to the breakpoint */ + /* And mark that timestep was set by temporary breakpoint */ + ckt->CKTsaveDelta = ckt->CKTdelta; + ckt->CKTdelta = g_mif_info.breakpoint.current - ckt->CKTtime; + g_mif_info.breakpoint.last = ckt->CKTtime + ckt->CKTdelta; + } + else { + /* Else, mark that timestep was not set by temporary breakpoint */ + g_mif_info.breakpoint.last = 1.0e30; + } + + /* gtri - end - wbk - Add Breakpoint stuff */ + + /* gtri - begin - wbk - Modify Breakpoint stuff */ + /* Throw out any permanent breakpoint times <= current time */ + while(1) { + printf(" brk_pt: %g ckt_time: %g ckt_min_break: %g\n",*(ckt->CKTbreaks), ckt->CKTtime, ckt->CKTminBreak); + if(AlmostEqualUlps(*(ckt->CKTbreaks),ckt->CKTtime, 100) || *(ckt->CKTbreaks) <= (ckt->CKTtime + ckt->CKTminBreak)) { + printf("throwing out permanent breakpoint times <= current time (brk pt: %g)\n",*(ckt->CKTbreaks)); + printf("ckt_time: %g ckt_min_break: %g\n",ckt->CKTtime, ckt->CKTminBreak); + CKTclrBreak(ckt); + } + else + break; + } + /* Force the breakpoint if appropriate */ + if((ckt->CKTtime + ckt->CKTdelta) > *(ckt->CKTbreaks)) { + ckt->CKTbreak = 1; + ckt->CKTsaveDelta = ckt->CKTdelta; + ckt->CKTdelta = *(ckt->CKTbreaks) - ckt->CKTtime; + } + + /* gtri - end - wbk - Modify Breakpoint stuff */ + #else /* !XSPICE */ + + /* don't want to get below delmin for no reason */ + ckt->CKTdelta = MAX(ckt->CKTdelta, ckt->CKTdelmin*2.0); + } + else if(ckt->CKTtime + ckt->CKTdelta >= *(ckt->CKTbreaks)) { + ckt->CKTsaveDelta = ckt->CKTdelta; + ckt->CKTdelta = *(ckt->CKTbreaks) - ckt->CKTtime; + /*(void)printf("delta cut to %g to hit breakpoint\n",ckt->CKTdelta);*/ + fflush(stdout); + ckt->CKTbreak = 1; /* why? the current pt. is not a bkpt. */ + } + #ifdef CLUSTER + if(!CLUsync(ckt->CKTtime,&ckt->CKTdelta,0)){ + printf("Sync error!\n"); + exit(0); + } + #endif + #ifdef PARALLEL_ARCH + DGOP_( &type, &(ckt->CKTdelta), &length, "min" ); + #endif /* PARALLEL_ARCH */ + + #endif /* XSPICE */ + + #ifdef XSPICE + /* gtri - begin - wbk - Do event solution */ + + if(ckt->evt->counts.num_insts > 0) { + + /* if time = 0 and op_alternate was specified as false during */ + /* dcop analysis, call any changed instances to let them */ + /* post their outputs with their associated delays */ + if((ckt->CKTtime == 0.0) && (! ckt->evt->options.op_alternate)) + EVTiter(ckt); + + /* while there are events on the queue with event time <= next */ + /* projected analog time, process them */ + while((g_mif_info.circuit.evt_step = EVTnext_time(ckt)) + <= (ckt->CKTtime + ckt->CKTdelta)) { + + /* Initialize temp analog bkpt to infinity */ + g_mif_info.breakpoint.current = 1e30; + + /* Pull items off queue and process them */ + EVTdequeue(ckt, g_mif_info.circuit.evt_step); + EVTiter(ckt); + + /* If any instances have forced an earlier */ + /* next analog time, cut the delta */ + if(*(ckt->CKTbreaks) < g_mif_info.breakpoint.current) + if(*(ckt->CKTbreaks) > (ckt->CKTtime + ckt->CKTminBreak)) + g_mif_info.breakpoint.current = *(ckt->CKTbreaks); + if(g_mif_info.breakpoint.current < (ckt->CKTtime + ckt->CKTdelta)) { + /* Breakpoint must be > last accepted timepoint */ + /* and >= current event time */ + if(g_mif_info.breakpoint.current > (ckt->CKTtime + ckt->CKTminBreak) + && (g_mif_info.breakpoint.current >= g_mif_info.circuit.evt_step)) { + ckt->CKTsaveDelta = ckt->CKTdelta; + ckt->CKTdelta = g_mif_info.breakpoint.current - ckt->CKTtime; + g_mif_info.breakpoint.last = ckt->CKTtime + ckt->CKTdelta; + } + } + + } /* end while next event time <= next analog time */ + } /* end if there are event instances */ + + /* gtri - end - wbk - Do event solution */ + #endif + for(i5=5;i5>=0;i5--) { + ckt->CKTdeltaOld[i5+1]=ckt->CKTdeltaOld[i5]; + } + ckt->CKTdeltaOld[0]=ckt->CKTdelta; + + temp = ckt->CKTstates[ckt->CKTmaxOrder+1]; + for(i5=ckt->CKTmaxOrder;i5>=0;i5--) { + ckt->CKTstates[i5+1] = ckt->CKTstates[i5]; + } + ckt->CKTstates[0] = temp; + + /* 600 */ + while (1) { + #ifdef CLUSTER + redostep = 1; + #endif + #ifdef XSPICE + /* gtri - add - wbk - 4/17/91 - Fix Berkeley bug */ + /* This is needed here to allow CAPask to output currents */ + /* during Transient analysis. A grep for CKTcurrentAnalysis */ + /* indicates that it should not hurt anything else ... */ + + ckt->CKTcurrentAnalysis = DOING_TRAN; + + /* gtri - end - wbk - 4/17/91 - Fix Berkeley bug */ + #endif + olddelta=ckt->CKTdelta; + /* time abort? */ + + /* ************************************ */ + /* ********** CKTtime update ********** */ + /* ************************************ */ + /* delta manipulation */ + if (!in_stabilization && !in_pss){ + if ( (ckt->CKTtime-(time_temp+1/ckt->CKTguessedFreq)<1/ckt->CKTguessedFreq/10) && (ckt->CKTbreak==0) ) + { + if ( !(flag_tu_1) ) flag_tu_2=ckt->CKTdelta; /* store previous delta */ + if ( (ckt->CKTtime-(time_temp+1/ckt->CKTguessedFreq)<1/ckt->CKTguessedFreq/1.0e5) && (ckt->CKTbreak==0) ) + { + if ( (ckt->CKTtime-(time_temp+1/ckt->CKTguessedFreq)<1/ckt->CKTguessedFreq/1.0e7) && (ckt->CKTbreak==0) ) + { + if (rr_history[shooting_cycle_counter-1]CKTdelta=1/ckt->CKTguessedFreq/1.0e5; /* get closer to accurate solution? */ + } else { + ckt->CKTdelta=1/ckt->CKTguessedFreq/1.0e4; + } + } else + { + ckt->CKTdelta=1/ckt->CKTguessedFreq/1.0e1; + } + } else + { + ckt->CKTdelta=1/ckt->CKTguessedFreq/0.25e1; + } + flag_tu_1=1; + } else + { + if (flag_tu_1) { + ckt->CKTdelta=flag_tu_2; /* restore prevoius delta */ + flag_tu_1=0; + } + } + } + if ( in_pss ) ckt->CKTdelta=1/ckt->CKTguessedFreq/((ckt->CKTpsspoints-1)); /* fixed delta in PSS */ + /* ************************************ */ + /* ******* END CKTtime update ********* */ + /* ************************************ */ + + ckt->CKTtime += ckt->CKTdelta; + + #ifdef CLUSTER + CLUinput(ckt); + #endif + ckt->CKTdeltaOld[0]=ckt->CKTdelta; + NIcomCof(ckt); + #ifdef PREDICTOR + error = NIpred(ckt); + #endif /* PREDICTOR */ + save_mode = ckt->CKTmode; + save_order = ckt->CKTorder; + #ifdef XSPICE + /* gtri - begin - wbk - Add Breakpoint stuff */ + + /* Initialize temporary breakpoint to infinity */ + g_mif_info.breakpoint.current = 1.0e30; + + /* gtri - end - wbk - Add Breakpoint stuff */ + + + /* gtri - begin - wbk - add convergence problem reporting flags */ + /* delta is forced to equal delmin on last attempt near line 650 */ + if(ckt->CKTdelta <= ckt->CKTdelmin) + ckt->enh->conv_debug.last_NIiter_call = MIF_TRUE; + else + ckt->enh->conv_debug.last_NIiter_call = MIF_FALSE; + /* gtri - begin - wbk - add convergence problem reporting flags */ + + + + /* gtri - begin - wbk - Call all hybrids */ + + /* gtri - begin - wbk - Set evt_step */ + + if(ckt->evt->counts.num_insts > 0) { + g_mif_info.circuit.evt_step = ckt->CKTtime; + } + /* gtri - end - wbk - Set evt_step */ + #endif + + converged = NIiter(ckt,ckt->CKTtranMaxIter); + + #ifdef XSPICE + if(ckt->evt->counts.num_insts > 0) + { + g_mif_info.circuit.evt_step = ckt->CKTtime; + EVTcall_hybrids(ckt); + } + /* gtri - end - wbk - Call all hybrids */ + + #endif + ckt->CKTstat->STATtimePts ++; + ckt->CKTmode = (ckt->CKTmode&MODEUIC)|MODETRAN | MODEINITPRED; + if(firsttime) { + for(i=0;iCKTnumStates;i++) { + *(ckt->CKTstate2+i) = *(ckt->CKTstate1+i); + *(ckt->CKTstate3+i) = *(ckt->CKTstate1+i); + } + } + /* txl, cpl addition */ + if (converged == 1111) { + return(converged); + } + + if(converged != 0) { + #ifndef CLUSTER + ckt->CKTtime = ckt->CKTtime -ckt->CKTdelta; + ckt->CKTstat->STATrejected ++; + #endif + ckt->CKTdelta = ckt->CKTdelta/8; + /*printf("delta cut to %g for non-convergance\n",ckt->CKTdelta);*/ + #ifdef STEPDEBUG + (void)printf("delta cut to %g for non-convergance\n",ckt->CKTdelta); + fflush(stdout); + #endif + if(firsttime) { + ckt->CKTmode = (ckt->CKTmode&MODEUIC)|MODETRAN | MODEINITTRAN; + } + ckt->CKTorder = 1; + + #ifdef XSPICE + /* gtri - begin - wbk - Add Breakpoint stuff */ + + /* Force backup if temporary breakpoint is < current time */ + } else if(g_mif_info.breakpoint.current < ckt->CKTtime) { + ckt->CKTsaveDelta = ckt->CKTdelta; + ckt->CKTtime -= ckt->CKTdelta; + ckt->CKTdelta = g_mif_info.breakpoint.current - ckt->CKTtime; + g_mif_info.breakpoint.last = ckt->CKTtime + ckt->CKTdelta; + + if(firsttime) { + ckt->CKTmode = (ckt->CKTmode&MODEUIC)|MODETRAN | MODEINITTRAN; + } + ckt->CKTorder = 1; + + /* gtri - end - wbk - Add Breakpoint stuff */ + #endif + + } else { + if (firsttime) { + #ifdef WANT_SENSE2 + if(ckt->CKTsenInfo && (ckt->CKTsenInfo->SENmode & TRANSEN)){ + save1 = ckt->CKTmode; + save2 = ckt->CKTorder; + ckt->CKTmode = save_mode; + ckt->CKTorder = save_order; + if(error = CKTsenDCtran(ckt)) return(error); + ckt->CKTmode = save1; + ckt->CKTorder = save2; + } + #endif + firsttime=0; + #ifndef CLUSTER + goto nextTime; /* no check on + * first time point + */ + #else + redostep = 0; + goto chkStep; + #endif + } + new = ckt->CKTdelta; + error = CKTtrunc(ckt,&new); + if(error) { + ckt->CKTcurrentAnalysis = DOING_TRAN; + ckt->CKTstat->STATtranTime += + SPfrontEnd->IFseconds() - startTime; + ckt->CKTstat->STATtranIter += + ckt->CKTstat->STATnumIter - startIters; + ckt->CKTstat->STATtranDecompTime += ckt->CKTstat->STATdecompTime + - startdTime; + ckt->CKTstat->STATtranSolveTime += ckt->CKTstat->STATsolveTime + - startsTime; + ckt->CKTstat->STATtranLoadTime += ckt->CKTstat->STATloadTime + - startlTime; + ckt->CKTstat->STATtranCombTime += ckt->CKTstat->STATcombineTime + - startcTime; + ckt->CKTstat->STATtranSyncTime += ckt->CKTstat->STATsyncTime + - startkTime; + return(error); + } + if(new>.9 * ckt->CKTdelta) { + if(ckt->CKTorder == 1) { + new = ckt->CKTdelta; + ckt->CKTorder = 2; + error = CKTtrunc(ckt,&new); + if(error) { + ckt->CKTcurrentAnalysis = DOING_TRAN; + ckt->CKTstat->STATtranTime += + SPfrontEnd->IFseconds() - startTime; + ckt->CKTstat->STATtranIter += + ckt->CKTstat->STATnumIter - startIters; + ckt->CKTstat->STATtranDecompTime += + ckt->CKTstat->STATdecompTime - startdTime; + ckt->CKTstat->STATtranSolveTime += + ckt->CKTstat->STATsolveTime - startsTime; + ckt->CKTstat->STATtranLoadTime += + ckt->CKTstat->STATloadTime - startlTime; + ckt->CKTstat->STATtranCombTime += + ckt->CKTstat->STATcombineTime - startcTime; + ckt->CKTstat->STATtranSyncTime += + ckt->CKTstat->STATsyncTime - startkTime; + return(error); + } + if(new <= 1.05 * ckt->CKTdelta) { + ckt->CKTorder = 1; + } + } + /* time point OK - 630*/ + ckt->CKTdelta = new; + #ifdef NDEV + /* show a time process indicator, by Gong Ding, gdiso@ustc.edu */ + if(ckt->CKTtime/ckt->CKTfinalTime*100<10.0) + printf("%%%3.2lf\b\b\b\b\b",ckt->CKTtime/ckt->CKTfinalTime*100); + else if(ckt->CKTtime/ckt->CKTfinalTime*100<100.0) + printf("%%%4.2lf\b\b\b\b\b\b",ckt->CKTtime/ckt->CKTfinalTime*100); + else + printf("%%%5.2lf\b\b\b\b\b\b\b",ckt->CKTtime/ckt->CKTfinalTime*100); + fflush(stdout); + #endif + + #ifdef STEPDEBUG + (void)printf( + "delta set to truncation error result: %g. Point accepted at CKTtime: %g\n", + ckt->CKTdelta,ckt->CKTtime); + fflush(stdout); + #endif + #ifdef WANT_SENSE2 + if(ckt->CKTsenInfo && (ckt->CKTsenInfo->SENmode & TRANSEN)){ + save1 = ckt->CKTmode; + save2 = ckt->CKTorder; + ckt->CKTmode = save_mode; + ckt->CKTorder = save_order; + if(error = CKTsenDCtran(ckt)) return(error); + ckt->CKTmode = save1; + ckt->CKTorder = save2; + } + #endif + #ifndef CLUSTER + /* go to 650 - trapezoidal */ + goto nextTime; + #else + redostep = 0; + goto chkStep; + #endif + } else { + #ifndef CLUSTER + ckt->CKTtime = ckt->CKTtime -ckt->CKTdelta; + ckt->CKTstat->STATrejected ++; + #endif + ckt->CKTdelta = new; + #ifdef STEPDEBUG + (void)printf( + "delta set to truncation error result:point rejected\n"); + #endif + } + } + #ifdef PARALLEL_ARCH + DGOP_( &type, &(ckt->CKTdelta), &length, "min" ); + #endif /* PARALLEL_ARCH */ + + if (ckt->CKTdelta <= ckt->CKTdelmin) { + if (olddelta > ckt->CKTdelmin) { + ckt->CKTdelta = ckt->CKTdelmin; + /*#ifdef STEPDEBUG*/ + (void)printf("delta at delmin\n"); + /*#endif*/ + } else { + ckt->CKTcurrentAnalysis = DOING_TRAN; + ckt->CKTstat->STATtranTime += + SPfrontEnd->IFseconds() - startTime; + ckt->CKTstat->STATtranIter += + ckt->CKTstat->STATnumIter - startIters; + ckt->CKTstat->STATtranDecompTime += + ckt->CKTstat->STATdecompTime - startdTime; + ckt->CKTstat->STATtranSolveTime += + ckt->CKTstat->STATsolveTime - startsTime; + ckt->CKTstat->STATtranLoadTime += + ckt->CKTstat->STATloadTime - startlTime; + ckt->CKTstat->STATtranCombTime += + ckt->CKTstat->STATcombineTime - startcTime; + ckt->CKTstat->STATtranSyncTime += + ckt->CKTstat->STATsyncTime - startkTime; + errMsg = CKTtrouble(ckt, "Timestep too small"); + return(E_TIMESTEP); + } + } + #ifdef XSPICE + /* gtri - begin - wbk - Do event backup */ + + if(ckt->evt->counts.num_insts > 0) + EVTbackup(ckt, ckt->CKTtime + ckt->CKTdelta); + + /* gtri - end - wbk - Do event backup */ + #endif + #ifdef CLUSTER + chkStep: + if(CLUsync(ckt->CKTtime,&ckt->CKTdelta,redostep)){ + goto nextTime; + } else { + ckt->CKTtime -= olddelta; + ckt->CKTstat->STATrejected ++; + } + #endif + } + /* NOTREACHED */ +} + +int +CKTfour(int ndata, /* number of entries in the Time and + Value arrays */ + int numFreq, /* number of harmonics to calculate */ + double *thd, /* total harmonic distortion (percent) + to be returned */ + double *Time, /* times at which the voltage/current + values were measured*/ + double *Value, /* voltage or current vector whose + transform is desired */ + double FundFreq, /* the fundamental frequency of the + analysis */ + double *Freq, /* the frequency value of the various + harmonics */ + double *Mag, /* the Magnitude of the fourier + transform */ + double *Phase, /* the Phase of the fourier transform */ + double *nMag, /* the normalized magnitude of the + transform: nMag(fund)=1*/ + double *nPhase) /* the normalized phase of the + transform: Nphase(fund)=0 */ +{ + /* Note: we can consider these as a set of arrays. The sizes are: + * Time[ndata], Value[ndata], Freq[numFreq], Mag[numfreq], + * Phase[numfreq], nMag[numfreq], nPhase[numfreq] + * + * The arrays must all be allocated by the caller. + * The Time and Value array must be reasonably distributed over at + * least one full period of the fundamental Frequency for the + * fourier transform to be useful. The function will take the + * last period of the frequency as data for the transform. + * + * We are assuming that the caller has provided exactly one period + * of the fundamental frequency. */ + int i; + int j; + double tmp; + + NG_IGNORE(Time); + + /* clear output/computation arrays */ + + for(i=0;i1) *thd += nMag[i]*nMag[i]; + } + *thd = 100*sqrt(*thd); + return(OK); +} diff --git a/src/spicelib/analysis/pssaskq.c b/src/spicelib/analysis/pssaskq.c new file mode 100755 index 000000000..f58b157f0 --- /dev/null +++ b/src/spicelib/analysis/pssaskq.c @@ -0,0 +1,52 @@ +/********** +Author: 2010-05 Stefano Perticaroli ``spertica'' +**********/ + +#include "ngspice.h" +#include "ifsim.h" +#include "iferrmsg.h" +#include "cktdefs.h" +#include "pssdefs.h" + +/* ARGSUSED */ +int +PSSaskQuest(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value) +{ + NG_IGNORE(ckt); + + switch(which) { + + case GUESSED_FREQ: + value->rValue = ((PSSan *)anal)->PSSguessedFreq; + break; + case OSC_NODE: + value->nValue = ((PSSan *)anal)->PSSoscNode; + break; + case STAB_TIME: + value->rValue = ((PSSan *)anal)->PSSstabTime; + break; + case PSS_UIC: + if(((PSSan *)anal)->PSSmode & MODEUIC) { + value->iValue = 1; + } else { + value->iValue = 0; + } + break; + case PSS_POINTS: + value->iValue = ((PSSan *)anal)->PSSpoints; + break; + case PSS_HARMS: + value->iValue = ((PSSan *)anal)->PSSharms; + break; + case SC_ITER: + value->iValue = ((PSSan *)anal)->sc_iter; + break; + case STEADY_COEFF: + value->rValue = ((PSSan *)anal)->steady_coeff; + break; + + default: + return(E_BADPARM); + } + return(OK); +} diff --git a/src/spicelib/analysis/pssinit.c b/src/spicelib/analysis/pssinit.c new file mode 100755 index 000000000..307f6b021 --- /dev/null +++ b/src/spicelib/analysis/pssinit.c @@ -0,0 +1,33 @@ +/********** +Author: 2010-05 Stefano Perticaroli ``spertica'' +**********/ + +#include "ngspice.h" +#include "cktdefs.h" +#include "trandefs.h" +#include "pssdefs.h" +#include "iferrmsg.h" + +int PSSinit(CKTcircuit *ckt, JOB *job) +{ + /* Final time depends on stabilization time requested for PSS + and on at least one more oscillation period */ + ckt->CKTfinalTime = ((PSSan*)job)->PSSstabTime + 2/((PSSan*)job)->PSSguessedFreq; + /* Step is chosen empirically to be 1% of PSSguessedFreq */ + ckt->CKTstep = 0.01 * (1/((PSSan*)job)->PSSguessedFreq); + /* Init time should be always zero */ + ckt->CKTinitTime = 0; + /* MaxStep should not exceed Nyquist criterion */ + ckt->CKTmaxStep = 0.5*(1/((PSSan*)job)->PSSguessedFreq); + ckt->CKTdelmin = 1e-9*ckt->CKTmaxStep; + ckt->CKTmode = ((PSSan*)job)->PSSmode; + /* modified CKTdefs.h for the following - 100609 */ + ckt->CKTstabTime = ((PSSan*)job)->PSSstabTime; + ckt->CKTguessedFreq = ((PSSan*)job)->PSSguessedFreq; + ckt->CKTharms = ((PSSan*)job)->PSSharms; + ckt->CKTpsspoints = ((PSSan*)job)->PSSpoints; + ckt->CKTsc_iter = ((PSSan*)job)->sc_iter; + ckt->CKTsteady_coeff = ((PSSan*)job)->steady_coeff; + + return OK; +} diff --git a/src/spicelib/analysis/psssetp.c b/src/spicelib/analysis/psssetp.c new file mode 100755 index 000000000..5bc12fac4 --- /dev/null +++ b/src/spicelib/analysis/psssetp.c @@ -0,0 +1,81 @@ +/********** +Author: 2010-05 Stefano Perticaroli ``spertica'' +**********/ + +#include "ngspice.h" +#include "ifsim.h" +#include "iferrmsg.h" +#include "cktdefs.h" +#include "pssdefs.h" + +#include "analysis.h" + +/* ARGSUSED */ +int +PSSsetParm(CKTcircuit *ckt, JOB *anal, int which, IFvalue *value) +{ + NG_IGNORE(ckt); + + switch(which) { + + case GUESSED_FREQ: + ((PSSan *)anal)->PSSguessedFreq = value->rValue; + break; + case OSC_NODE: + ((PSSan *)anal)->PSSoscNode = value->nValue; + break; + case STAB_TIME: + ((PSSan *)anal)->PSSstabTime = value->rValue; + break; + case PSS_POINTS: + ((PSSan *)anal)->PSSpoints = value->iValue; + break; + case PSS_HARMS: + ((PSSan *)anal)->PSSharms = value->iValue; + break; + case PSS_UIC: + if(value->iValue) { + ((PSSan *)anal)->PSSmode |= MODEUIC; + } + break; + case SC_ITER: + ((PSSan *)anal)->sc_iter = value->iValue; + break; + case STEADY_COEFF: + ((PSSan *)anal)->steady_coeff = value->rValue; + break; + + default: + return(E_BADPARM); + } + return(OK); +} + + +static IFparm PSSparms[] = { + { "fguess", GUESSED_FREQ, IF_SET|IF_REAL, "guessed frequency" }, + { "oscnode", OSC_NODE, IF_SET|IF_STRING, "oscillation node" }, + { "stabtime", STAB_TIME, IF_SET|IF_REAL, "stabilization time" }, + { "points", PSS_POINTS, IF_SET|IF_INTEGER, "pick equispaced number of time points in PSS" }, + { "harmonics", PSS_HARMS, IF_SET|IF_INTEGER, "consider only given number of harmonics in PSS from DC" }, + { "uic", PSS_UIC, IF_SET|IF_INTEGER, "use initial conditions (1 true - 0 false)" }, + { "sc_iter", SC_ITER, IF_SET|IF_INTEGER, "maxmimum number of shooting cycle iterations" }, + { "steady_coeff", STEADY_COEFF, IF_SET|IF_INTEGER, "set steady coefficient for convergence test" } +}; + +SPICEanalysis PSSinfo = { + { + "PSS", + "Periodic Steady State analysis", + + sizeof(PSSparms)/sizeof(IFparm), + PSSparms + }, + sizeof(PSSan), + TIMEDOMAIN, + 1, + PSSsetParm, + PSSaskQuest, + PSSinit, + DCpss +}; diff --git a/src/spicelib/parser/inp2dot.c b/src/spicelib/parser/inp2dot.c index 71aed63cb..2a268c143 100644 --- a/src/spicelib/parser/inp2dot.c +++ b/src/spicelib/parser/inp2dot.c @@ -627,6 +627,73 @@ dot_sens2(char *line, CKTcircuit *ckt, INPtables *tab, card *current, } #endif +/*SP: Steady State Analyis */ +static int +dot_pss(char *line, void *ckt, INPtables *tab, card *current, + void *task, void *gnode, JOB *foo) +{ + int error; /* error code temporary */ + IFvalue ptemp; /* a value structure to package resistance into */ + IFvalue *parm; /* a pointer to a value struct for function returns */ + char *nname; /* the oscNode name */ + CKTnode *nnode; /* the oscNode node */ + int which; /* which analysis we are performing */ + int i; /* generic loop variable */ + char *word; /* something to stick a word of input into */ + + NG_IGNORE(gnode); + + /* .pss Fguess StabTime OscNode */ + which = -1; + for (i = 0; i < ft_sim->numAnalyses; i++) { + if (strcmp(ft_sim->analyses[i]->name, "PSS") == 0) { + which = i; + break; + } + } + if (which == -1) { + LITERR("Periodic steady state analysis unsupported.\n"); + return (0); + } + IFC(newAnalysis, (ckt, which, "Periodic steady state Analysis", &foo, task)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* Fguess */ + GCA(INPapName, (ckt, which, foo, "fguess", parm)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* StabTime */ + GCA(INPapName, (ckt, which, foo, "stabtime", parm)); + + INPgetNetTok(&line, &nname, 0); + INPtermInsert(ckt, &nname, tab, &nnode); + ptemp.nValue = nnode; + GCA(INPapName, (ckt, which, foo, "oscnode", &ptemp)) /* OscNode given as string */ + + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* PSS points */ + GCA(INPapName, (ckt, which, foo, "points", parm)); + + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* PSS harmonics */ + GCA(INPapName, (ckt, which, foo, "harmonics", parm)); + + if (*line) { + INPgetTok(&line, &word, 1); /* uic? */ + if (strcmp(word, "uic") == 0) { + ptemp.iValue = 1; + GCA(INPapName, (ckt, which, foo, "uic", &ptemp)); + } else { + LITERR(" Error: unknown parameter on .pss - ignored\n"); + } + } + + parm = INPgetValue(ckt, &line, IF_INTEGER, tab); /* SC iterations */ + GCA(INPapName, (ckt, which, foo, "sc_iter", parm)); + + parm = INPgetValue(ckt, &line, IF_REAL, tab); /* Steady coefficient */ + GCA(INPapName, (ckt, which, foo, "steady_coeff", parm)); + + return (0); +} +/* SP */ + static int dot_options(char *line, CKTcircuit *ckt, INPtables *tab, card *current, TSKtask *task, CKTnode *gnode, JOB *foo) @@ -706,6 +773,11 @@ INP2dot(CKTcircuit *ckt, INPtables *tab, card *current, TSKtask *task, CKTnode * } else if ((strcmp(token, ".tran") == 0)) { rtn = dot_tran(line, ckt, tab, current, task, gnode, foo); goto quit; +/* SP: Steady State Analysis */ + } else if ((strcmp(token, ".pss") == 0)) { + rtn = dot_pss(line, ckt, tab, current, task, gnode, foo); + goto quit; +/* SP */ } else if ((strcmp(token, ".subckt") == 0) || (strcmp(token, ".ends") == 0)) { /* not yet implemented - warn & ignore */