Enable optional voltage based truncation error correction.

This is selectable as 'option newtrunc' (--enable-Predictor is required)

Remove --enable-NEWTRUNC

Add three option parameters (default):
lteTrtol (500), lteReltol (1e-3), lteAbstol (1e-6)

Add new error function for TRAP, similar to GEAR
 (tentative, to be improved), in ckttrun.c

Enable PREDICTOR as default with Visual Studio.

configure.ac

@@ -229,10 +229,6 @@ AC_ARG_ENABLE([predictor],
 AC_ARG_ENABLE([newpred],
     [AS_HELP_STRING([--enable-newpred], [Enable NEWPRED whatever it is(?)])])
 
-# --enable-newtrunc: define NEWTRUNC for the code
-AC_ARG_ENABLE([newtrunc],
-    [AS_HELP_STRING([--enable-newtrunc], [Enable, how we want extrapolate capacitances.])])
-
 # --enable-sense2: define WANT_SENSE2 for the code
 AC_ARG_ENABLE([sense2],
     [AS_HELP_STRING([--enable-sense2], [Use spice2 sensitivity analysis.])])
@@ -1044,10 +1040,6 @@ if test "x$enable_newpred" = xyes; then
     AC_DEFINE([NEWPRED], [], [Define if you want to discover :)])
     AC_MSG_RESULT([NEWPRED enabled])
 fi
-if test "x$enable_newtrunc" = xyes; then
-    AC_DEFINE([NEWTRUNC], [], [Do not trigger unwanted traps by default])
-    AC_MSG_RESULT([New truncation error calculation enabled])
-fi
 if test "x$enable_experimental" = xyes; then
     AC_DEFINE([EXPERIMENTAL_CODE], [], [Define if we want some experimental code])
     AC_MSG_RESULT([EXPERIMENTAL_CODE enabled])

src/frontend/com_option.c

@@ -77,13 +77,16 @@ com_option(wordlist *wl)
         printf("gminsteps = %d\n", circuit->CKTnumGminSteps);
         printf("srcsteps = %d\n", circuit->CKTnumSrcSteps);
 
-        printf("\nTruncation error correction:\n");
+        if (circuit->CKTnewtrunc) {
-        printf("trtol = %f\n", circuit->CKTtrtol);
+            printf("\nTruncation error correction, voltage based, is selected:\n");
-#ifdef NEWTRUNC
             printf("ltereltol = %g\n", circuit->CKTlteReltol);
             printf("lteabstol = %g\n", circuit->CKTlteAbstol);
-#endif /* NEWTRUNC */
+            printf("ltetrtol = %g\n", circuit->CKTlteTrtol);
-
+        }
+        else {
+            printf("\nTruncation error correction, charge based, is selected:\n");
+            printf("trtol = %f\n", circuit->CKTtrtol);
+        }
         printf("\nConductances:\n");
         printf("gmin     (devices)  = %g\n", circuit->CKTgmin);
         printf("diaggmin (stepping) = %g\n", circuit->CKTdiagGmin);

src/frontend/misccoms.c

@@ -274,9 +274,6 @@ com_version(wordlist *wl)
 #ifdef PREDICTOR
             fprintf(cp_out, "** --enable-predictor\n");
 #endif
-#ifdef NEWTRUNC
-            fprintf(cp_out, "** --enable-newtrunc\n");
-#endif
 #ifdef WANT_SENSE2
             fprintf(cp_out, "** --enable-sense2\n");
 #endif

src/include/ngspice/cktdefs.h

@@ -216,11 +216,10 @@ struct CKTcircuit {
     double CKTreltol;           /* --- */
     double CKTchgtol;           /* --- */
     double CKTvoltTol;          /* --- */
-    /* What is this define for  ? */
+    double CKTlteReltol;        /* relative error in voltage based truncation error estimation */
-#ifdef NEWTRUNC
+    double CKTlteAbstol;        /* absolute error in voltage based truncation error estimation */
-    double CKTlteReltol;
+    double CKTlteTrtol;         /* scaling time step in voltage based truncation error estimation */
-    double CKTlteAbstol;
+    int CKTnewtrunc;            /* enable lte (local truncation error) based on voltages */
-#endif /* NEWTRUNC */
     double CKTgmin;             /* .options GMIN */
     double CKTgshunt;           /* .options RSHUNT */
     double CKTcshunt;           /* .options CSHUNT */

src/include/ngspice/optdefs.h

@@ -132,6 +132,10 @@ enum {
     OPT_KLU_MEMGROW_FACTOR,
 #endif
 
+    OPT_LTERELTOL,
+    OPT_LTEABSTOL,
+    OPT_LTETRTOL,
+    OPT_NEWTRUNC,
 };
 
 #ifdef XSPICE

src/include/ngspice/tskdefs.h

@@ -48,10 +48,10 @@ struct TSKtask {
     double TSKreltol;
     double TSKchgtol;
     double TSKvoltTol;
-#ifdef NEWTRUNC
     double TSKlteReltol;
     double TSKlteAbstol;
-#endif /* NEWTRUNC */
+    double TSKlteTrtol;
+    unsigned int TSKnewtrunc:1; /* voltage controlled truncation */
     double TSKgmin;
     double TSKgshunt;   /* shunt conductance (CKTdiagGmin) */
     double TSKcshunt;   /* shunt capacitor to ground */

src/spicelib/analysis/cktdojob.c

@@ -114,14 +114,17 @@ CKTdoJob(CKTcircuit* ckt, int reset, TSKtask* task)
     ckt->CKTkluMemGrowFactor = task->TSKkluMemGrowFactor ;
 #endif
 
-#ifdef NEWTRUNC
     ckt->CKTlteReltol = task->TSKlteReltol;
     ckt->CKTlteAbstol = task->TSKlteAbstol;
-#endif /* NEWTRUNC */
+    ckt->CKTlteTrtol = task->TSKlteTrtol;
+    ckt->CKTnewtrunc = task->TSKnewtrunc;
 
     fprintf(stdout, "Doing analysis at TEMP = %f and TNOM = %f\n\n",
         ckt->CKTtemp - CONSTCtoK, ckt->CKTnomTemp - CONSTCtoK);
 
+    if (ckt->CKTnewtrunc)
+        fprintf(stdout, "Note: Voltage based truncation error correction selected\n");
+
     /* call altermod and alter on device and model parameters assembled in
        devtlist and modtlist (if using temper) because we have a new temperature */
     inp_evaluate_temper(ft_curckt);

src/spicelib/analysis/cktntask.c

@@ -82,11 +82,10 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr)
         tsk->TSKkluMemGrowFactor = def->TSKkluMemGrowFactor ;
 #endif
 
-#ifdef NEWTRUNC
         tsk->TSKlteReltol       = def->TSKlteReltol;
         tsk->TSKlteAbstol       = def->TSKlteAbstol;
-#endif
+        tsk->TSKlteTrtol       = def->TSKlteTrtol;
-
+        tsk->TSKnewtrunc       = def->TSKnewtrunc;
     } else {
 #endif /*CDHW*/
 
@@ -98,11 +97,11 @@ CKTnewTask(CKTcircuit *ckt, TSKtask **taskPtr, IFuid taskName, TSKtask **defPtr)
         tsk->TSKreltol          = 1e-3;
         tsk->TSKchgtol          = 1e-14;
         tsk->TSKvoltTol         = 1e-6;
-#ifdef NEWTRUNC
         tsk->TSKlteReltol       = 1e-3;
         tsk->TSKlteAbstol       = 1e-6;
-#endif
+        tsk->TSKlteTrtol        = 500.;
-        tsk->TSKtrtol           = 7;
+        tsk->TSKnewtrunc        = 0;
+        tsk->TSKtrtol           = 7.;
         tsk->TSKbypass          = 0;
         tsk->TSKtranMaxIter     = 10;
         tsk->TSKdcMaxIter       = 100;

src/spicelib/analysis/cktsopt.c

@@ -187,7 +187,24 @@ CKTsetOpt(CKTcircuit *ckt, JOB *anal, int opt, IFvalue *val)
         task->TSKkluMemGrowFactor = (val->rValue == 1.2);
         break;
 #endif
-
+    case OPT_LTERELTOL:
+        task->TSKlteReltol = val->rValue;
+        break;
+    case OPT_LTEABSTOL:
+        task->TSKlteAbstol = val->rValue;
+        break;
+    case OPT_LTETRTOL:
+        task->TSKlteTrtol = val->rValue;
+        break;
+    case OPT_NEWTRUNC:
+#ifdef PREDICTOR
+        task->TSKnewtrunc = (val->iValue != 0);
+#else
+        task->TSKnewtrunc = 0;
+        fprintf(stderr, "Warning: Option 'newtrunc' ignored,\n"
+            "    compilation with preprocessor flag 'PREDICTOR' is required.\n");
+#endif
+        break;
 /* gtri - begin - wbk - add new options */
 #ifdef XSPICE
     case OPT_EVT_MAX_OP_ALTER:
@@ -230,7 +247,7 @@ CKTsetOpt(CKTcircuit *ckt, JOB *anal, int opt, IFvalue *val)
           ckt->enh->rshunt_data.gshunt = 1.0 / val->rValue;
         }
         else {
-          printf("WARNING - Rshunt option too small.  Ignored.\n");
+          fprintf(stderr, "WARNING - Rshunt option too small.  Ignored.\n");
         }
         break;
 #endif
@@ -353,9 +370,14 @@ static IFparm OPTtbl[] = {
  { "klu", OPT_KLU, IF_SET|IF_FLAG,
         "Set KLU as Direct Linear Solver" },
  { "klu_memgrow_factor", OPT_KLU_MEMGROW_FACTOR, IF_SET|IF_REAL,
-        "KLU Memory Grow Factor (default is 1.2)" }
+        "KLU Memory Grow Factor (default is 1.2)" },
 #endif
 
+ { "ltereltol", OPT_LTERELTOL,IF_SET | IF_REAL ,"Relative error tolerence" },
+ { "lteabstol", OPT_LTEABSTOL,IF_SET | IF_REAL,"Absolute error tolerence" },
+ { "ltetrtol", OPT_LTETRTOL,IF_SET | IF_REAL,"Truncation error overestimation factor" },
+ { "newtrunc", OPT_NEWTRUNC,IF_SET | IF_FLAG,"voltage controlled truncation" }
+
 };
 
 int OPTcount = NUMELEMS(OPTtbl);

src/spicelib/analysis/ckttrunc.c

@@ -15,11 +15,11 @@ Author: 1985 Thomas L. Quarles
 #include "ngspice/devdefs.h"
 #include "ngspice/sperror.h"
 
-
+//#define STEPDEBUG
 int
 CKTtrunc(CKTcircuit *ckt, double *timeStep)
 {
-#ifndef NEWTRUNC
+    if (ckt->CKTnewtrunc == 0) {
         int i;
         double timetemp;
 #ifdef STEPDEBUG
@@ -54,7 +54,8 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
 
         ckt->CKTstat->STATtranTruncTime += SPfrontEnd->IFseconds() - startTime;
         return(OK);
-#else /* NEWTRUNC */
+    }
+    else {
         int i;
         CKTnode* node;
         double timetemp;
@@ -69,12 +70,15 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
         timetemp = HUGE;
         size = SMPmatSize(ckt->CKTmatrix);
 #ifdef STEPDEBUG
-    printf("at time %g, delta %g\n",ckt->CKTtime,ckt->CKTdeltaOld[0]);
+        printf("\nNEWTRUNC at time %g, delta %g\n", ckt->CKTtime, ckt->CKTdeltaOld[0]);
 #endif
         node = ckt->CKTnodes;
         switch (ckt->CKTintegrateMethod) {
 
         case TRAPEZOIDAL:
+#ifdef STEPDEBUG
+            printf("TRAP, Order is %d\n", ckt->CKTorder);
+#endif
             switch (ckt->CKTorder) {
             case 1:
                 for (i = 1; i < size; i++) {
@@ -88,13 +92,15 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
                         ckt->CKTrhs[i], ckt->CKTpred[i]);
 #endif
                     if (diff != 0) {
-                    tmp = ckt->CKTtrtol * tol * 2 /diff;
+                        //                if (!AlmostEqualUlps(diff, 0, 10)) {
-                    tmp = ckt->CKTdeltaOld[0]*sqrt(fabs(tmp));
+                        tmp = ckt->CKTlteTrtol * tol * 2 / diff;
+                        tmp = ckt->CKTdeltaOld[0] * fabs(tmp);
                         timetemp = MIN(timetemp, tmp);
 #ifdef STEPDEBUG
                         printf("tol = %g, diff = %g, h->%g\n", tol, diff, tmp);
 #endif
-                } else {
+                    }
+                    else {
 #ifdef STEPDEBUG
                         printf("diff is 0\n");
 #endif
@@ -108,14 +114,44 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
                     node = node->next;
                     if (node->type != SP_VOLTAGE) continue;
                     diff = ckt->CKTrhs[i] - ckt->CKTpred[i];
+#ifdef STEPDEBUG
+                    printf("%s: cor=%g, pred=%g ", node->name,
+                        ckt->CKTrhs[i], ckt->CKTpred[i]);
+#endif
+                    if (diff != 0) {
+                        //                if (!AlmostEqualUlps(diff, 0, 10)) {
+                        tmp = ckt->CKTlteTrtol * tol * 2 / diff;
+                        tmp = ckt->CKTdeltaOld[0] * sqrt(fabs(tmp));
+                        timetemp = MIN(timetemp, tmp);
+#ifdef STEPDEBUG
+                        printf("tol = %g, diff = %g, h->%g\n", tol, diff, tmp);
+#endif
+                    }
+                    else {
+#ifdef STEPDEBUG
+                        printf("diff is 0\n");
+#endif
+                    }
+                }
+                break;
+                /*        case 2:
+                            for(i=1;i<size;i++) {
+                                tol = MAX( fabs(ckt->CKTrhs[i]),fabs(ckt->CKTpred[i]))*
+                                        ckt->CKTlteReltol+ckt->CKTlteAbstol;
+                                node = node->next;
+                                if(node->type!= SP_VOLTAGE) continue;
+                                diff = ckt->CKTrhs[i]-ckt->CKTpred[i];
                 #ifdef STEPDEBUG
                                 printf("%s: cor=%g, pred=%g ",node->name,ckt->CKTrhs[i],
                                         ckt->CKTpred[i]);
                 #endif
                                 if(diff != 0) {
-                    tmp = ckt->CKTdeltaOld[0]*ckt->CKTtrtol * tol * 3 * 
+                //                if(!AlmostEqualUlps(diff, 0, 10)) {
-                            (ckt->CKTdeltaOld[0]+ckt->CKTdeltaOld[1])/diff;
+                                    tmp = tol * ckt->CKTlteTrtol * (ckt->CKTdeltaOld[0] + ckt->CKTdeltaOld[1])
+                                        / (diff * ckt->CKTdelta);
                                     tmp = fabs(tmp);
+                                    tmp = exp(log(tmp) / 3);
+                                    tmp *= ckt->CKTdelta;
                                     timetemp = MIN(timetemp,tmp);
                 #ifdef STEPDEBUG
                                     printf("tol = %g, diff = %g, h->%g\n",tol,diff,tmp);
@@ -126,7 +162,7 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
                 #endif
                                 }
                             }
-            break;
+                            break; */
             default:
                 return(E_ORDER);
                 break;
@@ -135,6 +171,9 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
             break;
 
         case GEAR: {
+#ifdef STEPDEBUG
+            printf("GEAR, Order is %d\n", ckt->CKTorder);
+#endif
             double delsum = 0;
             for (i = 0; i <= ckt->CKTorder; i++) {
                 delsum += ckt->CKTdeltaOld[i];
@@ -150,7 +189,7 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
                     ckt->CKTpred[i]);
 #endif
                 if (diff != 0) {
-                tmp = tol*ckt->CKTtrtol*delsum/(diff*ckt->CKTdelta);
+                    tmp = tol * ckt->CKTlteTrtol * delsum / (diff * ckt->CKTdelta);
                     tmp = fabs(tmp);
                     switch (ckt->CKTorder) {
                     case 0:
@@ -167,7 +206,8 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
 #ifdef STEPDEBUG
                     printf("tol = %g, diff = %g, h->%g\n", tol, diff, tmp);
 #endif
-            } else {
+                }
+                else {
 #ifdef STEPDEBUG
                     printf("diff is 0\n");
 #endif
@@ -183,5 +223,5 @@ CKTtrunc(CKTcircuit *ckt, double *timeStep)
         *timeStep = MIN(2 * *timeStep, timetemp);
         ckt->CKTstat->STATtranTruncTime += SPfrontEnd->IFseconds() - startTime;
         return(OK);
-#endif /* NEWTRUNC */
+    }
 }

src/spicelib/devices/devsup.c

@@ -814,13 +814,16 @@ DEVcap(CKTcircuit *ckt, double vgd, double vgs, double vgb, double covlgd,
 double
 DEVpred(CKTcircuit *ckt, int loct)
 {
-#ifndef NEWTRUNC
+    if (!ckt->CKTnewtrunc) {
         double xfact;
 
         xfact = ckt->CKTdelta / ckt->CKTdeltaOld[1];
         return(((1 + xfact) * *(ckt->CKTstate1 + loct)) -
             (xfact * *(ckt->CKTstate2 + loct)));
-#endif /* NEWTRUNC */
+    }
+    else {
+        return 1.;
+    }
 }
 
 

src/xspice/icm/analog/delay/cfunc.mod

@@ -142,11 +142,10 @@ struct CKTcircuitmin {
     double CKTreltol;           /* --- */
     double CKTchgtol;           /* --- */
     double CKTvoltTol;          /* --- */
-    /* What is this define for  ? */
+    double CKTlteReltol;        /* relative error in voltage based truncation error estimation */
-#ifdef NEWTRUNC
+    double CKTlteAbstol;        /* absolute error in voltage based truncation error estimation */
-    double CKTlteReltol;
+    double CKTlteTrtol;         /* scaling time step in voltage based truncation error estimation */
-    double CKTlteAbstol;
+    int CKTnewtrunc;            /* enable lte (local truncation error) based on voltages */
-#endif /* NEWTRUNC */
     double CKTgmin;             /* .options GMIN */
     double CKTgshunt;           /* .options RSHUNT */
     double CKTcshunt;           /* .options CSHUNT */

visualc/src/include/ngspice/config.h

@@ -464,9 +464,6 @@
 /* Define if you want to discover :) */
 /* #undef NEWPRED */
 
-/* Do not trigger unwanted traps by default */
-/* #undef NEWTRUNC */
-
 /* Define if we want NOBYPASS */
 /* #undef NOBYPASS */
 
@@ -490,7 +487,7 @@
 #define PACKAGE_VERSION VERSION
 
 /* Define if we want predictor algorithm */
-/* #undef PREDICTOR */
+#define PREDICTOR
 
 /* Define to 1 if the C compiler supports function prototypes. */
 /* #undef PROTOTYPES */