On the road to a better scheme for handling the deletion of

multispecies phases.

Cantera/src/equil/vcs_VolPhase.cpp

@@ -635,9 +635,9 @@ namespace VCSnonideal {
       // This is where we will start to store a better approximation 
       // for the mole fractions, when the phase doesn't exist.
       // This is currently unimplemented.
-      for (int k = 0; k < m_numSpecies; k++) {
+      //for (int k = 0; k < m_numSpecies; k++) {
-	Xmol[k] = 1.0 / m_numSpecies;
+      //	Xmol[k] = 1.0 / m_numSpecies;
-      }
+      //}
       m_existence = VCS_PHASE_EXIST_NO;
     }
     /*

Cantera/src/equil/vcs_elem.cpp

@@ -74,52 +74,54 @@ namespace VCSnonideal {
      * Require 12 digits of accuracy on non-zero constraints.
      */
     for (i = 0; i < top; ++i) {
-      if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) {
+      if (m_elementActive[i]) {
-	/*
+	if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > (fabs(m_elemAbundancesGoal[i]) * 1.0e-12)) {
-	 * This logic is for charge neutrality condition
-	 */
-	if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) {
-	  AssertThrowVCS(m_elemAbundancesGoal[i] == 0.0, "vcs_elabcheck");
-	}
-	if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) {
-	  scale = VCS_DELETE_MINORSPECIES_CUTOFF;
 	  /*
-	   * Find out if the constraint is a multisign constraint.
+	   * This logic is for charge neutrality condition
-	   * If it is, then we have to worry about roundoff error
-	   * in the addition of terms. We are limited to 13
-	   * digits of finite arithmetic accuracy.
 	   */
-	  numNonZero = 0;
+	  if (m_elType[i] == VCS_ELEM_TYPE_CHARGENEUTRALITY) {
-	  multisign = false;
+	    AssertThrowVCS(m_elemAbundancesGoal[i] == 0.0, "vcs_elabcheck");
-	  for (int kspec = 0; kspec < m_numSpeciesTot; kspec++) {
-	    eval = m_formulaMatrix[i][kspec];
-	    if (eval < 0.0) {
-	      multisign = true;
-	    }
-	    if (eval != 0.0) {
-	      scale = MAX(scale, fabs(eval * m_molNumSpecies_old[kspec]));
-	      numNonZero++;
-	    }
 	  }
-	  if (multisign) {
+	  if (m_elemAbundancesGoal[i] == 0.0 || (m_elType[i] == VCS_ELEM_TYPE_ELECTRONCHARGE)) {
-	    if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) {
+	    scale = VCS_DELETE_MINORSPECIES_CUTOFF;
-	      return FALSE;
+	    /*
+	     * Find out if the constraint is a multisign constraint.
+	     * If it is, then we have to worry about roundoff error
+	     * in the addition of terms. We are limited to 13
+	     * digits of finite arithmetic accuracy.
+	     */
+	    numNonZero = 0;
+	    multisign = false;
+	    for (int kspec = 0; kspec < m_numSpeciesTot; kspec++) {
+	      eval = m_formulaMatrix[i][kspec];
+	      if (eval < 0.0) {
+		multisign = true;
+	      }
+	      if (eval != 0.0) {
+		scale = MAX(scale, fabs(eval * m_molNumSpecies_old[kspec]));
+		numNonZero++;
+	      }
+	    }
+	    if (multisign) {
+	      if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > 1e-11 * scale) {
+		return FALSE;
+	      }
+	    } else {
+	      if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) {
+		return FALSE;
+	      }
 	    }
 	  } else {
-	    if (fabs(m_elemAbundances[i] - m_elemAbundancesGoal[i]) > VCS_DELETE_MINORSPECIES_CUTOFF) {
+	    /*
+	     * For normal element balances, we require absolute compliance
+	     * even for rediculously small numbers.
+	     */
+	    if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
+	      return FALSE;
+	    } else {
 	      return FALSE;
 	    }
 	  }
-	} else {
-	  /*
-	   * For normal element balances, we require absolute compliance
-	   * even for rediculously small numbers.
-	   */
-	  if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) {
-	    return FALSE;
-	  } else {
-	    return FALSE;
-	  }
 	}
       }
     }

Cantera/src/equil/vcs_rxnadj.cpp

@@ -303,8 +303,8 @@ namespace VCSnonideal {
 	      if (k == kspec) {
 		m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDSS;
 		if (m_SSPhase[kspec] != 1) {
-		  printf("we shouldn't be here!\n");
+		  printf("vcs_RxnStepSizes:: we shouldn't be here!\n");
-		  exit(-1);
+		  std::exit(-1);
 		}
 	      }
 #ifdef DEBUG_MODE
@@ -586,7 +586,7 @@ namespace VCSnonideal {
     double diag = hessianDiag_Ideal;
     double hessActCoef = vcs_Hessian_actCoeff_diag(irxn);
     if (hessianDiag_Ideal <= 0.0) {
-      plogf("We shouldn't be here\n");
+      plogf("vcs_Hessian_diag_adj::We shouldn't be here\n");
       exit(-1);
     }
     if (hessActCoef >= 0.0) {

Cantera/src/equil/vcs_solve.cpp

@@ -62,6 +62,7 @@ namespace VCSnonideal {
   }
 
   //   Initialize the sizes within the VCS_SOLVE object
+
   /*
    *    This resizes all of the internal arrays within the object. This routine
    *    operates in two modes. If all of the parameters are the same as it
@@ -691,7 +692,7 @@ namespace VCSnonideal {
     m_elementActive.resize(nelements, 1);
   
     /*
-     *      Copy over the element names
+     *      Copy over the element names and types
      */
     for (i = 0; i < nelements; i++) {
       m_elementName[i] = pub->ElName[i];
@@ -705,6 +706,23 @@ namespace VCSnonideal {
 	}
       }
     }
+
+    for (i = 0; i < nelements; i++) {
+      if (m_elType[i] ==  VCS_ELEM_TYPE_CHARGENEUTRALITY) {
+	if (m_elemAbundancesGoal[i] != 0.0) {
+	  if (fabs(m_elemAbundancesGoal[i]) > 1.0E-9) {
+	    plogf("Charge neutrality condition %s is signicantly nonzero, %g. Giving up\n",
+		  m_elementName[i].c_str(), m_elemAbundancesGoal[i]);
+	    std::exit(-1);
+	  } else {
+	    plogf("Charge neutrality condition %s not zero, %g. Setting it zero\n",
+		  m_elementName[i].c_str(), m_elemAbundancesGoal[i]);
+	    m_elemAbundancesGoal[i] = 0.0;
+	  }
+	  
+	}
+      }
+    }
  
     /*
      *      Copy over the species names

Cantera/src/equil/vcs_solve.h

@@ -507,7 +507,9 @@ public:
    *  Reconciles Phase existence flags with total moles in each phase.
    */
   double vcs_tmoles();
-
+#ifdef DEBUG_MODE
+  void check_tmoles() const;
+#endif
 
   //! This subroutine calculates reaction free energy changes for 
   //! all noncomponent formation reactions.
@@ -1415,7 +1417,6 @@ public:
    *
    *   length = [nspecies0][nelements0]
    */
-  //DoubleStarStar sc;
   DoubleStarStar m_stoichCoeffRxnMatrix;
 
   //! Absolute size of the stoichiometric coefficients
@@ -1718,7 +1719,7 @@ public:
 
   //! Species string name for the kth species
   /*!
-   *  SpName[k] = Species string name for the kth species
+   *   Species string name for the kth species
    */
   std::vector<std::string> m_speciesName; 
 

Cantera/src/equil/vcs_solve_TP.cpp

@@ -404,7 +404,7 @@ namespace VCSnonideal {
       if (uptodate_minors == FALSE) {
 	vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
 	vcs_dfe(VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc);
-	vcs_deltag(1, false, VCS_STATECALC_NEW);
+	vcs_deltag(1, false, VCS_STATECALC_OLD);
       }
       uptodate_minors = TRUE;
     } else {
@@ -422,7 +422,15 @@ namespace VCSnonideal {
 	plogf(" (only major species)\n");
       }
     }
-
+    /*
+     *  Calculate the total moles in each phase -> old solution
+     *   -> Needed for numerical stability when phases disappear.
+     *   -> the phase moles tend to drift off without this step.
+     */
+#ifdef DEBUG_MODE
+    check_tmoles();
+#endif
+    vcs_tmoles();
     /*
      *  Copy the old solution into the new solution as an initial guess
      */
@@ -431,6 +439,7 @@ namespace VCSnonideal {
     vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new), 
 	      VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc);
     vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_new), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc);
+    
  
     /*        Go find a new reaction adjustment -> 
      *         i.e., change in extent of reaction for each reaction. 
@@ -1278,25 +1287,32 @@ namespace VCSnonideal {
 	if (m_tPhaseMoles_old[iph] != 0.0 &&
 	    m_tPhaseMoles_old[iph]/m_totalMolNum <= VCS_DELETE_PHASE_CUTOFF) {
 	  soldel = 1;
+
 	  for (kspec = 0; kspec < m_numSpeciesRdc; kspec++) {
 	    if (m_phaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) {
 	      irxn = kspec - m_numComponents;
-	      if (kspec < m_numComponents) {
+
-		if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF) {
+	      // Both of these conditions are false and should be discarded.
-		  soldel = 0;
+	      // I think a proper special case would be if a species in a small phase
-		  break;
+	      // had a significant contribution to total element total of a single
-		}
+	      // element constraint. That's it.
-	      } else {
+
-		for (k = 0; k < m_numComponents; k++) {
+	      // if (kspec < m_numComponents) {
-		  if (m_stoichCoeffRxnMatrix[irxn][k] != 0.0) {
+	      //		if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF) {
-		    if (m_molNumSpecies_old[kspec]/m_molNumSpecies_old[k] >
+	      //soldel = 0;
-			VCS_DELETE_PHASE_CUTOFF) {
+	      //break;
-		      soldel = 0;
+	      //}
-		      break;
+	      //} else {
-		    }
+	      //for (k = 0; k < m_numComponents; k++) {
-		  }
+	      //if (m_stoichCoeffRxnMatrix[irxn][k] != 0.0) {
-		}
+	      //  if (m_molNumSpecies_old[kspec]/m_molNumSpecies_old[k] >
-	      }
+	      //VCS_DELETE_PHASE_CUTOFF) {
+	      //    soldel = 0;
+	      //    break;
+	      //  }
+	      //}
+	      //}
+	      //}
 	    }
 	  }
 	  if (soldel) {
@@ -1323,23 +1339,19 @@ namespace VCSnonideal {
       retn = vcs_basopt(FALSE, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa),
 			VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test, 
 			&usedZeroedSpecies);
-      if (retn != VCS_SUCCESS) return retn;
+      if (retn != VCS_SUCCESS) {
+#ifdef DEBUG_MODE
+	plogf("   --- BASOPT returned with an error condition\n");
+#endif
+        std::exit(-1);
+      }
       vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD);
       vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc);
+
       vcs_deltag(0, true, VCS_STATECALC_OLD);
-      uptodate_minors = TRUE;
+      iti = 0;
-      if (conv) {
+      goto  L_MAINLOOP_ALL_SPECIES ;
-	/*
+    
-	 *      HKM -> I don't understand why the code would just give
-	 *             up here in some cases.
-	 *             This should probably be taken out
-	 */	    
-	plogf(" DELETION OF MULTISPECIES PHASE. ");
-	plogf("Convergence to number of positive n(i) less than C.\n");
-	plogf("Check results to follow carefully.   \n");
-	plogendl();
-	goto L_RETURN_BLOCK;
-      }
     }
     /*************************************************************************/
     /***************** CHECK FOR ELEMENT ABUNDANCE****************************/
@@ -2418,7 +2430,7 @@ namespace VCSnonideal {
   /*
    *   When they are deleted, all of their species become active
    *   species, even though their mole numbers are set to zero.
-   *   The routine does not make the decision to eliminate multiphases.
+   *   This routine does not make the decision to eliminate multiphases.
    *
    *   Note, species in phases with zero mole numbers are still
    *   considered active. Whether the phase pops back into
@@ -2441,23 +2453,17 @@ namespace VCSnonideal {
       plogf("   --- delete_multiphase %d, %s\n", iph, Vphase->PhaseName.c_str());
     }
 #endif
-    /*
+  
-     * Zero out the total moles counters for the phase
-     */
-    m_tPhaseMoles_old[iph] = 0.0;
-    m_tPhaseMoles_new[iph] = 0.0;
-    m_deltaPhaseMoles[iph] = 0.0;
    
 
 
     /*
-     * Loop over all of the active species in the phase.
+     * Loop over all of the active noncomponent species in the phase.
      */
-    for (kspec = 0; kspec < m_numSpeciesRdc; ++kspec) {
+    for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
       if (m_phaseID[kspec] == iph) {
 	if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
 	  irxn = kspec - m_numComponents;
-	  if (irxn >= 0) {
 	    /*
 	     * calculate an extent of rxn, dx, that zeroes out the species.
 	     */
@@ -2493,40 +2499,66 @@ namespace VCSnonideal {
 	     * zeroed phase
 	     */
 	    m_rxnStatus[irxn] = VCS_SPECIES_ZEROEDMS;
-	    /*
-	     *  Changed the component mole numbers to account for the
-	     *  final extent of reaction. Make sure to keep component
-	     *  mole numbers above zero 
-	     * 
-	     */
-	    // for (j = 0; j < m_numComponents; ++j) {
-	    // m_molNumSpecies_old[j] += m_stoichCoeffRxnMatrix[irxn][j] * dx;
-	    //if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) {
-	    //if (m_molNumSpecies_old[j] < 0.0) {
-	    //  m_molNumSpecies_old[j] = 0.0;
-	    //}
-	    //}
-	    //}
-	  } 
-	  else {
-#ifdef DEBUG_MODE
-	    if (m_debug_print_lvl >= 2) {
-	      plogf("   --- delete_multiphase   One of the species is a component %d - %s with mole number %g\n",
-		    kspec, m_speciesName[kspec].c_str(),   m_molNumSpecies_old[kspec]);
-	    }
-#endif
-	    if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF * VCS_DELETE_PHASE_CUTOFF ) {
-	      plogf("   --- delete_multiphase   unknown situation error exit");
-	      plogendl();
-	      std::exit(-1);
-	    } else {
-	      m_molNumSpecies_old[kspec] = 0.0;
-	      m_molNumSpecies_new[kspec] = 0.0;
-	    }
-	  }
 	}
       }
     }
+    int jcomp;
+    double deltaLarge, dj, dxWant, dxPerm, dxPerm2;
+    for (int kcomp = 0; kcomp < m_numComponents; ++kcomp) {
+      if (m_phaseID[kcomp] == iph) {
+#ifdef DEBUG_MODE
+	if (m_debug_print_lvl >= 2) {
+	  plogf("   --- delete_multiphase   One of the species is a component %d - %s with mole number %g\n",
+		kcomp, m_speciesName[kcomp].c_str(),   m_molNumSpecies_old[kcomp]);
+	}
+#endif
+	if (m_molNumSpecies_old[kcomp] != 0.0) {
+	  deltaLarge = 0.0;
+	  for (kspec = m_numComponents; kspec < m_numSpeciesRdc; ++kspec) {
+	    irxn = kspec - m_numComponents;
+	    if (m_phaseID[kspec] != iph) {
+	      if (m_stoichCoeffRxnMatrix[irxn][kcomp] != 0.0) {
+		dxWant = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][kcomp];
+		if (dxWant + m_molNumSpecies_old[kspec]  < 0.0) {
+		  dxPerm = -m_molNumSpecies_old[kspec];
+		}
+		for (jcomp = 0;  kcomp < m_numComponents; ++kcomp) {
+		  if (jcomp != kcomp) {
+		    if (m_phaseID[jcomp] == iph) {
+		      dxPerm = 0.0;
+		    } else {
+		      dj = dxWant * m_stoichCoeffRxnMatrix[irxn][jcomp];
+		      if (dj + m_molNumSpecies_old[kcomp]  < 0.0) {
+			dxPerm2 = -m_molNumSpecies_old[kcomp] / m_stoichCoeffRxnMatrix[irxn][jcomp];
+		      }
+		      if (fabs(dxPerm2) < fabs(dxPerm)) {
+			dxPerm = dxPerm2;
+		      }
+		    }
+		  }
+		}
+	      }
+	      if (dxPerm != 0.0) {
+		delta_species(kspec, &dxPerm);
+	      }
+	    }
+	  }
+
+	}
+	if (m_molNumSpecies_old[kcomp] != 0.0) {
+#ifdef DEBUG_MODE
+	  if (m_debug_print_lvl >= 2) {
+	    plogf("   --- delete_multiphase   One of the species is a component %d - %s still with mole number %g\n",
+		  kcomp, m_speciesName[kcomp].c_str(),   m_molNumSpecies_old[kcomp]);
+	    plogf("   ---                     zeroing it \n");
+	  }
+#endif
+	  m_molNumSpecies_old[kcomp] = 0.0;
+	}
+      }
+    }
+  
+
     /*
      *   Loop over all of the inactive species in the phase:
      *   Right now we reinstate all species in a deleted multiphase.
@@ -2561,12 +2593,18 @@ namespace VCSnonideal {
 	}
       }
     }
+
+    /*
+     * Zero out the total moles counters for the phase
+     */
+    m_tPhaseMoles_old[iph] = 0.0;
+    m_tPhaseMoles_new[iph] = 0.0;
+    m_deltaPhaseMoles[iph] = 0.0;
     /*
      * Upload the state to the VP object
      */
-    Vphase->setMolesFromVCSCheck(VCS_STATECALC_OLD, 
+    Vphase->setTotalMoles(0.0);
-				 VCS_DATA_PTR(m_molNumSpecies_old),
+    
-				 VCS_DATA_PTR(m_tPhaseMoles_old));
     return successful;
   } 
   /**********************************************************************************/
@@ -3275,8 +3313,8 @@ namespace VCSnonideal {
 	  ncTrial = m_numComponents;
 	  int numPreDeleted = m_numRxnTot - m_numRxnRdc;
 	  if (numPreDeleted != (m_numSpeciesTot - m_numSpeciesRdc)) {
-	    plogf("we shouldn't be here\n");
+	    plogf("vcs_basopt:: we shouldn't be here\n");
-	    exit(-1);
+	    std::exit(-1);
 	  }
 	  m_numRxnTot = m_numSpeciesTot - ncTrial;
 	  m_numRxnRdc = m_numRxnTot - numPreDeleted;
@@ -4434,6 +4472,33 @@ namespace VCSnonideal {
     m_totalMolNum = sum;
     return m_totalMolNum;
   }
+
+#ifdef DEBUG_MODE
+  void VCS_SOLVE::check_tmoles() const {
+    int i;
+    double sum, m_tPhaseMoles_old_a;
+    //vcs_VolPhase *Vphase;
+    for (i = 0; i < m_numPhases; i++) {
+      m_tPhaseMoles_old_a = TPhInertMoles[i];
+    
+      for (int k = 0; k < m_numSpeciesTot; k++) {
+	if (m_speciesUnknownType[k] == VCS_SPECIES_TYPE_MOLNUM) {
+	  if (m_phaseID[k] == i) {
+	    m_tPhaseMoles_old_a += m_molNumSpecies_old[k];
+	  }
+	}
+      }
+      sum += m_tPhaseMoles_old_a;
+ 
+      double denom = m_tPhaseMoles_old[i]+ m_tPhaseMoles_old_a + 1.0E-19;
+      if (!vcs_doubleEqual(m_tPhaseMoles_old[i]/denom, m_tPhaseMoles_old_a/denom)) {
+	plogf("check_tmoles: we have found a problem with phase  %d: %20.15g, %20.15g\n",
+	      i, m_tPhaseMoles_old[i], m_tPhaseMoles_old_a);
+	//std::exit(-1);
+      }
+    }
+  }
+#endif
   /*****************************************************************************/
 
   //  This routine uploads the state of the system into all of the 
@@ -5052,8 +5117,8 @@ namespace VCSnonideal {
     // If it isn't, we don't know what's happening
     if (m_molNumSpecies_old[kspec] != 0.0) {
       w_kspec = 0.0;
-      plogf("we shouldn't be here\n");
+      plogf("vcs_birthGuess:: we shouldn't be here\n");
-      exit(-1);
+      std::exit(-1);
     }
 #endif
     int ss = m_SSPhase[kspec];