From eb35f9e2532769e1f356f5400bce5cbad1256ed4 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Fri, 31 Dec 2010 23:20:55 +0000
Subject: [PATCH] Some changes to reduce extra printing Some changes due to a
 bad algorithm - Known - need to fix this.

---
 Cantera/src/equil/vcs_MultiPhaseEquil.cpp | 19 +++++++++----------
 Cantera/src/equil/vcs_phaseStability.cpp  | 12 ++++++++----
 2 files changed, 17 insertions(+), 14 deletions(-)

diff --git a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp
index aa7eac6e6..5446bb00c 100644
--- a/Cantera/src/equil/vcs_MultiPhaseEquil.cpp
+++ b/Cantera/src/equil/vcs_MultiPhaseEquil.cpp
@@ -1547,8 +1547,7 @@ namespace VCSnonideal {
     if (!m_vsolvePtr) {
       m_vsolvePtr = new VCS_SOLVE();
     }
-    double feStable;
-    int iStable = m_vsolvePtr->vcs_PS(m_vprob, iph, printLvl, feStable);
+    int iStable = m_vsolvePtr->vcs_PS(m_vprob, iph, printLvl, funcStab);
 
     /*
      * Transfer the information back to the MultiPhase object.
@@ -1559,13 +1558,13 @@ namespace VCSnonideal {
      * states.
      */
     m_mix->uploadMoleFractionsFromPhases();
-   for (int i = 0; i < m_vprob->nspecies; i++) {
-     plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
-   }
+    //  for (int i = 0; i < m_vprob->nspecies; i++) {
+    // plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
+    //}
     m_mix->getChemPotentials(DATA_PTR(m_vprob->m_gibbsSpecies));
-   for (int i = 0; i < m_vprob->nspecies; i++) {
-     plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
-   }
+    //for (int i = 0; i < m_vprob->nspecies; i++) {
+    // plogf("%d %15.3e\n", m_vprob->m_gibbsSpecies[i]);
+    //}
 
     double te = tickTock.secondsWC();
     if (printLvl > 0) {
@@ -1576,9 +1575,9 @@ namespace VCSnonideal {
       plogf("Pressure    = %g Pa\n", m_vprob->PresPA);
       std::string sss = m_mix->phaseName(iph);
       if (iStable) {
-	plogf("Phase %d named %s is     stable, function value = %g > 0\n", iph, sss.c_str(), feStable);
+	plogf("Phase %d named %s is     stable, function value = %g > 0\n", iph, sss.c_str(), funcStab);
       } else {
-	plogf("Phase %d named %s is not stable + function value = %g < 0\n", iph, sss.c_str(), feStable);
+	plogf("Phase %d named %s is not stable + function value = %g < 0\n", iph, sss.c_str(), funcStab);
       }
       plogf("\n");
       plogf("----------------------------------------"
diff --git a/Cantera/src/equil/vcs_phaseStability.cpp b/Cantera/src/equil/vcs_phaseStability.cpp
index 6d323da8c..3307f6954 100644
--- a/Cantera/src/equil/vcs_phaseStability.cpp
+++ b/Cantera/src/equil/vcs_phaseStability.cpp
@@ -97,7 +97,7 @@ namespace VCSnonideal {
 	 * We loop through the regular reaction looking for a reaction that can pop the
 	 * component.
 	 */
-        printf("WE are here at new logic - CHECK\n");
+        //printf("WE are here at new logic - CHECK\n");
 	for (int jrxn = 0; jrxn < m_numRxnRdc; jrxn++) {
 	  bool foundJrxn = false;
 	  // First, if the component is a product of the reaction
@@ -111,7 +111,7 @@ namespace VCSnonideal {
 	      }
 	    }
             if (foundJrxn) {
-	      printf("We have found a component phase pop! CHECK1 \n");
+	      //printf("We have found a component phase pop! CHECK1 \n");
 	      return true;
 	    }
 	  }
@@ -131,7 +131,7 @@ namespace VCSnonideal {
 	      }
 	    }
             if (foundJrxn) {
-	      printf("We have found a component phase pop! CHECK2 \n");
+	      //printf("We have found a component phase pop! CHECK2 \n");
 	      return true;
 	    }
 	  }
@@ -717,7 +717,7 @@ namespace VCSnonideal {
 	      "  normUpdate     damp     FuncPhaseStability\n", KP, KP, KP, KP);
 	plogf("   --------------------------------------------------------------"
 	      "--------------------------------------------------------\n");
-      } else {
+      } else if (m_debug_print_lvl == 1) {
 	plogf("   --- vcs_phaseStabilityTest() called for phase %d\n", iph);
       }
 #endif
@@ -757,6 +757,10 @@ namespace VCSnonideal {
 	for (k = 0; k < Vphase->nSpecies(); k++) {
 	  sumFrac += fracDelta_old[k];
 	}
+        // Necessary because this can be identically zero. -> we need to fix this algorithm!
+	if (sumFrac <= 0.0) {
+	  sumFrac = 1.0;
+	}
 	double sum_Xcomp = 0.0;
 	for (k = 0; k < Vphase->nSpecies(); k++) {
 	  X_est[k] = fracDelta_old[k] / sumFrac;