/********** Author: 2010-05 Stefano Perticaroli ``spertica'' First Review: 2012-02 Francesco Lannutti and Stefano Perticaroli ``spertica'' Second Review: 2012-10 Stefano Perticaroli ``spertica'' and Francesco Lannutti **********/ /* Include files for the PSS analysis */ #include "ngspice/ngspice.h" #include "ngspice/cktdefs.h" #include "cktaccept.h" #include "ngspice/pssdefs.h" #include "ngspice/sperror.h" #include "ngspice/fteext.h" #ifdef XSPICE /* gtri - add - wbk - Add headers */ #include "ngspice/miftypes.h" #include "ngspice/evt.h" #include "ngspice/mif.h" #include "ngspice/evtproto.h" #include "ngspice/ipctiein.h" /* gtri - end - wbk - Add headers */ #endif #ifdef CLUSTER #include "ngspice/cluster.h" #endif #define INIT_STATS() \ do { \ startTime = SPfrontEnd->IFseconds(); \ startIters = ckt->CKTstat->STATnumIter; \ startdTime = ckt->CKTstat->STATdecompTime; \ startsTime = ckt->CKTstat->STATsolveTime; \ startlTime = ckt->CKTstat->STATloadTime; \ startkTime = ckt->CKTstat->STATsyncTime; \ } while(0) #define UPDATE_STATS(analysis) \ do { \ ckt->CKTcurrentAnalysis = analysis; \ 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->STATtranSyncTime += ckt->CKTstat->STATsyncTime - startkTime; \ } while(0) /* Define some useful macro */ #define HISTORY 1024 #define GF_LAST 313 static int DFT(long int, int, double *, double *, double *, double, double *, double *, double *, double *, double *); int DCpss(CKTcircuit *ckt, int restart) /* forced restart flag */ { PSSan *job = (PSSan *) ckt->CKTcurJob; int i; double olddelta; double delta; double newdelta; double *temp; double startdTime; double startsTime; double startlTime; double startkTime; double startTime; int startIters; int converged; int firsttime; int error; int save_order; long save_mode; IFuid timeUid; IFuid *nameList; int numNames; double maxstepsize = 0.0; int ltra_num; CKTnode *node; #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 */ int j, oscnNode; IFuid freqUid; enum {STABILIZATION, SHOOTING, PSS} pss_state = STABILIZATION; double err = 0, predsum = 0 ; double time_temp = 0, gf_history [HISTORY], rr_history [HISTORY], predsum_history [HISTORY], nextstep ; int msize, shooting_cycle_counter = 0; double *RHS_copy_se, *RHS_copy_der, *RHS_derivative, *pred, err_0 = HUGE_VAL ; double time_err_min_1 = 0, time_err_min_0 = 0, err_min_0 = HUGE_VAL, err_min_1 = 0 ; double err_1 = 0, err_max = HUGE_VAL ; int pss_points_cycle = 0, dynamic_test = 0 ; double gf_last_0 = HUGE_VAL, gf_last_1 = GF_LAST ; double thd = 0 ; double *psstimes, *pssvalues; double *RHS_max, *RHS_min, *err_conv ; /* Francesco Lannutti's MOD */ /* Stuff needed by frequency estimation reiteration, based on the DFT result */ int position; double max_freq; /* Print some useful information */ fprintf (stderr, "Periodic Steady State Analysis Started\n\n") ; fprintf (stderr, "PSS Guessed Frequency %g\n", ckt->CKTguessedFreq) ; fprintf (stderr, "PSS Points %ld\n", ckt->CKTpsspoints) ; fprintf (stderr, "PSS Harmonics number %d\n", ckt->CKTharms) ; fprintf (stderr, "PSS Steady Coefficient %g\n", ckt->CKTsteady_coeff) ; fprintf (stderr, "PSS sc_iter %d\n", ckt->CKTsc_iter) ; fprintf (stderr, "PSS Stabilization Time %g\n", ckt->CKTstabTime) ; oscnNode = job->PSSoscNode->number ; /* Variables and memory initialization */ for (i = 0 ; i < HISTORY ; i++) { rr_history [i] = 0.0 ; gf_history [i] = 0.0 ; } msize = SMPmatSize (ckt->CKTmatrix) ; RHS_copy_se = TMALLOC (double, msize) ; /* Set the current RHS reference for next Shooting Evaluation */ RHS_copy_der = TMALLOC (double, msize) ; /* Used to compute current Derivative */ RHS_derivative = TMALLOC (double, msize) ; pred = TMALLOC (double, msize) ; RHS_max = TMALLOC (double, msize) ; RHS_min = TMALLOC (double, msize) ; err_conv = TMALLOC (double, msize) ; for (i = 0 ; i < msize ; i++) { RHS_copy_se [i] = 0.0 ; RHS_copy_der [i] = 0.0 ; RHS_derivative [i] = 0.0 ; pred [i] = 0.0 ; } psstimes = TMALLOC (double, ckt->CKTpsspoints + 1) ; pssvalues = TMALLOC (double, msize * (ckt->CKTpsspoints + 1)) ; for (i = 0 ; i < ckt->CKTpsspoints + 1 ; i++) psstimes [i] = 0.0 ; for (i = 0 ; i < msize * (ckt->CKTpsspoints + 1) ; i++) pssvalues [i] = 0.0 ; /* Delta timestep and circuit time setup */ delta = ckt->CKTstep ; ckt->CKTtime = ckt->CKTinitTime ; ckt->CKTfinalTime = ckt->CKTstabTime ; /* Starting PSS Algorithm, based on Transient Analysis */ if(restart || ckt->CKTtime == 0) { delta = MIN (1 / ckt->CKTguessedFreq / 100, ckt->CKTstep) ; #ifdef STEPDEBUG fprintf (stderr, "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 */ /* Breakpoints initialization */ if(ckt->CKTbreaks) FREE(ckt->CKTbreaks); ckt->CKTbreaks = TMALLOC(double, 2); if(ckt->CKTbreaks == NULL) return(E_NOMEM); ckt->CKTbreaks[0] = 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 /* Minimum Breakpoint Setup */ 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 /* Time Domain plot start and prepared to be filled in later */ 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 Analysis", timeUid, IF_REAL, numNames, nameList, IF_REAL, &(job->PSSplot_td)); tfree(nameList); if(error) return(error); /* Time initialization for Transient Analysis */ 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 /* Looking for a working Operating Point */ 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); /* CKTnode *node; double new, old, tol; int i=1; fprintf(stdout,"\nTransient solution failed -\n\n"); fprintf(stdout,"Last Node Voltages\n"); fprintf(stdout,"------------------\n\n"); fprintf(stdout,"%-30s %20s %20s\n", "Node", "Last Voltage", "Previous Iter"); fprintf(stdout,"%-30s %20s %20s\n", "----", "------------", "-------------"); for(node=ckt->CKTnodes->next;node;node=node->next) { if (strstr(node->name, "#branch") || !strstr(node->name, "#")) { new = ckt->CKTrhsOld [i] ; old = ckt->CKTrhs [i] ; fprintf(stdout,"%-30s %20g %20g", node->name, new, old); if(node->type == SP_VOLTAGE) { tol = ckt->CKTreltol * (MAX(fabs(old),fabs(new))) + ckt->CKTvoltTol; } else { tol = ckt->CKTreltol * (MAX(fabs(old),fabs(new))) + ckt->CKTabstol; } if (fabs(new-old) >tol ) { fprintf(stdout," *"); } fprintf(stdout,"\n"); } i++; } */ 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 no convergence reached - NO valid Operating Point */ 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, 0.0, job->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 = HUGE_VAL; g_mif_info.breakpoint.last = HUGE_VAL; /* gtri - end - wbk - Add Breakpoint stuff */ #endif ckt->CKTstat->STATtimePts ++; /* Setting Integration Order to Backward Euler */ ckt->CKTorder = 1; /* Copying the maxStep to every deltaOld */ for(i=0;i<7;i++) { ckt->CKTdeltaOld[i]=ckt->CKTmaxStep; } /* Setting DELTA */ ckt->CKTdelta = delta; #ifdef STEPDEBUG fprintf (stderr, "delta initialized to %g\n", ckt->CKTdelta); #endif ckt->CKTsaveDelta = ckt->CKTfinalTime/50; ckt->CKTmode = (ckt->CKTmode&MODEUIC) | MODETRAN | MODEINITTRAN; /* Changing Circuit MODE */ /* modeinittran set here */ ckt->CKTag[0]=ckt->CKTag[1]=0; /* State0 copied into State1 - DEPRECATED LEGACY function - to be replaced with memmove() */ bcopy(ckt->CKTstate0, ckt->CKTstate1, (size_t) ckt->CKTnumStates * sizeof(double)); /* Statistics Initialization using a macro at the beginning of this code */ INIT_STATS(); #ifdef CLUSTER CLUsetup(ckt); #endif } else { /* saj As traninit resets CKTmode */ ckt->CKTmode = (ckt->CKTmode&MODEUIC) | MODETRAN | MODEINITPRED; /* saj */ INIT_STATS(); 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, &(job->PSSplot_td)); if(error) { fprintf(stderr, "Couldn't relink rawfile\n"); return error; } /* end saj*/ /* Skip nextTime if it isn't the firsttime! :) */ goto resume; } /* 650 */ nextTime: /* begin LTRA code addition */ if (ckt->CKTtimePoints) { ckt->CKTtimeIndex++; if (ckt->CKTtimeIndex >= ckt->CKTtimeListSize) { /* need more space */ int need; if (pss_state == STABILIZATION) need = (int) ceil((ckt->CKTstabTime - ckt->CKTtime) / maxstepsize ) ; else need = (int) ceil((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 */ /* Check for the timepoint acceptance */ error = CKTaccept(ckt); /* check if current breakpoint is outdated; if so, clear */ if (ckt->CKTtime > ckt->CKTbreaks[0]) 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 fprintf (stderr, "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) { UPDATE_STATS(DOING_TRAN); return(error); } #ifdef XSPICE /* gtri - modify - wbk - 12/19/90 - Send IPC stuff */ if(g_ipc.enabled) { if (pss_state == PSS) { /* 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, job->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 (pss_state == PSS) CLUoutput(ckt); #endif if (pss_state == PSS) { nextstep = time_temp + 1 / ckt->CKTguessedFreq * ((double)(pss_points_cycle) / (double)ckt->CKTpsspoints) ; /* If in_pss, store data for Time Domain Plot and gather ordered data for FFT computing */ if ((AlmostEqualUlps (ckt->CKTtime, nextstep, 10)) || (ckt->CKTtime > time_temp + 1 / ckt->CKTguessedFreq)) { #ifdef STEPDEBUG fprintf (stderr, "IN_PSS: time point accepted in evolution for FFT calculations.\n") ; fprintf (stderr, "Circuit time %1.15g, final time %1.15g, point index %d and total requested points %ld\n", ckt->CKTtime, nextstep, pss_points_cycle, ckt->CKTpsspoints) ; #endif CKTdump (ckt, ckt->CKTtime, job->PSSplot_td) ; /* Store the Time Base for the DFT */ psstimes [pss_points_cycle] = ckt->CKTtime ; /* Store values for the FFT calculation */ for (i = 1 ; i <= msize ; i++) pssvalues [i - 1 + pss_points_cycle * msize] = ckt->CKTrhsOld [i] ; /* Update PSS counter cycle, used to stop the entire algorithm */ pss_points_cycle++ ; /* Set the next BreakPoint for PSS */ CKTsetBreak (ckt, time_temp + (1 / ckt->CKTguessedFreq) * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ; #ifdef STEPDEBUG fprintf (stderr, "Next breakpoint set in: %1.15g\n", time_temp + 1 / ckt->CKTguessedFreq * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ; #endif } else { /* Algo can enter here but should do nothing */ #ifdef STEPDEBUG fprintf (stderr, "IN_PSS: time point accepted in evolution but dropped for FFT calculations\n") ; #endif } } #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 *******/ /* ***********************************/ switch(pss_state) { case STABILIZATION: { /* Test if stabTime has been reached */ if (AlmostEqualUlps (ckt->CKTtime, ckt->CKTstabTime, 100)) { time_temp = ckt->CKTtime ; /* Set the new Final Time - This is important because the last breakpoint is always CKTfinalTime */ ckt->CKTfinalTime = time_temp + 2 / ckt->CKTguessedFreq ; fprintf (stderr, "Exiting from stabilization\n") ; fprintf (stderr, "Time of first shooting evaluation will be %1.10g\n", time_temp + 1 / ckt->CKTguessedFreq) ; /* Next time is no more in stabilization - Unset the flag */ pss_state = SHOOTING; /* Save the RHS_copy_der as the NEW CKTrhsOld */ for (i = 1 ; i <= msize ; i++) RHS_copy_der [i - 1] = ckt->CKTrhsOld [i] ; /* Print RHS on exiting from stabilization */ fprintf (stderr, "RHS on exiting from stabilization: ") ; for (i = 1 ; i <= msize ; i++) { RHS_copy_se [i - 1] = ckt->CKTrhsOld [i] ; fprintf (stderr, "%-15g ", RHS_copy_se [i - 1]) ; } fprintf (stderr, "\n") ; /* RHS_max and RHS_min initialization - HUGE_VAL is the maximum machine error */ for (i = 0 ; i < msize ; i++) { RHS_max [i] = -HUGE_VAL ; RHS_min [i] = HUGE_VAL ; } } } break; case SHOOTING: { double offset, interval, nextBreak ; int i ; /* Calculation of error norms of RHS solution of every accepted nextTime */ err = 0 ; for (i = 0 ; i < msize ; i++) { /* Save max per node or branch of every estimated period */ if (RHS_max [i] < ckt->CKTrhsOld [i + 1]) RHS_max [i] = ckt->CKTrhsOld [i + 1] ; /* Save min per node or branch of every estimated period */ if (RHS_min [i] > ckt->CKTrhsOld [i + 1]) RHS_min [i] = ckt->CKTrhsOld [i + 1] ; /* CKTrhsOld is the last CORRECT value of RHS */ err_conv [i] = ckt->CKTrhsOld [i + 1] - RHS_copy_se [i] ; err += err_conv [i] * err_conv [i] ; /* Compute and store derivative */ RHS_derivative [i] = (ckt->CKTrhsOld [i + 1] - RHS_copy_der [i]) / ckt->CKTdelta ; /* Save the RHS_copy_der as the NEW CKTrhsOld */ RHS_copy_der [i] = ckt->CKTrhsOld [i + 1] ; #ifdef STEPDEBUG fprintf (stderr, "Pred is so high or so low! Diff is: %g\n", err_conv [i]) ; #endif } err = sqrt (err) ; /* Start frequency estimation */ if ((err < err_0) && (ckt->CKTtime >= time_temp + 0.5 / ckt->CKTguessedFreq)) /* far enough from time temp... */ { if (err < err_min_0) { err_min_1 = err_min_0 ; /* store previous minimum of RHS vector error */ err_min_0 = err ; /* store minimum of RHS vector error */ time_err_min_1 = time_err_min_0 ; /* store previous minimum of RHS vector error time */ time_err_min_0 = ckt->CKTtime ; /* store minimum of RHS vector error time */ } } 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 */ } err_1 = err ; /* *************************************** */ /* ********** Breakpoint update ********** */ /* *************************************** */ /* Force the tran analysis to evaluate requested breakpoints. Breakpoints are even more closer as the next occurence of guessed period is approaching. La lunga notte dei robot viventi... */ /* double offset, interval, nextBreak ; int i ; */ if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.995) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq))) { offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.995 ; interval = (1 / ckt->CKTguessedFreq) * (1 - 0.995) * (ckt->CKTsteady_coeff / 10) ; i = (int)((ckt->CKTtime - offset) / interval) ; nextBreak = offset + (i + 1) * interval ; CKTsetBreak (ckt, nextBreak) ; } else if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.8) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.995)) { offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.8 ; interval = (1 / ckt->CKTguessedFreq) * (0.995 - 0.8) * (ckt->CKTsteady_coeff / 5) ; i = (int)((ckt->CKTtime - offset) / interval) ; nextBreak = offset + (i + 1) * interval ; CKTsetBreak (ckt, nextBreak) ; } else if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.5) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.8)) { offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.5 ; interval = (1 / ckt->CKTguessedFreq) * (0.8 - 0.5) * (ckt->CKTsteady_coeff / 3) ; i = (int)((ckt->CKTtime - offset) / interval) ; nextBreak = offset + (i + 1) * interval ; CKTsetBreak (ckt, nextBreak) ; } else if ((ckt->CKTtime > time_temp + (1 / ckt->CKTguessedFreq) * 0.1) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.5)) { offset = time_temp + (1 / ckt->CKTguessedFreq) * 0.1 ; interval = (1 / ckt->CKTguessedFreq) * (0.5 - 0.1) * (ckt->CKTsteady_coeff / 2) ; i = (int)((ckt->CKTtime - offset) / interval) ; nextBreak = offset + (i + 1) * interval ; CKTsetBreak (ckt, nextBreak) ; } else if ((ckt->CKTtime > time_temp) && (ckt->CKTtime <= time_temp + (1 / ckt->CKTguessedFreq) * 0.1)) { offset = time_temp ; interval = (1 / ckt->CKTguessedFreq) * (0.1) * (ckt->CKTsteady_coeff) ; i = (int)((ckt->CKTtime - offset) / interval) ; nextBreak = offset + (i + 1) * interval ; CKTsetBreak (ckt, nextBreak) ; } else { fprintf (stderr, "Strange behavior\n\n") ; fprintf (stderr, "CKTtime: %g\ntime_temp: %g\n\n", ckt->CKTtime, time_temp) ; } /* *************************************** */ /* ******* END Breakpoint update ********* */ /* *************************************** */ /* If evolution is near shooting... */ if ((AlmostEqualUlps (ckt->CKTtime, time_temp + 1 / ckt->CKTguessedFreq, 10)) || (ckt->CKTtime > time_temp + 1 / ckt->CKTguessedFreq)) { int excessive_err_nodes = 0 ; /* Calculation of error norms of RHS solution of every accepted nextTime */ predsum = 0 ; for (i = 0 ; i < msize ; i++) { /* Pitagora ha sempre ragione!!! :))) */ /* pred is treated as FREQUENCY to avoid numerical overflow when derivative is close to ZERO */ pred [i] = RHS_derivative [i] / err_conv [i] ; #ifdef STEPDEBUG fprintf (stderr, "Pred is so high or so low! Diff is: %g\n", err_conv [i]) ; #endif if ((fabs (pred [i]) > 1.0e6 * ckt->CKTguessedFreq) || (err_conv [i] == 0)) { if (pred [i] > 0) pred [i] = 1.0e6 * ckt->CKTguessedFreq ; else pred [i] = -1.0e6 * ckt->CKTguessedFreq ; } predsum += pred [i] ; #ifdef STEPDEBUG fprintf (stderr, "Predsum in time before to be divided by dynamic_test has value %g\n", 1 / predsum) ; fprintf (stderr, "Current Diff: %g, Derivative: %g, Frequency Projection: %g\n", err_conv [i], RHS_derivative [i], pred [i]) ; #endif } // int excessive_err_nodes = 0 ; if (shooting_cycle_counter == 0) { /* If first time in shooting we warn about that ! */ fprintf (stderr, "In shooting...\n") ; } //#ifdef STEPDEBUG /* For debugging purpose */ fprintf (stderr, "\n----------------\n") ; fprintf (stderr, "Shooting cycle iteration number: %3d ||", shooting_cycle_counter) ; if (shooting_cycle_counter > 0) fprintf (stderr, " rr: %g || predsum: %g\n", rr_history [shooting_cycle_counter - 1], 1 / predsum) ; else fprintf (stderr, " rr: %g || predsum: %g\n", 0.0, 1 / predsum) ; // fprintf (stderr, "Print of dynamically consistent nodes voltages or branches currents:\n") ; /* --------------------- */ //#endif for (i = 0, node = ckt->CKTnodes->next ; node ; i++, node = node->next) { /* Voltage Node */ if (!strstr (node->name, "#")) { if (fabs (err_conv [i]) > (fabs (RHS_max [i] - RHS_min [i]) * ckt->CKTreltol + ckt->CKTvoltTol) * ckt->CKTtrtol * ckt->CKTsteady_coeff) { excessive_err_nodes++ ; } /* If the dynamic is below 10uV, it's dropped */ if (fabs (RHS_max [i] - RHS_min [i]) > 10 * 1e-6) { dynamic_test++ ; /* test on voltage dynamic consistence */ } /* Current Node */ } else { if (fabs (err_conv [i]) > (fabs (RHS_max [i] - RHS_min [i]) * ckt->CKTreltol + ckt->CKTabstol) * ckt->CKTtrtol * ckt->CKTsteady_coeff) { excessive_err_nodes++ ; } /* If the dynamic is below 10nA, it's dropped */ if (fabs (RHS_max [i] - RHS_min [i]) > 10 * 1e-9) { dynamic_test++ ; /* test on current dynamic consistence */ } } } if (dynamic_test == 0) { /* Test for dynamic existence */ fprintf (stderr, "No detectable dynamic on voltages nodes or currents branches. PSS analysis aborted\n") ; /* Terminates plot in Time Domain and frees the allocated memory */ SPfrontEnd->OUTendPlot (job->PSSplot_td) ; FREE (RHS_copy_se) ; FREE (RHS_copy_der) ; FREE (RHS_max) ; FREE (RHS_min) ; FREE (err_conv) ; FREE (psstimes) ; FREE (pssvalues) ; return (E_PANIC) ; /* to be corrected with definition of new error macro in iferrmsg.h */ } else if ((time_err_min_0 - time_temp) < 0) { /* Something has gone wrong... */ fprintf (stderr, "Cannot find a minimum for error vector in estimated period. Try to adjust tstab! PSS analysis aborted\n") ; /* Terminates plot in Time Domain and frees the allocated memory */ SPfrontEnd->OUTendPlot (job->PSSplot_td) ; FREE (RHS_copy_se) ; FREE (RHS_copy_der) ; FREE (RHS_max) ; FREE (RHS_min) ; FREE (err_conv) ; FREE (psstimes) ; FREE (pssvalues) ; return (E_PANIC) ; /* to be corrected with definition of new error macro in iferrmsg.h */ } //#ifdef STEPDEBUG // fprintf (stderr, "Global Convergence Error reference: %g, Time Projection: %g.\n", // err_conv_ref / dynamic_test, predsum) ; //#endif /* Take the mean value of time prediction trough the dynamic test variable - predsum becomes TIME */ predsum = 1 / (predsum * dynamic_test) ; /* Store the predsum history as absolute value */ predsum_history [shooting_cycle_counter] = fabs (predsum) ; /***********************************/ /*** FREQUENCY ESTIMATION UPDATE ***/ /***********************************/ if ((err_min_0 == err) || (err_min_0 == HUGE_VAL)) { /* Enters here if guessed frequency is higher than the 'real' value */ ckt->CKTguessedFreq = 1 / (1 / ckt->CKTguessedFreq + fabs (predsum)) ; #ifdef STEPDEBUG fprintf (stderr, "Frequency DOWN: est per %g, err min %g, err min 1 %g, err max %g, err %g\n", time_err_min_0 - time_temp, err_min_0, err_min_1, err_max, err) ; #endif } else { /* Enters here if guessed frequency is lower than the 'real' value */ ckt->CKTguessedFreq = 1 / (time_err_min_0 - time_temp) ; #ifdef STEPDEBUG fprintf (stderr, "Frequency UP: est per %g, err min %g, err min 1 %g, err max %g, err %g\n", time_err_min_0 - time_temp, err_min_0, err_min_1, err_max, err) ; #endif } /* Temporary variables to store previous occurrence of guessed frequency */ gf_last_1 = gf_last_0 ; gf_last_0 = ckt->CKTguessedFreq ; /* Next evaluation of shooting will be updated time (time_temp) summed to updated guessed period */ time_temp = ckt->CKTtime ; /* Store error history */ rr_history [shooting_cycle_counter] = err ; gf_history [shooting_cycle_counter] = ckt->CKTguessedFreq ; shooting_cycle_counter++ ; fprintf (stderr, "Updated guessed frequency: %1.10lg .\n", ckt->CKTguessedFreq) ; fprintf (stderr, "Next shooting evaluation time is %1.10g and current time is %1.10g.\n", time_temp + 1 / ckt->CKTguessedFreq, ckt->CKTtime) ; /* Restore maximum and minimum error for next search */ err_min_0 = HUGE_VAL ; err_max = -HUGE_VAL ; err_0 = HUGE_VAL ; err_1 = -HUGE_VAL ; dynamic_test = 0 ; /* Reset actual RHS reference for next shooting evaluation */ for (i = 1 ; i <= msize ; i++) RHS_copy_se [i - 1] = ckt->CKTrhsOld [i] ; #ifdef STEPDEBUG fprintf (stderr, "RHS on new shooting cycle: ") ; for (i = 0 ; i < msize ; i++) fprintf (stderr, "%-15g ", RHS_copy_se [i]) ; fprintf (stderr, "\n") ; #endif for (i = 0 ; i < msize ; i++) { /* Reset max and min per node or branch on every shooting cycle */ RHS_max [i] = -HUGE_VAL ; RHS_min [i] = HUGE_VAL ; } fprintf (stderr, "----------------\n\n") ; /* Shooting Exit Condition */ if ((shooting_cycle_counter > ckt->CKTsc_iter) || (excessive_err_nodes == 0)) { int k ; double minimum ; pss_state = PSS ; #ifdef STEPDEBUG fprintf (stderr, "\nFrequency estimation (FE) and RHS period residual (PR) evolution\n") ; #endif // minimum = rr_history [0] ; minimum = predsum_history [0] ; k = 0 ; for (i = 0 ; i < shooting_cycle_counter ; i++) { /* Print some statistics */ fprintf (stderr, "%-3d -> FE: %-15.10g || RR: %15.10g", i, gf_history [i], rr_history [i]) ; /* Take the minimum residual iteration */ // if (minimum > rr_history [i]) if (minimum > predsum_history [i]) { // minimum = rr_history [i] ; minimum = predsum_history [i] ; k = i ; } fprintf (stderr, " || predsum/dynamic_test: %15.10g || minimum: %15.10g\n", predsum_history [i], minimum) ; } if (excessive_err_nodes == 0) /* SHOOTING has converged */ ckt->CKTguessedFreq = gf_history [shooting_cycle_counter - 1] ; else ckt->CKTguessedFreq = gf_history [k] ; /* Save the current Time */ time_temp = ckt->CKTtime ; /* Set the new Final Time - This is important because the last breakpoint is always CKTfinalTime */ ckt->CKTfinalTime = time_temp + 1 / ckt->CKTguessedFreq ; /* Dump the first PSS point for the FFT */ CKTdump (ckt, ckt->CKTtime, job->PSSplot_td) ; psstimes [pss_points_cycle] = ckt->CKTtime ; for (i = 1 ; i <= msize ; i++) pssvalues [i - 1 + pss_points_cycle * msize] = ckt->CKTrhsOld [i] ; /* Update the PSS points counter and set the next Breakpoint */ pss_points_cycle++ ; CKTsetBreak (ckt, time_temp + (1 / ckt->CKTguessedFreq) * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ; if (excessive_err_nodes == 0) fprintf (stderr, "\nConvergence reached. Final circuit time is %1.10g seconds (iteration n° %d) and predicted fundamental frequency is %15.10g Hz\n", ckt->CKTtime, shooting_cycle_counter - 1, ckt->CKTguessedFreq) ; else fprintf (stderr, "\nConvergence not reached. However the most near convergence iteration has predicted (iteration %d) a fundamental frequency of %15.10g Hz\n", k, ckt->CKTguessedFreq) ; #ifdef STEPDEBUG fprintf (stderr, "time_temp %g\n", time_temp) ; fprintf (stderr, "IN_PSS: FIRST time point accepted in evolution for FFT calculations\n") ; fprintf (stderr, "Circuit time %1.15g, final time %1.15g, point index %d and total requested points %ld\n", ckt->CKTtime, time_temp + 1 / ckt->CKTguessedFreq * ((double)pss_points_cycle / (double)ckt->CKTpsspoints), pss_points_cycle, ckt->CKTpsspoints) ; fprintf (stderr, "Next breakpoint set in: %1.15g\n", time_temp + 1 / ckt->CKTguessedFreq * ((double)pss_points_cycle / (double)ckt->CKTpsspoints)) ; #endif } else { /* Set the new Final Time - This is important because the last breakpoint is always CKTfinalTime */ ckt->CKTfinalTime = time_temp + 1 / ckt->CKTguessedFreq ; /* Set next the breakpoint */ CKTsetBreak (ckt, time_temp + 1 / ckt->CKTguessedFreq) ; } } } break; case PSS: { /* The algorithm enters here when in_pss is set */ #ifdef STEPDEBUG fprintf (stderr, "ttemp %1.15g, final_time %1.15g, current_time %1.15g\n", time_temp, time_temp + 1 / ckt->CKTguessedFreq, ckt->CKTtime) ; #endif if ((pss_points_cycle == ckt->CKTpsspoints + 1) || (ckt->CKTtime > ckt->CKTfinalTime)) { double *pssfreqs = TMALLOC (double, ckt->CKTharms); double *pssmags = TMALLOC (double, ckt->CKTharms); double *pssphases = TMALLOC (double, ckt->CKTharms); double *pssnmags = TMALLOC (double, ckt->CKTharms); double *pssnphases = TMALLOC (double, ckt->CKTharms); double *pssValues = TMALLOC (double, ckt->CKTpsspoints + 1); double *pssResults = TMALLOC (double, msize * ckt->CKTharms); /* End plot in Time Domain */ SPfrontEnd->OUTendPlot (job->PSSplot_td) ; /* Frequency Plot Creation */ 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 Analysis", freqUid, IF_REAL, numNames, nameList, IF_REAL, &(job->PSSplot_fd)) ; tfree (nameList) ; SPfrontEnd->OUTattributes (job->PSSplot_fd, NULL, PLOT_COMB, NULL) ; /* ******************** */ /* Starting DFT on data */ /* ******************** */ for (i = 0 ; i < msize ; i++) { for (j = 0 ; j < ckt->CKTpsspoints ; j++) pssValues [j] = pssvalues [j * msize + i] ; DFT (ckt->CKTpsspoints, ckt->CKTharms, &thd, psstimes, pssValues, ckt->CKTguessedFreq, pssfreqs, pssmags, pssphases, pssnmags, pssnphases) ; for (j = 0 ; j < ckt->CKTharms ; j++) pssResults [j * msize + i] = pssmags [j] ; } for (j = 0 ; j < ckt->CKTharms ; j++) { for (i = 0 ; i < msize ; i++) ckt->CKTrhsOld [i + 1] = pssResults [j * msize + i] ; CKTdump (ckt, pssfreqs [j], job->PSSplot_fd) ; } /* ****************** */ /* Ending DFT on data */ /* ****************** */ /* Terminates plot in Frequency Domain and frees the allocated memory */ SPfrontEnd->OUTendPlot (job->PSSplot_fd) ; /* Francesco Lannutti's MOD */ /* Verify the frequency found */ max_freq = pssResults [msize] ; /* max_freq = pssResults [1 * msize + 0] ; */ position = 1 ; for (j = 1 ; j < ckt->CKTharms ; j++) { for (i = 0 ; i < msize ; i++) { if (max_freq < pssResults [j * msize + i]) { max_freq = pssResults [j * msize + i] ; position = j ; } } } if (pssfreqs [position] != ckt->CKTguessedFreq) { ckt->CKTguessedFreq = pssfreqs [position] ; fprintf (stderr, "The predicted fundamental frequency is incorrect.\nRelaunching the analysis...\n\n") ; fprintf (stderr, "The new guessed fundamental frequency is: %.6g\n\n", ckt->CKTguessedFreq) ; DCpss (ckt, 1) ; } /****************************/ FREE (pssResults) ; FREE (pssValues) ; FREE (pssnphases) ; FREE (pssnmags) ; FREE (pssphases) ; FREE (pssmags) ; FREE (pssfreqs) ; FREE (RHS_copy_se) ; FREE (RHS_copy_der) ; FREE (RHS_max) ; FREE (RHS_min) ; FREE (err_conv) ; FREE (psstimes) ; FREE (pssvalues) ; return (OK) ; } } break; } /* switch(pss_state) */ /* ********************************** */ /* **** END SHOOTING CODE BLOCK ***** */ /* ********************************** */ if(SPfrontEnd->IFpauseTest()) { /* user requested pause... */ UPDATE_STATS(DOING_TRAN); return(E_PAUSE); } /* RESUME */ resume: #ifdef STEPDEBUG if( (ckt->CKTdelta <= ckt->CKTfinalTime/50) && (ckt->CKTdelta <= ckt->CKTmaxStep)) { ; } else { if(ckt->CKTfinalTime/50CKTmaxStep) { fprintf (stderr, "limited by Tstop/50\n"); } else { fprintf (stderr, "limited by Tmax == %g\n", ckt->CKTmaxStep); } } #endif #ifdef HAS_PROGREP if (ckt->CKTtime == 0.) SetAnalyse( "tran init", 0); else if ((pss_state != PSS) && (shooting_cycle_counter > 0)) SetAnalyse("shooting", shooting_cycle_counter) ; 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[0]) || (ckt->CKTbreaks[0] - */ if (ckt->CKTbreaks [0] - ckt->CKTtime <= ckt->CKTdelmin) { /*if ( AlmostEqualUlps( ckt->CKTtime, ckt->CKTbreaks[0], 100 ) || (ckt->CKTbreaks[0] - * (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; #ifdef STEPDEBUG if( (ckt->CKTdelta > .1*ckt->CKTsaveDelta) || (ckt->CKTdelta > .1*(ckt->CKTbreaks[1] - ckt->CKTbreaks[0])) ) { if(ckt->CKTsaveDelta < (ckt->CKTbreaks[1] - ckt->CKTbreaks[0])) { fprintf (stderr, "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 [0], ckt->CKTtime) ; } else { fprintf (stderr, "limited by next breakpoint\n") ; fprintf (stderr, "(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] == 0) ckt->CKTdelta = ckt->CKTdelmin ; else ckt->CKTdelta = MIN (ckt->CKTdelta, .1 * MIN (ckt->CKTsaveDelta, ckt->CKTbreaks[1] - ckt->CKTbreaks[0])); if(firsttime) { ckt->CKTdelta /= 10; #ifdef STEPDEBUG fprintf(stderr, "delta cut for initial timepoint\n"); #endif } #ifndef XSPICE /* don't want to get below delmin for no reason */ ckt->CKTdelta = MAX(ckt->CKTdelta, ckt->CKTdelmin*2.0); #endif } #ifndef XSPICE else if(ckt->CKTtime + ckt->CKTdelta >= ckt->CKTbreaks[0]) { ckt->CKTsaveDelta = ckt->CKTdelta; ckt->CKTdelta = ckt->CKTbreaks[0] - ckt->CKTtime; /* fprintf (stderr, "delta cut to %g to hit breakpoint\n" ,ckt->CKTdelta) ; */ fflush(stdout); ckt->CKTbreak = 1; /* why? the current pt. is not a bkpt. */ } #endif /* !XSPICE */ #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 = HUGE_VAL; } /* gtri - end - wbk - Add Breakpoint stuff */ /* gtri - begin - wbk - Modify Breakpoint stuff */ /* Throw out any permanent breakpoint times <= current time */ for (;;) { #ifdef STEPDEBUG fprintf (stderr, " brk_pt: %g ckt_time: %g ckt_min_break: %g\n", ckt->CKTbreaks [0], ckt->CKTtime, ckt->CKTminBreak) ; #endif if(AlmostEqualUlps(ckt->CKTbreaks[0], ckt->CKTtime, 100) || ckt->CKTbreaks[0] <= ckt->CKTtime + ckt->CKTminBreak) { #ifdef STEPDEBUG fprintf (stderr, "throwing out permanent breakpoint times <= current time (brk pt: %g)\n", ckt->CKTbreaks [0]) ; fprintf (stderr, "ckt_time: %g ckt_min_break: %g\n", ckt->CKTtime, ckt->CKTminBreak) ; #endif CKTclrBreak(ckt); } else { break; } } /* Force the breakpoint if appropriate */ if(ckt->CKTtime + ckt->CKTdelta > ckt->CKTbreaks[0]) { ckt->CKTbreak = 1; ckt->CKTsaveDelta = ckt->CKTdelta; ckt->CKTdelta = ckt->CKTbreaks[0] - ckt->CKTtime; } /* gtri - end - wbk - Modify Breakpoint stuff */ /* 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 = HUGE_VAL; /* 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[0] < g_mif_info.breakpoint.current) if(ckt->CKTbreaks[0] > ckt->CKTtime + ckt->CKTminBreak) g_mif_info.breakpoint.current = ckt->CKTbreaks[0]; 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 */ #else #ifdef CLUSTER if(!CLUsync(ckt->CKTtime,&ckt->CKTdelta,0)) { fprintf (stderr, "Sync error!\n"); exit(0); } #endif /* CLUSTER */ #endif /* What is that??? */ for(i=5; i>=0; i--) ckt->CKTdeltaOld[i+1] = ckt->CKTdeltaOld[i]; ckt->CKTdeltaOld[0] = ckt->CKTdelta; temp = ckt->CKTstates[ckt->CKTmaxOrder+1]; for(i=ckt->CKTmaxOrder;i>=0;i--) { ckt->CKTstates[i+1] = ckt->CKTstates[i]; } ckt->CKTstates[0] = temp; /* 600 */ for (;;) { #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? */ 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 = HUGE_VAL; /* 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); } #ifdef STEPDEBUG if (pss_state == PSS) fprintf (stderr, "pss_state: %d, converged: %d\n", pss_state, converged) ; #endif if(converged != 0) { #ifndef CLUSTER ckt->CKTtime = ckt->CKTtime -ckt->CKTdelta; ckt->CKTstat->STATrejected ++; #endif ckt->CKTdelta = ckt->CKTdelta/8; #ifdef STEPDEBUG fprintf (stderr, "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) { firsttime = 0; #ifndef CLUSTER goto nextTime; /* no check on * first time point */ #else redostep = 0; goto chkStep; #endif } newdelta = ckt->CKTdelta; error = CKTtrunc(ckt,&newdelta); if(error) { UPDATE_STATS(DOING_TRAN); return(error); } if(newdelta > .9 * ckt->CKTdelta) { if((ckt->CKTorder == 1) && (ckt->CKTmaxOrder > 1)) { /* don't rise the order for backward Euler */ newdelta = ckt->CKTdelta; ckt->CKTorder = 2; error = CKTtrunc(ckt,&newdelta); if(error) { UPDATE_STATS(DOING_TRAN); return(error); } if(newdelta <= 1.05 * ckt->CKTdelta) { ckt->CKTorder = 1; } } /* time point OK - 630 */ ckt->CKTdelta = newdelta; #ifdef NDEV /* show a time process indicator, by Gong Ding, gdiso@ustc.edu */ if(ckt->CKTtime/ckt->CKTfinalTime*100<10.0) fprintf (stderr, "%%%3.2lf\b\b\b\b\b", ckt->CKTtime / ckt->CKTfinalTime * 100) ; else if(ckt->CKTtime/ckt->CKTfinalTime*100<100.0) fprintf (stderr, "%%%4.2lf\b\b\b\b\b\b", ckt->CKTtime / ckt->CKTfinalTime * 100) ; else fprintf (stderr, "%%%5.2lf\b\b\b\b\b\b\b", ckt->CKTtime / ckt->CKTfinalTime * 100) ; fflush(stdout); #endif #ifdef STEPDEBUG fprintf (stderr, "delta set to truncation error result: %g. Point accepted at CKTtime: %g\n", ckt->CKTdelta, ckt->CKTtime) ; fflush(stdout); #endif #ifndef CLUSTER /* go to 650 - trapezoidal */ goto nextTime; #else redostep = 0; goto chkStep; #endif } else { /* newdelta <= .9 * ckt->CKTdelta */ #ifndef CLUSTER ckt->CKTtime = ckt->CKTtime -ckt->CKTdelta; ckt->CKTstat->STATrejected ++; #endif ckt->CKTdelta = newdelta; #ifdef STEPDEBUG fprintf (stderr, "delta set to truncation error result:point rejected\n") ; #endif } } if (ckt->CKTdelta <= ckt->CKTdelmin) { if (olddelta > ckt->CKTdelmin) { ckt->CKTdelta = ckt->CKTdelmin; #ifdef STEPDEBUG fprintf (stderr, "delta at delmin\n"); #endif } else { UPDATE_STATS(DOING_TRAN); 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 */ } static int DFT ( long 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, j; double tmp; NG_IGNORE (Time); /* clear output/computation arrays */ for (i = 0; i < numFreq; i++) { Mag [i] = 0; Phase [i] = 0; } for (i = 0; i < ndata; i++) { for (j = 0; j < numFreq; j++) { Mag [j] += (Value [i] * sin (j * 2.0 * M_PI * i / ((double)ndata))); Phase [j] += (Value [i] * cos (j * 2.0 * M_PI * i / ((double)ndata))); } } Mag [0] = Phase [0] / (double)ndata; Phase [0] = 0; nMag [0] = 0; nPhase [0] = 0; Freq [0] = 0; *thd = 0; for (i = 1; i < numFreq; i++) { tmp = Mag [i] * 2.0 / (double)ndata; Phase [i] *= 2.0 / (double)ndata; Freq [i] = i * FundFreq; Mag [i] = hypot (tmp, Phase [i]); Phase [i] = atan2 (Phase [i], tmp) * 180.0 / M_PI; nMag [i] = Mag [i] / Mag [1]; nPhase [i] = Phase [i] - Phase [1]; if (i > 1) *thd += nMag [i] * nMag [i]; } *thd = 100 * sqrt (*thd); return (OK); }