Implemented the KCL verification in place of the classical SPICE3 check for the circuit currents

This version modifies every model directly to obtain fvk in a separate vector, instead of using a matrix-vector product

src/include/ngspice/cktdefs.h

@@ -110,6 +110,7 @@ struct CKTcircuit {
     double *CKTrhs;             /* current rhs value - being loaded */
     double *CKTrhsOld;          /* previous rhs value for convergence
                                    testing */
+    double *CKTfvk ;            /* KCL Verification array */
     double *CKTrhsSpare;        /* spare rhs value for reordering */
     double *CKTirhs;            /* current rhs value - being loaded
                                    (imag) */
@@ -429,6 +430,10 @@ extern void NInzIter(CKTcircuit *, int, int);
 extern int NIpred(CKTcircuit *ckt);
 #endif
 
+/* KCL Verification */
+extern int NIkclVerification (CKTcircuit *) ;
+/********************/
+
 extern IFfrontEnd *SPfrontEnd;
 
 #endif

src/maths/ni/niconv.c

@@ -45,18 +45,6 @@ NIconvTest(CKTcircuit *ckt)
 #ifdef STEPDEBUG
                 printf(" non-convergence at node (type=3) %s (fabs(new-old)>tol --> fabs(%g-%g)>%g)\n",CKTnodName(ckt,i),new,old,tol);
 		printf("    reltol: %g    voltTol: %g   (tol=reltol*(MAX(fabs(old),fabs(new))) + voltTol)\n",ckt->CKTreltol,ckt->CKTvoltTol);
-#endif /* STEPDEBUG */
-		ckt->CKTtroubleNode = i;
-		ckt->CKTtroubleElt = NULL;
-                return(1);
-            }
-        } else {
-            tol =  ckt->CKTreltol * (MAX(fabs(old),fabs(new))) +
-                    ckt->CKTabstol;
-            if (fabs(new-old) >tol ) {
-#ifdef STEPDEBUG
-                printf(" non-convergence at node (type=%d) %s (fabs(new-old)>tol --> fabs(%g-%g)>%g)\n",node->type,CKTnodName(ckt,i),new,old,tol);
-		printf("    reltol: %g    abstol: %g   (tol=reltol*(MAX(fabs(old),fabs(new))) + abstol)\n",ckt->CKTreltol,ckt->CKTabstol);
 #endif /* STEPDEBUG */
 		ckt->CKTtroubleNode = i;
 		ckt->CKTtroubleElt = NULL;
@@ -74,3 +62,42 @@ NIconvTest(CKTcircuit *ckt)
     return(0);
 #endif /* NEWCONV */
 }
+
+
+/**
+ * Routine to Verify the KCL
+ */
+
+int NIkclVerification (CKTcircuit *ckt)
+{
+    int i, size ;
+    double maximum = 0 ;
+    CKTnode *node ;
+
+    size = SMPmatSize (ckt->CKTmatrix) ;
+
+    node = ckt->CKTnodes ;
+    for (i = 1 ; i <= size ; i++)
+    {
+        if (node->type == SP_CURRENT)
+        {
+            if (maximum < fabs (ckt->CKTrhsOld [i]))
+                maximum = fabs (ckt->CKTrhsOld [i]) ;
+        }
+        node = node->next ;
+    }
+
+    node = ckt->CKTnodes ;
+    for (i = 1 ; i <= size ; i++)
+    {
+        if (node->type == SP_VOLTAGE)
+        {
+            if ((ckt->CKTfvk [i] + ckt->CKTdiagGmin * ckt->CKTrhsOld [i]) > (ckt->CKTreltol * maximum + ckt->CKTabstol))
+//            if ((ckt->CKTfvk [i]) > (ckt->CKTreltol * maximum + ckt->CKTabstol))
+                return 1 ;
+        }
+        node = node->next ;
+    }
+
+    return 0 ;
+}

src/maths/ni/niiter.c

@@ -37,6 +37,9 @@ NIiter(CKTcircuit *ckt, int maxIter)
     CKTnode *node; /* current matrix entry */
     double diff, maxdiff, damp_factor, *OldCKTstate0=NULL;
 
+    /* KCL_verified */
+    int KCL_verified ;
+
     if ( maxIter < 100 ) maxIter = 100; /* some convergence issues that get resolved by increasing max iter */
 
     iterno=0;
@@ -110,6 +113,16 @@ NIiter(CKTcircuit *ckt, int maxIter)
                 ckt->CKTniState |= NISHOULDREORDER;
             }
 
+            /* KCL Verification */
+            if ((ckt->CKTnoncon == 0) && (iterno != 1))
+            {
+                KCL_verified = NIkclVerification (ckt) ;
+                if (KCL_verified) /* NOT VERIFIED */
+                    ckt->CKTnoncon = 1 ;
+                else
+                    ckt->CKTnoncon = 0 ;
+            }
+
             if(ckt->CKTniState & NISHOULDREORDER) {
                 startTime = SPfrontEnd->IFseconds();
                 error = SMPreorder(ckt->CKTmatrix,ckt->CKTpivotAbsTol,

src/maths/ni/nireinit.c

@@ -29,6 +29,7 @@ NIreinit( CKTcircuit *ckt)
     size = SMPmatSize(ckt->CKTmatrix);
     CKALLOC(CKTrhs,size+1,double);
     CKALLOC(CKTrhsOld,size+1,double);
+    CKALLOC(CKTfvk,size+1,double);
     CKALLOC(CKTrhsSpare,size+1,double);
     CKALLOC(CKTirhs,size+1,double);
     CKALLOC(CKTirhsOld,size+1,double);

src/spicelib/analysis/cktload.c

@@ -52,6 +52,7 @@ CKTload(CKTcircuit *ckt)
     size = SMPmatSize(ckt->CKTmatrix);
     for (i = 0; i <= size; i++) {
         ckt->CKTrhs[i] = 0;
+        ckt->CKTfvk[i] = 0;
     }
     SMPclear(ckt->CKTmatrix);
 #ifdef STEPDEBUG