From 789574ac8030a3686b352991f373f3101c057745 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Tue, 27 May 2014 02:52:21 +0000 Subject: [PATCH] Remove commented-out code --- src/apps/cti2ctml.cpp | 3 - src/base/application.cpp | 1 - src/base/global.cpp | 1 - src/base/stringUtils.cpp | 1 - src/base/xml.cpp | 10 ++- src/clib/ct.cpp | 2 - src/equil/ChemEquil.cpp | 15 ----- src/equil/MultiPhase.cpp | 14 +--- src/equil/MultiPhaseEquil.cpp | 28 +------- src/equil/vcs_MultiPhaseEquil.cpp | 8 --- src/equil/vcs_TP.cpp | 21 ------ src/equil/vcs_VolPhase.cpp | 19 ------ src/equil/vcs_equilibrate.cpp | 1 - src/equil/vcs_inest.cpp | 2 - src/equil/vcs_phaseStability.cpp | 3 - src/equil/vcs_prob.cpp | 3 - src/equil/vcs_rxnadj.cpp | 6 -- src/equil/vcs_setMolesLinProg.cpp | 4 -- src/equil/vcs_solve.cpp | 1 - src/equil/vcs_solve_TP.cpp | 32 --------- src/equil/vcs_solve_phaseStability.cpp | 18 ----- src/fortran/fctxml.cpp | 1 - src/kinetics/AqueousKinetics.cpp | 1 - src/kinetics/InterfaceKinetics.cpp | 40 ------------ src/kinetics/ReactionPath.cpp | 49 ++------------ src/kinetics/importKinetics.cpp | 14 +--- src/kinetics/solveSP.cpp | 6 -- src/matlab/ctmethods.cpp | 2 - src/matlab/flowdevicemethods.cpp | 2 - src/matlab/reactornetmethods.cpp | 2 - src/matlab/wallmethods.cpp | 2 - src/numerics/BEulerInt.cpp | 35 ---------- src/numerics/BandMatrix.cpp | 1 - src/numerics/CVodeInt.cpp | 2 - src/numerics/CVodesIntegrator.cpp | 13 ---- src/numerics/DenseMatrix.cpp | 2 +- src/numerics/Func1.cpp | 8 --- src/numerics/IDA_Solver.cpp | 1 - src/numerics/NonlinearSolver.cpp | 30 +-------- src/numerics/SquareMatrix.cpp | 3 - src/numerics/solveProb.cpp | 13 ---- src/oneD/Domain1D.cpp | 1 - src/oneD/MultiNewton.cpp | 2 +- src/oneD/OneDim.cpp | 1 - src/oneD/Sim1D.cpp | 11 +--- src/oneD/StFlow.cpp | 2 - src/oneD/boundaries1D.cpp | 18 +---- src/thermo/Elements.cpp | 2 - src/thermo/GeneralSpeciesThermo.cpp | 4 -- src/thermo/GibbsExcessVPSSTP.cpp | 3 - src/thermo/HMWSoln.cpp | 85 ------------------------ src/thermo/HMWSoln_input.cpp | 20 ------ src/thermo/IdealGasPhase.cpp | 2 - src/thermo/IdealMolalSoln.cpp | 1 - src/thermo/IonsFromNeutralVPSSTP.cpp | 32 --------- src/thermo/LatticePhase.cpp | 11 ---- src/thermo/LatticeSolidPhase.cpp | 4 -- src/thermo/MargulesVPSSTP.cpp | 40 +----------- src/thermo/MineralEQ3.cpp | 6 -- src/thermo/MixedSolventElectrolyte.cpp | 36 +--------- src/thermo/MixtureFugacityTP.cpp | 32 +-------- src/thermo/Nasa9Poly1.cpp | 2 - src/thermo/PDSS_ConstVol.cpp | 1 - src/thermo/PDSS_HKFT.cpp | 9 --- src/thermo/PDSS_IdealGas.cpp | 2 - src/thermo/PDSS_Water.cpp | 6 -- src/thermo/Phase.cpp | 4 -- src/thermo/PhaseCombo_Interaction.cpp | 7 +- src/thermo/PureFluidPhase.cpp | 1 - src/thermo/RedlichKisterVPSSTP.cpp | 9 +-- src/thermo/RedlichKwongMFTP.cpp | 21 ------ src/thermo/SemiconductorPhase.cpp | 14 ---- src/thermo/SingleSpeciesTP.cpp | 1 - src/thermo/SpeciesThermoFactory.cpp | 3 - src/thermo/SpeciesThermoInterpType.cpp | 1 - src/thermo/StatMech.cpp | 2 - src/thermo/VPSSMgr.cpp | 32 --------- src/thermo/VPSSMgr_ConstVol.cpp | 1 - src/thermo/VPSSMgr_General.cpp | 2 - src/thermo/VPSSMgr_IdealGas.cpp | 1 - src/thermo/VPStandardStateTP.cpp | 1 - src/thermo/WaterProps.cpp | 43 ------------ src/thermo/WaterPropsIAPWS.cpp | 7 -- src/thermo/WaterPropsIAPWSphi.cpp | 1 - src/thermo/WaterSSTP.cpp | 4 -- src/tpx/CarbonDioxide.cpp | 3 - src/tpx/HFC134a.cpp | 13 ---- src/tpx/Hydrogen.cpp | 1 - src/tpx/Oxygen.cpp | 1 - src/tpx/RedlichKwong.cpp | 1 - src/tpx/Water.cpp | 15 ----- src/transport/AqueousTransport.cpp | 28 -------- src/transport/LiquidTranInteraction.cpp | 58 ----------------- src/transport/LiquidTransport.cpp | 53 --------------- src/transport/MultiTransport.cpp | 17 ----- src/transport/PecosTransport.cpp | 43 ------------ src/transport/SimpleTransport.cpp | 87 ------------------------- src/transport/SolidTransport.cpp | 9 --- src/transport/TransportFactory.cpp | 18 ----- src/zeroD/FlowReactor.cpp | 1 - src/zeroD/ReactorNet.cpp | 2 +- 101 files changed, 33 insertions(+), 1225 deletions(-) diff --git a/src/apps/cti2ctml.cpp b/src/apps/cti2ctml.cpp index 00a93dcde..9289a1514 100644 --- a/src/apps/cti2ctml.cpp +++ b/src/apps/cti2ctml.cpp @@ -62,9 +62,6 @@ int main(int argc, char** argv) XML_Node* xc = new XML_Node(); std::string path = findInputFile(infile); ctml::get_CTML_Tree(xc, path, 0); - //XML_Node *xd = new XML_Node(); - //xc->copy(xd); - } catch (CanteraError& err) { std::cout << err.what() << std::endl; } diff --git a/src/base/application.cpp b/src/base/application.cpp index 0a2b5fa98..18abccee7 100644 --- a/src/base/application.cpp +++ b/src/base/application.cpp @@ -167,7 +167,6 @@ Application::Application() : // install a default logwriter that writes to standard // output / standard error - // logwriter = new Logger(); setDefaultDirectories(); #if defined(THREAD_SAFE_CANTERA) Unit::units() ; diff --git a/src/base/global.cpp b/src/base/global.cpp index 517f98df9..732a1db38 100644 --- a/src/base/global.cpp +++ b/src/base/global.cpp @@ -104,7 +104,6 @@ void thread_complete() XML_Node* get_XML_File(const std::string& file, int debug) { XML_Node* xtmp = app()->get_XML_File(file, debug) ; - //writelog("get_XML_File: returned from app:get_XML_FILE " + int2str(xtmp) + "\n"); return xtmp; } diff --git a/src/base/stringUtils.cpp b/src/base/stringUtils.cpp index 0f88ef01a..2aff3031b 100644 --- a/src/base/stringUtils.cpp +++ b/src/base/stringUtils.cpp @@ -132,7 +132,6 @@ compositionMap parseCompString(const std::string& ss, s = s.substr(ibegin,s.size()); size_t icolon = s.find(':'); size_t iend = s.find_first_of(", ;\n\t"); - //icomma = s.find(','); if (icolon != std::string::npos) { std::string name = stripws(s.substr(0, icolon)); if (iend != std::string::npos) { diff --git a/src/base/xml.cpp b/src/base/xml.cpp index d53a3f9ef..d829dc1c8 100644 --- a/src/base/xml.cpp +++ b/src/base/xml.cpp @@ -1119,12 +1119,10 @@ XML_Node* findXMLPhase(XML_Node* root, } for (size_t n = 0; n < root->nChildren(); n++) { sc = vsc[n]; - //if (sc->name() != "phase") { - scResult = findXMLPhase(sc, idtarget); - if (scResult) { - return scResult; - } - //} + scResult = findXMLPhase(sc, idtarget); + if (scResult) { + return scResult; + } } return scResult; } diff --git a/src/clib/ct.cpp b/src/clib/ct.cpp index cb6d965d9..d999a596d 100644 --- a/src/clib/ct.cpp +++ b/src/clib/ct.cpp @@ -1394,7 +1394,6 @@ extern "C" { try { string s; writelog("function readlog is deprecated!"); - //getlog(s); int nlog = static_cast(s.size()); if (n < 0) { return nlog; @@ -1403,7 +1402,6 @@ extern "C" { copy(s.begin(), s.begin() + nn, buf); buf[min(nlog, n-1)] = '\0'; - //clearlog(); return 0; } catch (...) { return handleAllExceptions(-1, ERR); diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index 8a2ffa5b5..c9c61e43f 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -249,14 +249,6 @@ int ChemEquil::setInitialMoles(thermo_t& s, vector_fp& elMoleGoal, int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, vector_fp& elMolesGoal, int loglevel) { - //for (k = 0; k < m_kk; k++) { - // if (m_molefractions[k] > 0.0) { - // m_molefractions[k] = fmaxx(m_molefractions[k], 0.05); - // } - //} - //s.setState_PX(s.pressure(), m_molefractions.begin()); - - vector_fp b(m_mm, -999.0); vector_fp mu_RT(m_kk, 0.0); vector_fp xMF_est(m_kk, 0.0); @@ -646,12 +638,6 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, int info = estimateEP_Brinkley(s, x, elMolesGoal); if (info == 0) { setToEquilState(s, x, s.temperature()); - // Tempting -> However, nonideal is a problem. Turn on if not worried - // about nonideality and you are having problems with the main - // algorithm. - //if (XY == TP) { - // return 0; - //} } /* @@ -1554,7 +1540,6 @@ int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, "Jacobian is singular. \nTry adding more species, " "changing the elemental composition slightly, \nor removing " "unused elements."); - //return -3; } /* diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp index 62ded410f..0cef809ee 100644 --- a/src/equil/MultiPhase.cpp +++ b/src/equil/MultiPhase.cpp @@ -218,7 +218,6 @@ void MultiPhase::init() esum += m_atoms(m,k) * m_atomicNumber[m]; } } - //m_atoms(m_eloc, k) += esum; } } @@ -615,10 +614,6 @@ doublereal MultiPhase::equilibrate(int XY, doublereal err, // if 'strt' is false, the current composition will be used as // the starting estimate; otherwise it will be estimated - // if (e) { - // cout << "e should be zero, but it is not!" << endl; - // delete e; - // } e = new MultiPhaseEquil(this, strt); // start with a loose error tolerance, but tighten it as we get // close to the final temperature @@ -654,19 +649,17 @@ doublereal MultiPhase::equilibrate(int XY, doublereal err, } else { tnew = sqrt(Tlow*Thigh); dt = tnew - m_temp; - //cpb = cp(); } herr = fabs((h0 - hnow)/h0); - if (herr < err) { // || dta < 1.0e-4) { + if (herr < err) { goto done; } tnew = m_temp + dt; if (tnew < 0.0) { tnew = 0.5*m_temp; } - //dta = fabs(tnew - m_temp); setTemperature(tnew); // if the size of Delta T is not too large, use @@ -696,8 +689,8 @@ doublereal MultiPhase::equilibrate(int XY, doublereal err, "No convergence for T"); } else if (XY == SP) { s0 = entropy(); - Tlow = 1.0; // m_Tmin; // lower bound on T - Thigh = 1.0e6; // m_Tmax; // upper bound on T + Tlow = 1.0; // lower bound on T + Thigh = 1.0e6; // upper bound on T for (n = 0; n < maxiter; n++) { delete e; e = new MultiPhaseEquil(this, strt); @@ -749,7 +742,6 @@ doublereal MultiPhase::equilibrate(int XY, doublereal err, throw CanteraError("MultiPhase::equilibrate", "No convergence for T"); } else if (XY == TV) { - // doublereal dt = 1.0e3; doublereal v0 = volume(); doublereal dVdP; int n; diff --git a/src/equil/MultiPhaseEquil.cpp b/src/equil/MultiPhaseEquil.cpp index 54e2c36fb..277d7e2fc 100644 --- a/src/equil/MultiPhaseEquil.cpp +++ b/src/equil/MultiPhaseEquil.cpp @@ -16,9 +16,6 @@ namespace Cantera MultiPhaseEquil::MultiPhaseEquil(MultiPhase* mix, bool start, int loglevel) : m_mix(mix) { - // the multi-phase mixture - // m_mix = mix; - // store some mixture parameters locally m_nel_mix = mix->nElements(); m_nsp_mix = mix->nSpecies(); @@ -524,7 +521,7 @@ doublereal MultiPhaseEquil::stepComposition(int loglevel) if (m_work[k] < 0.0 && m_moles[k] > 0.0) { omax = -m_moles[k]/m_work[k]; if (omax < omegamax) { - omegamax = omax; //*1.000001; + omegamax = omax; if (omegamax < 1.0e-5) { m_force = true; } @@ -622,8 +619,7 @@ doublereal MultiPhaseEquil::computeReactionSteps(vector_fp& dxi) for (k = 0; k < m_nsp; k++) { kc = m_species[k]; if (m_mix->speciesPhaseIndex(kc) == ip) { - // bug fixed 7/12/06 DGG - stoich = nu[k]; // nu[kc]; + stoich = nu[k]; psum += stoich * stoich; } } @@ -758,7 +754,6 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) VolPM.resize(nSpecies, 0.0); tref.getMoleFractions(&mf[istart]); tref.getPartialMolarVolumes(DATA_PTR(VolPM)); - //vcs_VolPhase *volP = m_vprob->VPhaseList[iphase]; double TMolesPhase = phaseMoles(iphase); double VolPhaseVolumes = 0.0; @@ -773,8 +768,6 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) fprintf(FP,"Temperature = %11.5g kelvin\n", Temp); fprintf(FP,"Pressure = %11.5g Pascal\n", pres); fprintf(FP,"Total Volume = %11.5g m**3\n", vol); - // fprintf(FP,"Number Basis optimizations = %d\n", m_vprob->m_NumBasisOptimizations); - // fprintf(FP,"Number VCS iterations = %d\n", m_vprob->m_Iterations); for (size_t iphase = 0; iphase < nphase; iphase++) { istart = m_mix->speciesIndex(0, iphase); @@ -783,9 +776,7 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) ThermoPhase* tp = &tref; tp->getMoleFractions(&mf[istart]); string phaseName = tref.name(); - // vcs_VolPhase *volP = m_vprob->VPhaseList[iphase]; double TMolesPhase = phaseMoles(iphase); - //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase)); nSpecies = tref.nSpecies(); activity.resize(nSpecies, 0.0); ac.resize(nSpecies, 0.0); @@ -859,21 +850,6 @@ void MultiPhaseEquil::reportCSV(const std::string& reportFile) VolPM[k], VolPhaseVolumes); } } -#ifdef DEBUG_MODE - /* - * Check consistency: These should be equal - */ - tp->getChemPotentials(&(fe[istart])); - for (k = 0; k < nSpecies; k++) { - //if (!vcs_doubleEqual(fe[istart+k], mu[k])) { - // fprintf(FP,"ERROR: incompatibility!\n"); - // fclose(FP); - // printf("ERROR: incompatibility!\n"); - // exit(EXIT_FAILURE); - // } - } -#endif - } fclose(FP); } diff --git a/src/equil/vcs_MultiPhaseEquil.cpp b/src/equil/vcs_MultiPhaseEquil.cpp index 3de7d7017..83f4c9746 100644 --- a/src/equil/vcs_MultiPhaseEquil.cpp +++ b/src/equil/vcs_MultiPhaseEquil.cpp @@ -58,7 +58,6 @@ int vcs_MultiPhaseEquil::equilibrate_TV(int XY, doublereal xtarget, int printLvl, doublereal err, int maxsteps, int loglevel) { - // doublereal dt = 1.0e3; doublereal Vtarget = m_mix->volume(); if ((XY != TV) && (XY != HV) && (XY != UV) && (XY != SV)) { throw CanteraError("vcs_MultiPhaseEquil::equilibrate_TV", @@ -546,7 +545,6 @@ int vcs_MultiPhaseEquil::equilibrate_TP(int estimateEquil, for (size_t k = 0; k < tref.nSpecies(); k++, kGlob++) { phaseMole += m_vprob.w[kGlob]; } - //phaseMole *= 1.0E-3; m_mix->setPhaseMoles(ip, phaseMole); } @@ -1409,13 +1407,7 @@ int vcs_MultiPhaseEquil::determine_PhaseStability(int iph, double& funcStab, int * states. */ m_mix->uploadMoleFractionsFromPhases(); - // 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]); - //} double te = tickTock.secondsWC(); if (printLvl > 0) { diff --git a/src/equil/vcs_TP.cpp b/src/equil/vcs_TP.cpp index f9e125917..7212485cd 100644 --- a/src/equil/vcs_TP.cpp +++ b/src/equil/vcs_TP.cpp @@ -59,14 +59,6 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) { - // int i; - //double R; - /* - * At this level of the program, we are still using values - * for the free energies that have units. - */ - // R = vcsUtil_gasConstant(m_VCS_UnitsFormat); - /* * We need to special case VCS_UNITS_UNITLESS, here. * cpc_ts_GStar_calc() returns units of Kelvin. Also, the temperature @@ -78,19 +70,6 @@ int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) * them back to unitless at the end of this routine. */ - - /* - * Loop over the species calculating the standard state Gibbs free - * energies. -> These are energies that only depend upon the Temperature - * and possibly on the pressure (i.e., ideal gas, etc). - */ - // HKM -> We can change this to looks over phases, calling the vcs_VolPhase - // object. Working to get rid of VCS_SPECIES_THERMO object - //for (i = 0; i < m_numSpeciesTot; ++i) { - // VCS_SPECIES_THERMO *spt = SpeciesThermo[i]; - // ff[i] = R * spt->GStar_R_calc(i, Temp, pres); - //} - for (size_t iph = 0; iph < m_numPhases; iph++) { vcs_VolPhase* vph = m_VolPhaseList[iph]; vph->setState_TP(m_temperature, m_pressurePA); diff --git a/src/equil/vcs_VolPhase.cpp b/src/equil/vcs_VolPhase.cpp index 95c657017..5b2da4c38 100644 --- a/src/equil/vcs_VolPhase.cpp +++ b/src/equil/vcs_VolPhase.cpp @@ -148,7 +148,6 @@ vcs_VolPhase& vcs_VolPhase::operator=(const vcs_VolPhase& b) * Do a shallow copy because we haven' figured this out. */ IndSpecies = b.IndSpecies; - //IndSpeciesContig = b.IndSpeciesContig; for (size_t k = 0; k < old_num; k++) { if (ListSpeciesPtr[k]) { @@ -539,9 +538,6 @@ void vcs_VolPhase::setMolesFromVCS(const int stateCalc, // 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++) { - // Xmol_[k] = 1.0 / m_numSpecies; - //} m_existence = VCS_PHASE_EXIST_NO; } /* @@ -815,21 +811,6 @@ void vcs_VolPhase::_updateLnActCoeffJac() */ _updateMoleFractionDependencies(); _updateActCoeff(); - /* - * Calculate the column of the matrix - */ - double* const np_lnActCoeffCol = np_dLnActCoeffdMolNumber[j]; - for (size_t k = 0; k < m_numSpecies; k++) { - double tmp; - tmp = (ActCoeff[k] - ActCoeff_Base[k]) / - ((ActCoeff[k] + ActCoeff_Base[k]) * 0.5 * deltaMoles_j); - if (fabs(tmp - np_lnActCoeffCol[k]) > 1.0E-4 * fabs(tmp) + fabs(np_lnActCoeffCol[k])) { - // printf(" we have an error\n"); - - } - //tmp = lnActCoeffCol[k]; - - } /* * Revert to the base case Xmol_, v_totalMoles */ diff --git a/src/equil/vcs_equilibrate.cpp b/src/equil/vcs_equilibrate.cpp index 25085ac8f..60243a886 100644 --- a/src/equil/vcs_equilibrate.cpp +++ b/src/equil/vcs_equilibrate.cpp @@ -158,7 +158,6 @@ int vcs_equilibrate_1(MultiPhase& s, int ixy, } } else { throw CanteraError("equilibrate","unsupported option"); - //return -1.0; } } else { throw CanteraError("vcs_equilibrate_1", "unknown solver"); diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index f2d8ee5c9..53db9c5f2 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -25,8 +25,6 @@ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, size_t nspecies = m_numSpeciesTot; size_t nrxn = m_numRxnTot; - // double *molNum = VCS_DATA_PTR(m_molNumSpecies_old); - /* * CALL ROUTINE TO SOLVE MAX(CC*molNum) SUCH THAT AX*molNum = BB * AND molNum(I) .GE. 0.0 diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index 8e7c3de46..4c7f8cb0e 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -90,7 +90,6 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const * We loop through the regular reaction looking for a reaction that can pop the * component. */ - //printf("WE are here at new logic - CHECK\n"); for (size_t jrxn = 0; jrxn < m_numRxnRdc; jrxn++) { bool foundJrxn = false; // First, if the component is a product of the reaction @@ -105,7 +104,6 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const } } if (foundJrxn) { - //printf("We have found a component phase pop! CHECK1 \n"); return true; } } @@ -126,7 +124,6 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const } } if (foundJrxn) { - //printf("We have found a component phase pop! CHECK2 \n"); return true; } } diff --git a/src/equil/vcs_prob.cpp b/src/equil/vcs_prob.cpp index 520b200ac..bc47b717c 100644 --- a/src/equil/vcs_prob.cpp +++ b/src/equil/vcs_prob.cpp @@ -261,7 +261,6 @@ void VCS_PROB::prob_report(int print_lvl) plogf(" Target_kmol ElemType ElActive\n"); double fac = 1.0; if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - //fac = 1.0E3; fac = 1.0; } for (size_t i = 0; i < ne; ++i) { @@ -406,7 +405,6 @@ void VCS_PROB::reportCSV(const std::string& reportFile) for (size_t iphase = 0; iphase < NPhase; iphase++) { size_t istart = iK; vcs_VolPhase* volP = VPhaseList[iphase]; - //const Cantera::ThermoPhase *tptr = volP->ptrThermoPhase(); size_t nSpeciesPhase = volP->nSpecies(); volPM.resize(nSpeciesPhase, 0.0); volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); @@ -439,7 +437,6 @@ void VCS_PROB::reportCSV(const std::string& reportFile) size_t nSpeciesPhase = volP->nSpecies(); volP->sendToVCS_VolPM(VCS_DATA_PTR(volPM)); double TMolesPhase = volP->totalMoles(); - //AssertTrace(TMolesPhase == m_mix->phaseMoles(iphase)); activity.resize(nSpeciesPhase, 0.0); ac.resize(nSpeciesPhase, 0.0); diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index 5ac1e65cb..66d66f547 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -707,9 +707,6 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) if (dx_orig > 0.0) { dx = 0.0; #ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - //plogf(" --- %s :Warning possible error dx>0 dg > 0\n", SpName[kspec]); - } sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx>0 dg > 0"); #endif return dx; @@ -718,9 +715,6 @@ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) if (dx_orig < 0.0) { dx = 0.0; #ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - //plogf(" --- %s :Warning possible error dx<0 dg < 0\n", SpName[kspec]); - } sprintf(ANOTE, "Rxn reduced to zero step size in line search: dx<0 dg < 0"); #endif return dx; diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index ee55bc6d4..2d067d6a7 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -187,10 +187,6 @@ int VCS_SOLVE::vcs_setMolesLinProg() } } - // set the moles of the phase objects to match - // updateMixMoles(); - // Update the phase objects with the contents of the m_molNumSpecies_old vector - // vcs_updateVP(0); #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { printProgress(m_speciesName, m_molNumSpecies_old, m_SSfeSpecies); diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index c911734b2..d2b0c5078 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -91,7 +91,6 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, " Number of species is nonpositive\n"); } - //vcs_priv_init(this); m_VCS_UnitsFormat = VCS_UNITS_UNITLESS; /* diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 62f97c765..2c7cb9ea8 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -602,7 +602,6 @@ L_MAINLOOP_ALL_SPECIES: if (resurrect) { bool phaseResurrected = false; if (Vphase->exists() == VCS_PHASE_EXIST_NO) { - //Vphase->setExistence(1); phaseResurrected = true; } @@ -1913,18 +1912,6 @@ L_RETURN_BLOCK_B: * with this information. */ vcs_updateVP(VCS_STATECALC_OLD); - /* - * Store the final Delta G values for each non-component species - * in the species slot rather than the reaction slot - */ - // kspec = m_numSpeciesTot; - // i = m_numRxnTot; - //for (irxn = 0; irxn < m_numRxnTot; ++irxn) { - // --kspec; - // --i; - // m_deltaGRxn_new[kspec] = m_deltaGRxn_new[i]; - //} - // vcs_dzero(VCS_DATA_PTR(m_deltaGRxn_new), m_numComponents); /* * Evaluate the final mole fractions * storing them in wt[] @@ -2015,7 +2002,6 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete * get the diagonal of the activity coefficient jacobian */ s = m_np_dLnActCoeffdMolNum[kspec][kspec] / (m_tPhaseMoles_old[iph]); - // s *= (m_tPhaseMoles_old[iph]); /* * We fit it to a power law approximation of the activity coefficient * @@ -2043,7 +2029,6 @@ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete tmp = 200.; } wTrial = w_kspec * exp(tmp); - // wTrial = w_kspec * exp(-dg_irxn); molNum_kspec_new = wTrial; @@ -2260,7 +2245,6 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec) void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) { - // int irxn = kspec - m_numComponents; size_t iph = m_phaseID[kspec]; #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { @@ -3321,7 +3305,6 @@ L_END_LOOP: for (size_t j = 0; j < ncTrial; j++) { plogf(" %10.10s", m_speciesName[j].c_str()); } - //plogf("| m_scSize"); plogf("\n"); for (size_t i = 0; i < m_numRxnTot; i++) { plogf(" --- %3d ", m_indexRxnToSpecies[i]); @@ -3334,7 +3317,6 @@ L_END_LOOP: for (size_t j = 0; j < ncTrial; j++) { plogf(" %+7.3f", m_stoichCoeffRxnMatrix[i][j]); } - //plogf(" | %6.2f", m_scSize[i]); plogf("\n"); } @@ -4162,8 +4144,6 @@ double VCS_SOLVE::vcs_tmoles() for (size_t i = 0; i < m_numPhases; i++) { sum += m_tPhaseMoles_old[i]; vcs_VolPhase* Vphase = m_VolPhaseList[i]; - // Took out because we aren't updating mole fractions in Vphase - // Vphase->TMoles = m_tPhaseMoles_old[i]; if (m_tPhaseMoles_old[i] == 0.0) { Vphase->setTotalMoles(0.0); } else { @@ -4469,7 +4449,6 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, * This can probably be solved by successive iteration. * This should be implemented. */ - //alterZeroedPhases = false; if (alterZeroedPhases && false) { for (size_t iph = 0; iph < m_numPhases; iph++) { bool lneed = false; @@ -4560,7 +4539,6 @@ void VCS_SOLVE::vcs_printDeltaG(const int stateCalc) for (size_t j = 0; j < m_numComponents; j++) { plogf("%-10.10s", m_speciesName[j].c_str()); } - //plogf("| m_scSize"); plogf("\n"); for (size_t i = 0; i < m_numRxnTot; i++) { plogf(" --- %3d ", m_indexRxnToSpecies[i]); @@ -4573,7 +4551,6 @@ void VCS_SOLVE::vcs_printDeltaG(const int stateCalc) for (size_t j = 0; j < m_numComponents; j++) { plogf(" %6.2f", m_stoichCoeffRxnMatrix[i][j]); } - //plogf(" | %6.2f", m_scSize[i]); plogf("\n"); } plogf(" "); @@ -4843,8 +4820,6 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k #endif pv1->setSpGlobalIndexVCS(kp1, k2); pv2->setSpGlobalIndexVCS(kp2, k1); - //pv1->IndSpecies[kp1] = k2; - //pv2->IndSpecies[kp2] = k1; std::swap(m_speciesName[k1], m_speciesName[k2]); std::swap(m_molNumSpecies_old[k1], m_molNumSpecies_old[k2]); std::swap(m_speciesUnknownType[k1], m_speciesUnknownType[k2]); @@ -4904,13 +4879,6 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k std::swap(m_deltaGRxn_new[i1], m_deltaGRxn_new[i2]); std::swap(m_deltaGRxn_old[i1], m_deltaGRxn_old[i2]); std::swap(m_deltaGRxn_tmp[i1], m_deltaGRxn_tmp[i2]); - - /* - * We don't want to swap ir[], because the values of ir should - * stay the same after the swap - * - * vcs_isw(ir, i1, i2); - */ } } diff --git a/src/equil/vcs_solve_phaseStability.cpp b/src/equil/vcs_solve_phaseStability.cpp index 792c493b1..b427f78d3 100644 --- a/src/equil/vcs_solve_phaseStability.cpp +++ b/src/equil/vcs_solve_phaseStability.cpp @@ -134,21 +134,6 @@ int VCS_SOLVE::vcs_PS(VCS_PROB* vprob, int iphase, int printLvl, double& feStabl */ vcs_redim_TP(); - /* - vcs_VolPhase *Vphase = m_VolPhaseList[iphase]; - - std::vector mfPop = Vphase->moleFractions(); - int nsp = Vphase->nSpecies(); - - vcs_VolPhase *VPphase = vprob->VPhaseList[iphase]; - int kstart = Vphase->spGlobalIndexVCS(0); - for (int k = 0; k < nsp; k++) { - vprob->mf[kstart + k] = mfPop[k]; - } - VPphase->setMoleFractionsState(Vphase->totalMoles(), - VCS_DATA_PTR(Vphase->moleFractions()), - VCS_STATECALC_TMP); - */ vcs_prob_update(vprob); /* * Return the convergence success flag. @@ -162,7 +147,6 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, { double test = -1.0E-10; bool usedZeroedSpecies; - // std::vector phasePopPhaseIDs(0); int iStab = 0; std::vector sm(m_numElemConstraints*m_numElemConstraints, 0.0); @@ -187,8 +171,6 @@ int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, vcs_printDeltaG(VCS_STATECALC_OLD); } vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_Deficient), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); - // phasePopPhaseIDs.clear(); - // vcs_popPhaseID(phasePopPhaseIDs); funcVal = vcs_phaseStabilityTest(iph); if (funcVal > 0.0) { iStab = 1; diff --git a/src/fortran/fctxml.cpp b/src/fortran/fctxml.cpp index 230cf71a3..c23bbb044 100644 --- a/src/fortran/fctxml.cpp +++ b/src/fortran/fctxml.cpp @@ -310,7 +310,6 @@ extern "C" { for (int i = 0; i < nv; i++) { data[i] = v[i]; } - //n = nv; } catch (...) { return handleAllExceptions(-1, ERR); } diff --git a/src/kinetics/AqueousKinetics.cpp b/src/kinetics/AqueousKinetics.cpp index f9f0329a5..a870c6e35 100644 --- a/src/kinetics/AqueousKinetics.cpp +++ b/src/kinetics/AqueousKinetics.cpp @@ -140,7 +140,6 @@ void AqueousKinetics::updateKc() // compute Delta G^0 for all reversible reactions m_rxnstoich.getRevReactionDelta(m_ii, &m_grt[0], &m_rkcn[0]); - //doublereal logStandConc = m_kdata->m_logStandConc; doublereal rrt = 1.0/(GasConstant * thermo().temperature()); for (size_t i = 0; i < m_nrev; i++) { size_t irxn = m_revindex[i]; diff --git a/src/kinetics/InterfaceKinetics.cpp b/src/kinetics/InterfaceKinetics.cpp index 7d9d14566..6b9b908f7 100644 --- a/src/kinetics/InterfaceKinetics.cpp +++ b/src/kinetics/InterfaceKinetics.cpp @@ -232,7 +232,6 @@ void InterfaceKinetics::updateKc() { fill(m_rkcn.begin(), m_rkcn.end(), 0.0); - //static vector_fp mu(nTotalSpecies()); if (m_nrev > 0) { /* * Get the vector of standard state electrochemical potentials for species in the Interfacial @@ -278,15 +277,12 @@ void InterfaceKinetics::checkPartialEquil() doublereal rt = GasConstant*thermo(0).temperature(); cout << "T = " << thermo(0).temperature() << " " << rt << endl; size_t nsp, ik=0; - //doublereal rt = GasConstant*thermo(0).temperature(); - // doublereal rrt = 1.0/rt; doublereal delta; for (size_t n = 0; n < nPhases(); n++) { thermo(n).getChemPotentials(DATA_PTR(dmu) + m_start[n]); nsp = thermo(n).nSpecies(); for (size_t k = 0; k < nsp; k++) { delta = Faraday * m_phi[n] * thermo(n).charge(k); - //cout << thermo(n).speciesName(k) << " " << (delta+dmu[ik])/rt << " " << dmu[ik]/rt << endl; dmu[ik] += delta; ik++; } @@ -431,7 +427,6 @@ void InterfaceKinetics::applyButlerVolmerCorrection(doublereal* const kf) for (size_t i = 0; i < m_beta.size(); i++) { size_t irxn = m_ctrxn[i]; eamod = m_beta[i]*m_rwork[irxn]; - // if (eamod != 0.0 && m_E[irxn] != 0.0) { if (eamod != 0.0) { #ifdef DEBUG_KIN_MODE ea = GasConstant * m_E[irxn]; @@ -529,16 +524,11 @@ void InterfaceKinetics::updateROP() // multiply ropf by concentration products m_rxnstoich.multiplyReactants(DATA_PTR(m_conc), DATA_PTR(m_ropf)); - //m_reactantStoich.multiply(m_conc.begin(), ropf.begin()); // for reversible reactions, multiply ropr by concentration // products m_rxnstoich.multiplyRevProducts(DATA_PTR(m_conc), DATA_PTR(m_ropr)); - //m_revProductStoich.multiply(m_conc.begin(), ropr.begin()); - - // do global reactions - //m_globalReactantStoich.power(m_conc.begin(), ropf.begin()); for (size_t j = 0; j != m_ii; ++j) { m_ropnet[j] = m_ropf[j] - m_ropr[j]; @@ -617,9 +607,6 @@ void InterfaceKinetics::getDeltaGibbs(doublereal* deltaG) for (size_t n = 0; n < nPhases(); n++) { thermo(n).getChemPotentials(DATA_PTR(m_grt) + m_start[n]); } - //for (n = 0; n < m_grt.size(); n++) { - // cout << n << "G_RT = " << m_grt[n] << endl; - //} /* * Use the stoichiometric manager to find deltaG for each @@ -843,33 +830,6 @@ void InterfaceKinetics::setIOFlag(int ioFlag) } } -// void InterfaceKinetics:: -// addGlobalReaction(const ReactionData& r) { - -// int iloc; -// // install rate coeff calculator -// vector_fp rp = r.rateCoeffParameters; -// int ncov = r.cov.size(); -// for (int m = 0; m < ncov; m++) rp.push_back(r.cov[m]); -// iloc = m_rates.install( reactionNumber(), -// r.rateCoeffType, rp.size(), -// rp.begin() ); -// // store activation energy -// m_E.push_back(r.rateCoeffParameters[2]); -// // add constant term to rate coeff value vector -// m_rfn.push_back(r.rateCoeffParameters[0]); - -// int nr = r.order.size(); -// vector_fp ordr(nr); -// for (int n = 0; n < nr; n++) { -// ordr[n] = r.order[n] - r.rstoich[n]; -// } -// m_globalReactantStoich.add( reactionNumber(), -// r.reactants, ordr); - -// registerReaction( reactionNumber(), GLOBAL_RXN, iloc); -// } - void InterfaceKinetics::installReagents(const ReactionData& r) { diff --git a/src/kinetics/ReactionPath.cpp b/src/kinetics/ReactionPath.cpp index b41d18406..6af6a0c38 100644 --- a/src/kinetics/ReactionPath.cpp +++ b/src/kinetics/ReactionPath.cpp @@ -148,12 +148,6 @@ vector_int ReactionPathDiagram::reactions() void ReactionPathDiagram::add(ReactionPathDiagram& d) { - // doublereal f1, f2; - // int nnodes = nNodes(); - // if (nnodes != d.nNodes()) { - // throw CanteraError("ReactionPathDiagram::add", - // "number of nodes must be the same"); - // } size_t np = nPaths(); size_t n, k1, k2; Path* p = 0; @@ -210,10 +204,8 @@ void ReactionPathDiagram::writeData(ostream& s) k2 = m_speciesNumber[i2]; f1 = flow(k1, k2); f2 = flow(k2, k1); - //if (f1 > 0.001 || f2 > 0.001) { s << m_nodes[k1]->name << " " << m_nodes[k2]->name << " " << f1 << " " << -f2 << endl; - //} } } } @@ -221,7 +213,6 @@ void ReactionPathDiagram::writeData(ostream& s) void ReactionPathDiagram::exportToDot(ostream& s) { doublereal flxratio, flmax = 0.0, lwidth; - //s.flags(std::ios_base::showpoint+std::ios_base::fixed); s.precision(3); // a directed graph @@ -238,7 +229,6 @@ void ReactionPathDiagram::exportToDot(ostream& s) << endl; } - //s << "color = white;" << endl; if (dot_options != "") { s << dot_options << endl; } @@ -380,17 +370,16 @@ void ReactionPathDiagram::exportToDot(ostream& s) lwidth = 1.0 - 4.0 * log10(flxratio/threshold)/log10(threshold) + 1.0; s << "[fontname=\""+m_font+"\", style=\"setlinewidth(" - //<< 1.0 - arrow_width*flxratio << lwidth << ")\""; s << ", arrowsize=" - << std::min(6.0, 0.5*lwidth); // 1 - arrow_width*flxratio; + << std::min(6.0, 0.5*lwidth); } else { s << ", style=\"setlinewidth(" << arrow_width << ")\""; s << ", arrowsize=" << flxratio + 1; } - doublereal hue = 0.7; //2.0/(1.0 + pow(log10(flxratio),2)) ; + doublereal hue = 0.7; doublereal bright = 0.9; s << ", color=" << "\"" << hue << ", " << flxratio + 0.5 << ", " << bright << "\"" << endl; @@ -412,12 +401,11 @@ void ReactionPathDiagram::exportToDot(ostream& s) for (; b != m_nodes.end(); ++b) { if (b->second->visible) { s << "s" << b->first << " [ fontname=\""+m_font+"\", label=\"" << b->second->name - //<< " \\n " << b->second->value << "\"];" << endl; } } s << " label = " << "\"" << "Scale = " - << flmax << "\\l " << title << "\";" << endl; //created with Cantera (www.cantera.org)\\l\";" + << flmax << "\\l " << title << "\";" << endl; s << " fontname = \""+m_font+"\";" << endl << "}" << endl; } @@ -675,24 +663,14 @@ void ReactionPathBuilder::findElements(Kinetics& kin) int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) { - //m_warn.clear(); m_transfer.clear(); - - //const Kinetics::thermo_t& ph = kin.thermo(); - m_elementSymbols.clear(); findElements(kin); - //m_nel = ph.nElements(); - m_ns = kin.nTotalSpecies(); //ph.nSpecies(); + m_ns = kin.nTotalSpecies(); m_nr = kin.nReactions(); - //for (m = 0; m < m_nel; m++) { - // m_elementSymbols.push_back(ph.elementName(m)); - //} - // all reactants / products, even ones appearing on both sides // of the reaction - // mod 8/18/01 dgg vector > allProducts; vector > allReactants; for (size_t i = 0; i < m_nr; i++) { @@ -749,7 +727,6 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) } size_t nrnet = m_reac[i].size(); - // int npnet = m_prod[i].size(); // compute number of atoms of each element in each reaction, // excluding molecules that appear on both sides of the @@ -759,7 +736,7 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) for (n = 0; n < nrnet; n++) { k = m_reac[i][n]; for (size_t m = 0; m < m_nel; m++) { - m_elatoms(m,i) += m_atoms(k,m); //ph.nAtoms(k,m); + m_elatoms(m,i) += m_atoms(k,m); } } } @@ -769,7 +746,7 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) m_sgroup.resize(m_ns); for (size_t j = 0; j < m_ns; j++) { for (size_t m = 0; m < m_nel; m++) { - comp[m] = int(m_atoms(j,m)); //ph.nAtoms(j,m)); + comp[m] = int(m_atoms(j,m)); } m_sgroup[j] = Group(comp); } @@ -792,7 +769,6 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) nar = 0; nap = 0; for (size_t j = 0; j < nr; j++) { - // if (ph.nAtoms(m_reac[i][j],m) > 0) nar++; if (m_atoms(m_reac[i][j],m) > 0) { nar++; } @@ -816,8 +792,6 @@ int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin) string reactionLabel(size_t i, size_t kr, size_t nr, const std::vector& slist, const Kinetics& s) { - - //int np = s.nPhases(); string label = ""; for (size_t l = 0; l < nr; l++) { if (l != kr) { @@ -842,8 +816,7 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, doublereal threshold = 0.0; bool fwd_incl, rev_incl, force_incl; - // const Kinetics::thermo_t& ph = s.thermo(); - size_t m = m_enamemap[element]-1; //ph.elementIndex(element); + size_t m = m_enamemap[element]-1; r.element = element; if (m == npos) { @@ -853,14 +826,6 @@ int ReactionPathBuilder::build(Kinetics& s, const string& element, s.getFwdRatesOfProgress(DATA_PTR(m_ropf)); s.getRevRatesOfProgress(DATA_PTR(m_ropr)); - //ph.getMoleFractions(m_x.begin()); - - //doublereal sum = 0.0; - //for (k = 0; k < kk; k++) { - // sum += m_x[k] * ph.nAtoms(k,m); - //} - //sum *= ph.molarDensity(); - // species explicitly included or excluded vector& in_nodes = r.included(); vector& out_nodes = r.excluded(); diff --git a/src/kinetics/importKinetics.cpp b/src/kinetics/importKinetics.cpp index 0eda3c5b5..47c24b166 100644 --- a/src/kinetics/importKinetics.cpp +++ b/src/kinetics/importKinetics.cpp @@ -91,7 +91,6 @@ void checkRxnElementBalance(Kinetics& kin, bal.clear(); balp.clear(); balr.clear(); - //cout << "checking " << rdata.equation << endl; size_t np = rdata.products.size(); // iterate over the products @@ -104,22 +103,17 @@ void checkRxnElementBalance(Kinetics& kin, for (size_t m = 0; m < ph.nElements(); m++) { bal[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); balp[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); - //cout << "product species " << ph.speciesName(klocal) << " has " << ph.nAtoms(klocal,m) - // << " atoms of " << ph.elementName(m) << " and kstoich = " << kstoich << endl; } } for (size_t index = 0; index < rdata.reactants.size(); index++) { size_t kr = rdata.reactants[index]; size_t n = kin.speciesPhaseIndex(kr); - //klocal = kr - kin.start(n); size_t klocal = kr - kin.kineticsSpeciesIndex(0,n); kstoich = rdata.rstoich[index]; const ThermoPhase& ph = kin.speciesPhase(kr); for (size_t m = 0; m < ph.nElements(); m++) { bal[ph.elementName(m)] -= kstoich*ph.nAtoms(klocal,m); balr[ph.elementName(m)] += kstoich*ph.nAtoms(klocal,m); - //cout << "reactant species " << ph.speciesName(klocal) << " has " << ph.nAtoms(klocal,m) - // << " atoms of " << ph.elementName(m) << " and kstoich = " << kstoich << endl; } } @@ -211,7 +205,6 @@ bool getReagents(const XML_Node& rxn, Kinetics& kin, int rp, stoich.push_back(stch); ord = doublereal(stch); order.push_back(ord); - //cout << key[n] << " " << isp << " " << stch << endl; /* * Needed to process reaction orders below. @@ -452,11 +445,11 @@ static void getEfficiencies(const XML_Node& eff, Kinetics& kin, getPairs(eff, key, val); string nm; string phse = kin.thermo(0).id(); - for (size_t n = 0; n < key.size(); n++) { // ; bb != ee; ++bb) { - nm = key[n];// bb->first; + for (size_t n = 0; n < key.size(); n++) { + nm = key[n]; size_t k = kin.kineticsSpeciesIndex(nm, phse); if (k != npos) { - rdata.thirdBodyEfficiencies[k] = fpValue(val[n]); // bb->second; + rdata.thirdBodyEfficiencies[k] = fpValue(val[n]); } else if (!rules.skipUndeclaredThirdBodies) { throw CanteraError("getEfficiencies", "Encountered third-body " "efficiency for undefined species \"" + nm + "\"\n" @@ -1018,7 +1011,6 @@ bool buildSolutionFromXML(XML_Node& root, const std::string& id, { XML_Node* x; x = get_XML_NameID(nm, string("#")+id, &root); - // x = get_XML_Node(string("#")+id, &root); if (!x) { return false; } diff --git a/src/kinetics/solveSP.cpp b/src/kinetics/solveSP.cpp index b542764c9..712e1972e 100644 --- a/src/kinetics/solveSP.cpp +++ b/src/kinetics/solveSP.cpp @@ -185,10 +185,6 @@ int solveSP::solveSurfProb(int ifunc, doublereal time_scale, doublereal TKelvin, // Calculate the largest species in each phase evalSurfLarge(DATA_PTR(m_CSolnSP)); - /* - * Get the net production rate of all species in the kinetics manager. - */ - // m_kin->getNetProductionRates(DATA_PTR(m_netProductionRatesSave)); if (m_ioflag) { print_header(m_ioflag, ifunc, time_scale, true, reltol, abstol, @@ -531,8 +527,6 @@ void solveSP::fun_eval(doublereal* resid, const doublereal* CSoln, kindexSP = m_numTotSurfSpecies; for (isp = 0; isp < m_numBulkPhasesSS; isp++) { doublereal* XBlk = DATA_PTR(m_numEqn1); - //ThermoPhase *THptr = m_bulkPhasePtrs[isp]; - //THptr->getMoleFractions(XBlk); nsp = m_nSpeciesSurfPhase[isp]; size_t surfPhaseIndex = m_indexKinObjSurfPhase[isp]; InterfaceKinetics* m_kin = m_objects[isp]; diff --git a/src/matlab/ctmethods.cpp b/src/matlab/ctmethods.cpp index efbf2bf19..f32b0e046 100644 --- a/src/matlab/ctmethods.cpp +++ b/src/matlab/ctmethods.cpp @@ -89,8 +89,6 @@ extern "C" { void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { - // mexPrintf("Number of lhs = %d\n", nlhs); - // mexPrintf("number of rhs = %d\n", nrhs); // create a log writer for error messages if this is the // first MATLAB function call initLogger(); diff --git a/src/matlab/flowdevicemethods.cpp b/src/matlab/flowdevicemethods.cpp index 722d8f6cd..46ecf7094 100644 --- a/src/matlab/flowdevicemethods.cpp +++ b/src/matlab/flowdevicemethods.cpp @@ -2,8 +2,6 @@ #include "clib/ctreactor.h" #include "clib/ct.h" -//const double Undef = -999.123; - void flowdevicemethods(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { diff --git a/src/matlab/reactornetmethods.cpp b/src/matlab/reactornetmethods.cpp index c6101bc10..aad6e235e 100644 --- a/src/matlab/reactornetmethods.cpp +++ b/src/matlab/reactornetmethods.cpp @@ -6,8 +6,6 @@ #include "clib/ct.h" #include "ctmatutils.h" -//const double Undef = -999.123; - void reactornetmethods(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { diff --git a/src/matlab/wallmethods.cpp b/src/matlab/wallmethods.cpp index b6ed83965..3aed32132 100644 --- a/src/matlab/wallmethods.cpp +++ b/src/matlab/wallmethods.cpp @@ -6,8 +6,6 @@ #include "clib/ct.h" #include "ctmatutils.h" -//const double Undef = -999.123; - void wallmethods(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]) { diff --git a/src/numerics/BEulerInt.cpp b/src/numerics/BEulerInt.cpp index 8cb1ec388..1075ba702 100644 --- a/src/numerics/BEulerInt.cpp +++ b/src/numerics/BEulerInt.cpp @@ -261,12 +261,6 @@ void BEulerInt::setSolnWeights() { int i; if (m_itol == 1) { - /* - * Adjust the atol vector if we are using vector - * atol conditions. - */ - // m_func->adjustAtol(m_abstol); - for (i = 0; i < m_neq; i++) { m_ewt[i] = m_abstol[i] + m_reltol * 0.5 * (fabs(m_y_n[i]) + fabs(m_y_pred_n[i])); @@ -532,10 +526,6 @@ void BEulerInt::beuler_jac(GeneralMatrix& J, double* const f, * Call the function to get a jacobian. */ m_func->evalJacobian(time_curr, delta_t_n, CJ, y, ydot, J, f); -#ifdef DEBUG_HKM - //double dddd = J(89, 89); - //checkFinite(dddd); -#endif m_nJacEval++; m_nfe++; } else { @@ -595,7 +585,6 @@ void BEulerInt::beuler_jac(GeneralMatrix& J, double* const f, col_j = (double*) J.ptrColumn(j); ysave = y[j]; dy = dyVector[j]; - //dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7); y[j] = ysave + dy; dy = y[j] - ysave; @@ -613,7 +602,6 @@ void BEulerInt::beuler_jac(GeneralMatrix& J, double* const f, for (i = 0; i < m_neq; i++) { diff = subtractRD(m_wksp[i], f[i]); col_j[i] = diff / dy; - //col_j[i] = (m_wksp[i] - f[i])/dy; } y[j] = ysave; @@ -1665,11 +1653,6 @@ int BEulerInt::dampStep(double time_curr, const double* y0, */ for (j = 0; j < m_neq; j++) { y1[j] = y0[j] + ff*step0[j]; - // HKM setting intermediate y's to zero was a tossup. - // slightly different, equivalent results - //#ifdef DEBUG_HKM - // y1[j] = MAX(0.0, y1[j]); - //#endif } calc_ydot(m_order, y1, ydot1); @@ -1866,24 +1849,6 @@ int BEulerInt::solve_nonlinear_problem(double* const y_comm, } } - // If we are converged, then let's use the best solution possible - // for an end result. We did a resolve in dampStep(). Let's update - // the solution to reflect that. - // HKM 5/16 -> Took this out, since if the last step was a - // damped step, then adding stp1[j] is undamped, and - // may lead to oscillations. It kind of defeats the - // purpose of dampStep() anyway. - // if (m == 1) { - // for (int j = 0; j < m_neq; j++) { - // y_new[j] += stp1[j]; - // HKM setting intermediate y's to zero was a tossup. - // slightly different, equivalent results - // #ifdef DEBUG_HKM - // y_new[j] = MAX(0.0, y_new[j]); - // #endif - // } - // } - bool m_filterIntermediate = false; if (m_filterIntermediate) { if (m == 0) { diff --git a/src/numerics/BandMatrix.cpp b/src/numerics/BandMatrix.cpp index 2e52eac75..599acce64 100644 --- a/src/numerics/BandMatrix.cpp +++ b/src/numerics/BandMatrix.cpp @@ -338,7 +338,6 @@ doublereal BandMatrix::rcond(doublereal a1norm) throw CanteraError("BandMatrix::rcond()", "matrix isn't factored correctly"); } - // doublereal anorm = oneNorm(); size_t ldab = (2 *m_kl + m_ku + 1); int rinfo = 0; rcond = ct_dgbcon('1', m_n, m_kl, m_ku, DATA_PTR(ludata), ldab, DATA_PTR(m_ipiv), a1norm, DATA_PTR(work_), diff --git a/src/numerics/CVodeInt.cpp b/src/numerics/CVodeInt.cpp index d3b68b2d5..100614f3f 100644 --- a/src/numerics/CVodeInt.cpp +++ b/src/numerics/CVodeInt.cpp @@ -255,8 +255,6 @@ void CVodeInt::reinitialize(double t0, FuncEval& func) m_iopt[MAXORD] = m_maxord; m_ropt[HMAX] = m_hmax; - //if (m_cvode_mem) CVodeFree(m_cvode_mem); - // pass a pointer to func in m_data m_data = (void*)&func; int result; diff --git a/src/numerics/CVodesIntegrator.cpp b/src/numerics/CVodesIntegrator.cpp index f86e6e40f..39338740c 100644 --- a/src/numerics/CVodesIntegrator.cpp +++ b/src/numerics/CVodesIntegrator.cpp @@ -111,9 +111,6 @@ CVodesIntegrator::CVodesIntegrator() : m_mupper(0), m_mlower(0), m_sens_ok(false) { - //m_ropt.resize(OPT_SIZE,0.0); - //m_iopt = new long[OPT_SIZE]; - //fill(m_iopt, m_iopt+OPT_SIZE,0); } CVodesIntegrator::~CVodesIntegrator() @@ -131,8 +128,6 @@ CVodesIntegrator::~CVodesIntegrator() N_VDestroy_Serial(m_abstol); } delete m_fdata; - - //delete[] m_iopt; } double& CVodesIntegrator::solution(size_t k) @@ -327,7 +322,6 @@ void CVodesIntegrator::initialize(double t0, FuncEval& func) delete m_fdata; m_fdata = new FuncData(&func, func.nparams()); - //m_data = (void*)&func; flag = CVodeSetUserData(m_cvode_mem, (void*)m_fdata); if (flag != CV_SUCCESS) { throw CVodesErr("CVodeSetUserData failed."); @@ -345,13 +339,7 @@ void CVodesIntegrator::reinitialize(double t0, FuncEval& func) { m_t0 = t0; m_time = t0; - //try { func.getInitialConditions(m_t0, m_neq, NV_DATA_S(m_y)); - //} - //catch (CanteraError) { - //showErrors(); - //error("Teminating execution"); - //} int result; @@ -432,7 +420,6 @@ int CVodesIntegrator::nEvals() const long int ne; CVodeGetNumRhsEvals(m_cvode_mem, &ne); return ne; - //return m_iopt[NFE]; } double CVodesIntegrator::sensitivity(size_t k, size_t p) diff --git a/src/numerics/DenseMatrix.cpp b/src/numerics/DenseMatrix.cpp index dd13191ed..7f5e93d81 100644 --- a/src/numerics/DenseMatrix.cpp +++ b/src/numerics/DenseMatrix.cpp @@ -94,7 +94,7 @@ void DenseMatrix::mult(const double* b, double* prod) const { ct_dgemv(ctlapack::ColMajor, ctlapack::NoTranspose, static_cast(nRows()), - static_cast(nRows()), 1.0, ptrColumn(0), //begin(), + static_cast(nRows()), 1.0, ptrColumn(0), static_cast(nRows()), b, 1, 0.0, prod, 1); } diff --git a/src/numerics/Func1.cpp b/src/numerics/Func1.cpp index 9b7395fad..73fbe1554 100644 --- a/src/numerics/Func1.cpp +++ b/src/numerics/Func1.cpp @@ -239,20 +239,17 @@ string Pow1::write(const std::string& arg) const string Const1::write(const std::string& arg) const { - //cout << "Const1" << endl; return fp2str(m_c); } string Ratio1::write(const std::string& arg) const { - //cout << "Ratio1" << endl; return "\\frac{" + m_f1->write(arg) + "}{" + m_f2->write(arg) + "}"; } string Product1::write(const std::string& arg) const { - //cout << "Product1" << endl; string s = m_f1->write(arg); if (m_f1->order() < order()) { s = "\\left(" + s + "\\right)"; @@ -266,7 +263,6 @@ string Product1::write(const std::string& arg) const string Sum1::write(const std::string& arg) const { - //cout << "Sum1" << endl; string s1 = m_f1->write(arg); string s2 = m_f2->write(arg); if (s2[0] == '-') { @@ -278,7 +274,6 @@ string Sum1::write(const std::string& arg) const string Diff1::write(const std::string& arg) const { - //cout << "Diff1" << endl; string s1 = m_f1->write(arg); string s2 = m_f2->write(arg); if (s2[0] == '-') { @@ -290,14 +285,12 @@ string Diff1::write(const std::string& arg) const string Composite1::write(const std::string& arg) const { - //cout << "Composite1" << endl; string g = m_f2->write(arg); return m_f1->write(g); } string TimesConstant1::write(const std::string& arg) const { - //cout << "TimesConstant1" << endl; string s = m_f1->write(arg); if (m_f1->order() < order()) { s = "\\left(" + s + "\\right)"; @@ -317,7 +310,6 @@ string TimesConstant1::write(const std::string& arg) const string PlusConstant1::write(const std::string& arg) const { - //cout << "PlusConstant1" << endl; if (m_c == 0.0) { return m_f1->write(arg); } diff --git a/src/numerics/IDA_Solver.cpp b/src/numerics/IDA_Solver.cpp index 2c686d802..3aaebbb8e 100644 --- a/src/numerics/IDA_Solver.cpp +++ b/src/numerics/IDA_Solver.cpp @@ -109,7 +109,6 @@ extern "C" { doublereal* const* colPts = Jac->cols; Cantera::IDA_Solver* s = d->m_solver; double delta_t = s->getCurrentStepFromIDA(); - // printf(" delta_t = %g 1/cj = %g\n", delta_t, 1.0/c_j); f->evalJacobianDP(t, delta_t, c_j, ydata, ydotdata, colPts, rdata); return 0; } diff --git a/src/numerics/NonlinearSolver.cpp b/src/numerics/NonlinearSolver.cpp index d5f635986..a8bd34ac8 100644 --- a/src/numerics/NonlinearSolver.cpp +++ b/src/numerics/NonlinearSolver.cpp @@ -188,8 +188,6 @@ NonlinearSolver::NonlinearSolver(ResidJacEval* func) : m_ewt[i] = atolk_[i]; } - - // jacCopyPtr_->resize(neq_, 0.0); deltaX_CP_.resize(neq_, 0.0); Jd_.resize(neq_, 0.0); deltaX_trust_.resize(neq_, 1.0); @@ -1133,7 +1131,6 @@ int NonlinearSolver::doAffineNewtonSolve(const doublereal* const y_curr, const return info; } - // doublereal *JTF = delta_y; vector_fp delyH(neq_); // First recalculate the scaled residual. It got wiped out doing the newton solve if (m_rowScaling) { @@ -1329,9 +1326,6 @@ doublereal NonlinearSolver::doCauchyPointSolve(GeneralMatrix& jac) RJd_norm_ += m_resid[i] * Jd_[i] / m_residWts[i]; JdJd_norm_ += Jd_[i] * Jd_[i]; } - //if (RJd_norm_ > -1.0E-300) { - // printf("we are here: zero residual\n"); - //} if (fabs(JdJd_norm_) < 1.0E-290) { if (fabs(RJd_norm_) < 1.0E-300) { lambdaStar_ = 0.0; @@ -1529,16 +1523,6 @@ void NonlinearSolver::setupDoubleDogleg() * --------------------------------------------- * (grad f)T H (grad f) (grad f)T H-1 (grad f) */ - // doublereal sumG = 0.0; - // doublereal sumH = 0.0; - // for (int i = 0; i < neq_; i++) { - // sumG = deltax_cp_[i] * deltax_cp_[i]; - // sumH = deltax_cp_[i] * newtDir[i]; - // } - // double fac1 = sumG / lambdaStar_; - // double fac2 = sumH / lambdaStar_; - // double gamma = fac1 / fac2; - // doublereal gamma = m_normDeltaSoln_CP / m_normDeltaSoln_Newton; /* * This hasn't worked. so will do it heuristically. One issue is that the newton * direction is not the inverse of the Hessian times the gradient. The Hession @@ -1944,8 +1928,6 @@ void NonlinearSolver::readjustTrustVector() doublereal deltaXSizeOld = trustNorm; doublereal trustNormGoal = trustNorm * trustDelta_; - // This is the size of each component. - // doublereal trustDeltaEach = trustDelta_ * trustNorm / neq_; doublereal oldVal; doublereal fabsy; // we use the old value of the trust region as an indicator @@ -1954,7 +1936,6 @@ void NonlinearSolver::readjustTrustVector() fabsy = fabs(m_y_n_curr[i]); // First off make sure that each trust region vector is 1/2 the size of each variable or smaller // unless overridden by the deltaStepMininum value. - // doublereal newValue = trustDeltaEach * m_ewt[i] / wtSum; doublereal newValue = trustNormGoal * m_ewt[i]; if (newValue > 0.5 * fabsy) { if (fabsy * 0.5 > m_deltaStepMinimum[i]) { @@ -2598,14 +2579,7 @@ int NonlinearSolver::decideStep(const doublereal time_curr, int leg, doublereal * * If we had to bounds delta the update, decrease the trust region */ - if (m_dampBound < 1.0) { - // trustDelta_ *= 0.5; - // NextTrustFactor_ *= 0.5; - // ll = trustRegionLength(); - // if (m_print_flag >= 5) { - // printf("\t\tdecideStep(): Trust region decreased from %g to %g due to bounds constraint\n", ll*2, ll); - //} - } else { + if (m_dampBound >= 1.0) { retn = 0; /* * Calculate the expected residual from the quadratic model @@ -3086,7 +3060,6 @@ int NonlinearSolver::solve_nonlinear_problem(int SolnType, doublereal* const y_c } if (m_print_flag == 2 || m_print_flag == 3) { - // printf("\t Iter Resid NewJac | Fbound | ResidBound | Fdamp DampIts | DeltaSolnNewton ResidFinal \n"); if (ResidWtsReevaluated_) { printf("\t*"); } else { @@ -3426,7 +3399,6 @@ int NonlinearSolver::beuler_jac(GeneralMatrix& J, doublereal* const f, col_j = (doublereal*) J.ptrColumn(j); ysave = y[j]; dy = dyVector[j]; - //dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7); y[j] = ysave + dy; dy = y[j] - ysave; diff --git a/src/numerics/SquareMatrix.cpp b/src/numerics/SquareMatrix.cpp index a66a92c57..4074b160d 100644 --- a/src/numerics/SquareMatrix.cpp +++ b/src/numerics/SquareMatrix.cpp @@ -246,9 +246,6 @@ doublereal SquareMatrix::rcond(doublereal anorm) throw CELapackError("SquareMatrix::rcond()", "matrix isn't factored correctly"); } - // doublereal anorm = ct_dlange('1', m_nrows, m_nrows, &(*(begin())), m_nrows, DATA_PTR(work)); - - int rinfo = 0; rcond = ct_dgecon('1', m_nrows, &(*(begin())), m_nrows, anorm, DATA_PTR(work), DATA_PTR(iwork_), rinfo); diff --git a/src/numerics/solveProb.cpp b/src/numerics/solveProb.cpp index 65e4b86cf..e52dda964 100644 --- a/src/numerics/solveProb.cpp +++ b/src/numerics/solveProb.cpp @@ -156,14 +156,6 @@ int solveProb::solve(int ifunc, doublereal time_scale, */ std::copy(m_CSolnSP.begin(), m_CSolnSP.end(), m_CSolnSPOld.begin()); - /* - * Evaluate the largest surface species for each surface phase every - * 5 iterations. - */ - // if (iter%5 == 4) { - // evalSurfLarge(DATA_PTR(m_CSolnSP)); - // } - /* * Calculate the value of the time step * - heuristics to stop large oscillations in deltaT @@ -468,10 +460,6 @@ doublereal solveProb::calc_damping(doublereal x[], doublereal dxneg[], size_t di damp = APPROACH * (x[i] - xbot) / dxneg[i]; *label = i; } - // else if (fabs(xnew) > 2.0*MAX(fabs(x[i]), 1.0E-10)) { - // damp = 0.5 * MAX(fabs(x[i]), 1.0E-9)/ fabs(xnew); - // *label = i; - // } double denom = fabs(x[i]) + 1.0E5 * m_atol[i]; if ((fabs(delta_x) / denom) > 0.3) { double newdamp = 0.3 * denom / fabs(delta_x); @@ -772,7 +760,6 @@ void solveProb::printIteration(int ioflag, doublereal damp, size_t label_d, for (int isp = 0; isp < m_numSurfPhases; isp++) { int nsp = m_nSpeciesSurfPhase[isp]; InterfaceKinetics* m_kin = m_objects[isp]; - //int surfPhaseIndex = m_kinObjPhaseIDSurfPhase[isp]; m_kin->getNetProductionRates(DATA_PTR(m_numEqn1)); for (int k = 0; k < nsp; k++, kindexSP++) { int kspIndex = m_kinSpecIndex[kindexSP]; diff --git a/src/oneD/Domain1D.cpp b/src/oneD/Domain1D.cpp index b982ca880..c9258ffb4 100644 --- a/src/oneD/Domain1D.cpp +++ b/src/oneD/Domain1D.cpp @@ -169,7 +169,6 @@ void Domain1D::showSolution(const doublereal* x) { size_t nn = m_nv/5; size_t i, j, n; - //char* buf = new char[100]; char buf[100]; doublereal v; for (i = 0; i < nn; i++) { diff --git a/src/oneD/MultiNewton.cpp b/src/oneD/MultiNewton.cpp index a9b07006b..df771b0b1 100644 --- a/src/oneD/MultiNewton.cpp +++ b/src/oneD/MultiNewton.cpp @@ -184,7 +184,7 @@ void MultiNewton::resize(size_t sz) doublereal MultiNewton::norm2(const doublereal* x, const doublereal* step, OneDim& r) const { - doublereal f, sum = 0.0;//, fmx = 0.0; + doublereal f, sum = 0.0; size_t nd = r.nDomains(); for (size_t n = 0; n < nd; n++) { f = norm_square(x + r.start(n), step + r.start(n), diff --git a/src/oneD/OneDim.cpp b/src/oneD/OneDim.cpp index 9269d2f38..1e7daa61f 100644 --- a/src/oneD/OneDim.cpp +++ b/src/oneD/OneDim.cpp @@ -338,7 +338,6 @@ doublereal OneDim::timeStep(int nsteps, doublereal dt, doublereal* x, if (m == 100) { dt *= 1.5; } - // else dt /= 1.5; if (dt > m_tmax) { dt = m_tmax; } diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index cd4ce55d0..51ab0e3a2 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -22,7 +22,6 @@ static void sim1D_drawline() Sim1D::Sim1D() : OneDim() { - //writelog("Sim1D default constructor\n"); } Sim1D::Sim1D(vector& domains) : @@ -39,7 +38,6 @@ Sim1D::Sim1D(vector& domains) : // set some defaults m_tstep = 1.0e-5; - //m_maxtimestep = 10.0; m_steps.push_back(1); m_steps.push_back(2); m_steps.push_back(5); @@ -387,7 +385,6 @@ int Sim1D::refine(int loglevel) zmid = 0.5*(d.grid(m) + d.grid(m+1)); znew.push_back(zmid); np++; - //writelog(string("refine: adding point at ")+fp2str(zmid)+"\n"); // for each component, linearly interpolate // the solution to this point @@ -411,8 +408,7 @@ int Sim1D::refine(int loglevel) size_t gridstart = 0, gridsize; for (size_t n = 0; n < m_nd; n++) { Domain1D& d = domain(n); - // Refiner& r = d.refiner(); - gridsize = dsize[n]; // d.nPoints() + r.nNewPoints(); + gridsize = dsize[n]; d.setupGrid(gridsize, DATA_PTR(znew) + gridstart); gridstart += gridsize; } @@ -424,8 +420,6 @@ int Sim1D::refine(int loglevel) // resize the work array m_xnew.resize(xnew.size()); - // copy(xnew.begin(), xnew.end(), m_xnew.begin()); - resize(); finalize(); return np; @@ -512,8 +506,7 @@ int Sim1D::setFixedTemperature(doublereal t) size_t gridstart = 0, gridsize; for (n = 0; n < m_nd; n++) { Domain1D& d = domain(n); - // Refiner& r = d.refiner(); - gridsize = dsize[n]; // d.nPoints() + r.nNewPoints(); + gridsize = dsize[n]; d.setupGrid(gridsize, DATA_PTR(znew) + gridstart); gridstart += gridsize; } diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index c4e53fda6..a1671eee1 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -456,7 +456,6 @@ void StFlow::showSolution(const doublereal* x) { size_t nn = m_nv/5; size_t i, j, n; - //char* buf = new char[100]; char buf[100]; // The mean molecular weight is needed to convert @@ -600,7 +599,6 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel) writelog(nd["title"]+": "+nd.value()+"\n"); } - //map params; double pp = -1.0; pp = getFloat(dom, "pressure", "pressure"); setPressure(pp); diff --git a/src/oneD/boundaries1D.cpp b/src/oneD/boundaries1D.cpp index 6f69610dd..7bf78a993 100644 --- a/src/oneD/boundaries1D.cpp +++ b/src/oneD/boundaries1D.cpp @@ -443,9 +443,7 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg, rb[0] = xb[3]; rb[2] = xb[2] - xb[2 + nc]; for (k = 4; k < nc; k++) { - //if (m_flow_right->doSpecies(k-4)) { rb[k] = xb[k] - xb[k + nc]; - //} } } @@ -456,11 +454,8 @@ void Outlet1D::eval(size_t jg, doublereal* xg, doublereal* rg, integer* diagg, db = diag - nc; // zero Lambda - - if (!m_flow_left->fixed_mdot()) { - ; // rb[0] = xb[0] - xb[0-nc]; //zero U gradient - } else { - rb[0] = xb[3]; // zero Lambda + if (m_flow_left->fixed_mdot()) { + rb[0] = xb[3]; } rb[2] = xb[2] - xb[2 - nc]; // zero T gradient @@ -596,7 +591,7 @@ void OutletRes1D::eval(size_t jg, doublereal* xg, doublereal* rg, } else { rb[0] = xb[3]; // zero Lambda } - rb[2] = xb[2] - m_temp; //xb[2] - xb[2 - nc]; // zero dT/dz + rb[2] = xb[2] - m_temp; // zero dT/dz for (k = 5; k < nc; k++) { rb[k] = xb[k] - m_yres[k-4]; // fixed Y db[k] = 0; @@ -695,7 +690,6 @@ void Surf1D::eval(size_t jg, doublereal* xg, doublereal* rg, XML_Node& Surf1D::save(XML_Node& o, const doublereal* const soln) { const doublereal* s = soln + loc(); - //XML_Node& inlt = o.addChild("inlet"); XML_Node& inlt = Domain1D::save(o, soln); inlt.addAttribute("type","surface"); for (size_t k = 0; k < nComponents(); k++) { @@ -768,8 +762,6 @@ void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg, } m_sphase->setTemperature(x[0]); m_sphase->setCoverages(DATA_PTR(m_work)); - //m_kin->advanceCoverages(1.0); - //m_sphase->getCoverages(m_fixed_cov.begin()); // set the left gas state to the adjacent point @@ -789,10 +781,6 @@ void ReactingSurf1D::eval(size_t jg, doublereal* xg, doublereal* rg, m_kin->getNetProductionRates(DATA_PTR(m_work)); doublereal rs0 = 1.0/m_sphase->siteDensity(); - - //scale(m_work.begin(), m_work.end(), m_work.begin(), m_mult[0]); - - // bool enabled = true; size_t ioffset = m_kin->kineticsSpeciesIndex(0, m_surfindex); if (m_enabled) { diff --git a/src/thermo/Elements.cpp b/src/thermo/Elements.cpp index 0af6c23b4..f925e7bc5 100644 --- a/src/thermo/Elements.cpp +++ b/src/thermo/Elements.cpp @@ -545,9 +545,7 @@ void Elements::addElementsFromXML(const XML_Node& phase) for (i = 0; i < nel; i++) { e = 0; if (local_db) { - //writelog("looking in local database."); e = local_db->findByAttr("name",enames[i]); - //if (!e) writelog(enames[i]+" not found."); } if (!e) { e = dbe->findByAttr("name",enames[i]); diff --git a/src/thermo/GeneralSpeciesThermo.cpp b/src/thermo/GeneralSpeciesThermo.cpp index c29931500..11a4e5202 100644 --- a/src/thermo/GeneralSpeciesThermo.cpp +++ b/src/thermo/GeneralSpeciesThermo.cpp @@ -112,7 +112,6 @@ void GeneralSpeciesThermo::install(const std::string& name, m_kk = index+1; } - //int nfreq = 3; /* * Create the necessary object */ @@ -224,9 +223,6 @@ void GeneralSpeciesThermo::update(doublereal t, doublereal* cp_R, if (sp_ptr) { sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R); } - // else { - // writelog("General::update: sp_ptr is NULL!\n"); - //} } } diff --git a/src/thermo/GibbsExcessVPSSTP.cpp b/src/thermo/GibbsExcessVPSSTP.cpp index de5bfb21c..063037563 100644 --- a/src/thermo/GibbsExcessVPSSTP.cpp +++ b/src/thermo/GibbsExcessVPSSTP.cpp @@ -121,9 +121,6 @@ void GibbsExcessVPSSTP::setPressure(doublereal p) void GibbsExcessVPSSTP::calcDensity() { vector_fp vbar = getPartialMolarVolumesVector(); - // double *vbar = &m_pp[0]; - // getPartialMolarVolumes(&vbar[0]); - doublereal vtotal = 0.0; for (size_t i = 0; i < m_kk; i++) { vtotal += vbar[i] * moleFractions_[i]; diff --git a/src/thermo/HMWSoln.cpp b/src/thermo/HMWSoln.cpp index 8ba147c9e..3d7be737c 100644 --- a/src/thermo/HMWSoln.cpp +++ b/src/thermo/HMWSoln.cpp @@ -692,7 +692,6 @@ void HMWSoln::calcDensity() double HMWSoln::density() const { - // calcDensity(); return Phase::density(); } @@ -814,10 +813,6 @@ void HMWSoln::getActivities(doublereal* ac) const double xmolSolvent = moleFraction(m_indexSolvent); ac[m_indexSolvent] = exp(m_lnActCoeffMolal_Scaled[m_indexSolvent]) * xmolSolvent; - /* - * Apply the pH scale - */ - //applyphScale(ac); } void HMWSoln::getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const @@ -1044,7 +1039,6 @@ double HMWSoln::dA_DebyedT_TP(double tempArg, double presArg) const break; case A_DEBYE_WATER: dAdT = m_waterProps->ADebye(T, P, 1); - //dAdT = WaterProps::ADebye(T, P, 1); break; default: printf("shouldn't be here\n"); @@ -1418,11 +1412,6 @@ void HMWSoln::calcMolalitiesCropped() const for (size_t j = (i+1); j < m_kk; j++) { double charge_j = charge(j); double abs_charge_j = fabs(charge_j); - /* - * Find the counterIJ for the symmetric binary interaction - */ - // n = m_kk*i + j; - // counterIJ = m_CounterIJ[n]; /* * Only loop over oppositely charge species */ @@ -1721,18 +1710,6 @@ void HMWSoln::s_updatePitzer_CoeffWRTemp(int doDerivs) const + Theta_coeff[3]*twoinvT3 - Theta_coeff[4]*invT2; } - -#ifdef DEBUG_HKM - /* - * Turn terms off for debugging - */ - //m_Beta0MX_ij_L[counterIJ] = 0; - //m_Beta0MX_ij_LL[counterIJ] = 0; - //m_Beta1MX_ij_L[counterIJ] = 0; - //m_Beta1MX_ij_LL[counterIJ] = 0; - //m_CphiMX_ij_L[counterIJ] = 0; - //m_CphiMX_ij_LL[counterIJ] = 0; -#endif break; } } @@ -2240,10 +2217,6 @@ void HMWSoln::s_updatePitzer_lnMolalityActCoeff() const printf(" Step 7: \n"); } #endif - // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) - // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data - // ( A * ln(10) = 1.17593) - // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; Aphi = A_Debye_TP() / 3.0; F = -Aphi * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); @@ -2907,10 +2880,6 @@ void HMWSoln::s_update_dlnMolalityActCoeff_dT() const */ s_updatePitzer_dlnMolalityActCoeff_dT(); - //double xmolSolvent = moleFraction(m_indexSolvent); - //double xx = MAX(m_xmolSolventMIN, xmolSolvent); - // double lnxs = log(xx); - for (size_t k = 1; k < m_kk; k++) { if (CROP_speciesCropped_[k] == 2) { m_dlnActCoeffMolaldT_Unscaled[k] = 0.0; @@ -3280,19 +3249,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const printf(" Step 7: \n"); } #endif - // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) - // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data - // ( A * ln(10) = 1.17593) - // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; - double dA_DebyedT = dA_DebyedT_TP(); double dAphidT = dA_DebyedT /3.0; -#ifdef DEBUG_HKM - //dAphidT = 0.0; -#endif - //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) - // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); - //dAphidT = Al / (4.0 * GasConstant * T * T); dFdT = -dAphidT * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); #ifdef DEBUG_MODE @@ -3734,7 +3692,6 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dT() const * We have just computed act_0. However, this routine returns * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). */ - //double xmolSolvent = moleFraction(m_indexSolvent); m_dlnActCoeffMolaldT_Unscaled[0] = d_lnwateract_dT; #ifdef DEBUG_MODE if (m_debugCalc) { @@ -3762,10 +3719,6 @@ void HMWSoln::s_update_d2lnMolalityActCoeff_dT2() const */ s_updatePitzer_d2lnMolalityActCoeff_dT2(); - //double xmolSolvent = moleFraction(m_indexSolvent); - //double xx = MAX(m_xmolSolventMIN, xmolSolvent); - //double lnxs = log(xx); - for (size_t k = 1; k < m_kk; k++) { if (CROP_speciesCropped_[k] == 2) { m_d2lnActCoeffMolaldT2_Unscaled[k] = 0.0; @@ -4098,13 +4051,8 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const if (charge(i)*charge(j) > 0) { z1 = (int) fabs(charge(i)); z2 = (int) fabs(charge(j)); - //Phi[counterIJ] = thetaij[counterIJ] + etheta[z1][z2]; - //Phi_L[counterIJ] = thetaij_L[counterIJ]; Phi_LL[counterIJ] = thetaij_LL[counterIJ]; - //Phiprime[counterIJ] = etheta_prime[z1][z2]; Phiprime[counterIJ] = 0.0; - //Phiphi[counterIJ] = Phi[counterIJ] + Is * Phiprime[counterIJ]; - //Phiphi_L[counterIJ] = Phi_L[counterIJ] + Is * Phiprime[counterIJ]; Phiphi_LL[counterIJ] = Phi_LL[counterIJ]; } else { Phi_LL[counterIJ] = 0.0; @@ -4115,9 +4063,6 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const if (m_debugCalc) { sni = speciesName(i); snj = speciesName(j); - //printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", - // sni.c_str(), snj.c_str(), - // Phi_L[counterIJ], Phiprime[counterIJ], Phiphi_L[counterIJ] ); printf(" %-16s %-16s %10.6f %10.6f %10.6f \n", sni.c_str(), snj.c_str(), Phi_LL[counterIJ], Phiprime[counterIJ], Phiphi_LL[counterIJ]); @@ -4134,23 +4079,7 @@ void HMWSoln::s_updatePitzer_d2lnMolalityActCoeff_dT2() const printf(" Step 7: \n"); } #endif - // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) - // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data - // ( A * ln(10) = 1.17593) - // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; - // Aphi = m_A_Debye / 3.0; - - //double dA_DebyedT = dA_DebyedT_TP(); - //double dAphidT = dA_DebyedT /3.0; double d2AphidT2 = d2A_DebyedT2_TP() / 3.0; -#ifdef DEBUG_HKM - //d2AphidT2 = 0.0; -#endif - //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) - // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); - //dAphidT = Al / (4.0 * GasConstant * T * T); - //dFdT = -dAphidT * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) - // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); d2FdT2 = -d2AphidT2 * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); #ifdef DEBUG_MODE @@ -4947,9 +4876,7 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const if (charge(i)*charge(j) > 0) { z1 = (int) fabs(charge(i)); z2 = (int) fabs(charge(j)); - //Phi[counterIJ] = thetaij_L[counterIJ] + etheta[z1][z2]; Phi_P[counterIJ] = thetaij_P[counterIJ]; - //Phiprime[counterIJ] = etheta_prime[z1][z2]; Phiprime[counterIJ] = 0.0; Phiphi_P[counterIJ] = Phi_P[counterIJ] + Is * Phiprime[counterIJ]; } else { @@ -4977,19 +4904,8 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const printf(" Step 7: \n"); } #endif - // A_Debye_Huckel = 0.5092; (units = sqrt(kg/gmol)) - // A_Debye_Huckel = 0.5107; <- This value is used to match GWB data - // ( A * ln(10) = 1.17593) - // Aphi = A_Debye_Huckel * 2.30258509 / 3.0; - double dA_DebyedP = dA_DebyedP_TP(currTemp, currPres); double dAphidP = dA_DebyedP /3.0; -#ifdef DEBUG_MODE - //dAphidT = 0.0; -#endif - //F = -Aphi * ( sqrt(Is) / (1.0 + 1.2*sqrt(Is)) - // + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); - //dAphidT = Al / (4.0 * GasConstant * T * T); dFdP = -dAphidP * (sqrt(Is) / (1.0 + 1.2*sqrt(Is)) + (2.0/1.2) * log(1.0+1.2*(sqrtIs))); #ifdef DEBUG_MODE @@ -5436,7 +5352,6 @@ void HMWSoln::s_updatePitzer_dlnMolalityActCoeff_dP() const * We have just computed act_0. However, this routine returns * ln(actcoeff[]). Therefore, we must calculate ln(actcoeff_0). */ - //double xmolSolvent = moleFraction(m_indexSolvent); m_dlnActCoeffMolaldP_Unscaled[0] = d_lnwateract_dP; #ifdef DEBUG_MODE if (m_debugCalc) { diff --git a/src/thermo/HMWSoln_input.cpp b/src/thermo/HMWSoln_input.cpp index be9c446e9..e43eb5a97 100644 --- a/src/thermo/HMWSoln_input.cpp +++ b/src/thermo/HMWSoln_input.cpp @@ -1337,9 +1337,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) m_speciesSize[k] << endl; #endif } else { - // throw CanteraError("HMWSoln::initThermoXML", - // "Solvent SS Model \"" + modelStringa + - // "\" is not allowed, name = " + sss[0]); m_waterSS = providePDSS(0); m_waterSS->setState_TP(300., OneAtm); double dens = m_waterSS->density(); @@ -1467,10 +1464,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) if (jmap != npos) { const XML_Node& sp = *xspecies[jmap]; getOptionalFloat(sp, "stoichIsMods", m_speciesCharge_Stoich[k]); - // if (sp.hasChild("stoichIsMods")) { - // double val = getFloat(sp, "stoichIsMods"); - //m_speciesCharge_Stoich[k] = val; - //} } } @@ -1672,19 +1665,6 @@ void HMWSoln::initThermoXML(XML_Node& phaseNode, const std::string& id_) notDone = false; } } while (notDone); - - - - - - - - - // if (phaseNode.hasChild("state")) { - // XML_Node& stateNode = phaseNode.child("state"); - // setStateFromXML(stateNode); - //} - } void HMWSoln::calcIMSCutoffParams_() diff --git a/src/thermo/IdealGasPhase.cpp b/src/thermo/IdealGasPhase.cpp index 698529724..b849f9627 100644 --- a/src/thermo/IdealGasPhase.cpp +++ b/src/thermo/IdealGasPhase.cpp @@ -168,10 +168,8 @@ void IdealGasPhase::getStandardChemPotentials(doublereal* muStar) const void IdealGasPhase::getChemPotentials(doublereal* mu) const { getStandardChemPotentials(mu); - //doublereal logp = log(pressure()/m_spthermo->refPressure()); doublereal xx; doublereal rt = temperature() * GasConstant; - //const vector_fp& g_RT = gibbs_RT_ref(); for (size_t k = 0; k < m_kk; k++) { xx = std::max(SmallNumber, moleFraction(k)); mu[k] += rt * (log(xx)); diff --git a/src/thermo/IdealMolalSoln.cpp b/src/thermo/IdealMolalSoln.cpp index 7b4acb85c..986a84629 100644 --- a/src/thermo/IdealMolalSoln.cpp +++ b/src/thermo/IdealMolalSoln.cpp @@ -226,7 +226,6 @@ void IdealMolalSoln::setState_TP(doublereal temp, doublereal pres) Phase::setTemperature(temp); m_Pcurrent = pres; updateStandardStateThermo(); - //m_densWaterSS = m_waterSS->density(); calcDensity(); } diff --git a/src/thermo/IonsFromNeutralVPSSTP.cpp b/src/thermo/IonsFromNeutralVPSSTP.cpp index 5a7e4c052..1938d7dc3 100644 --- a/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -68,10 +68,6 @@ IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(const std::string& inputFile, } constructPhaseFile(inputFile, id_); geThermo = dynamic_cast(neutralMoleculePhase_); - //y.resize(numNeutralMoleculeSpecies_,0.0); - //size_t numNeutMolSpec = geThermo->nSpecies(); - //dlnActCoeff_NeutralMolecule.resize(numNeutMolSpec); - //dX_NeutralMolecule.resize(numNeutMolSpec); } IonsFromNeutralVPSSTP::IonsFromNeutralVPSSTP(XML_Node& phaseRoot, @@ -213,7 +209,6 @@ void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::strin * The phase object automatically constructs an XML object. * Use this object to store information. */ - //XML_Node& phaseNode_XML = xml(); XML_Node* fxml = new XML_Node(); fxml->build(fin); XML_Node* fxml_phase = findXMLPhase(fxml, id_); @@ -223,7 +218,6 @@ void IonsFromNeutralVPSSTP::constructPhaseFile(std::string inputFile, std::strin id_ + " in file named " + inputFile); } setXMLdata(*fxml_phase); - //fxml_phase->copy(&phaseNode_XML); constructPhaseXML(*fxml_phase, id_); delete fxml; } @@ -338,8 +332,6 @@ doublereal IonsFromNeutralVPSSTP::cv_mole() const // Need to revisit this, as it is wrong getPartialMolarCp(DATA_PTR(m_pp)); return mean_X(DATA_PTR(m_pp)); - //err("not implemented"); - //return 0.0; } /* @@ -356,12 +348,6 @@ void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs, void IonsFromNeutralVPSSTP::getActivityCoefficients(doublereal* ac) const { - - // This stuff has moved to the setState routines - // calcNeutralMoleculeMoleFractions(); - // neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_)); - // neutralMoleculePhase_->getStandardChemPotentials(DATA_PTR(muNeutralMolecule_)); - /* * Update the activity coefficients */ @@ -383,13 +369,6 @@ void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const { size_t icat, jNeut; doublereal xx, fact2; - /* - * Transfer the mole fractions to the slave neutral molecule - * phase - * Note we may move this in the future. - */ - //calcNeutralMoleculeMoleFractions(); - //neutralMoleculePhase_->setState_TPX(temperature(), pressure(), DATA_PTR(NeutralMolecMoleFractions_)); /* * Get the standard chemical potentials of netural molecules @@ -403,7 +382,6 @@ void IonsFromNeutralVPSSTP::getChemPotentials(doublereal* mu) const neutralMoleculePhase_->getChemPotentials(mu); break; case cIonSolnType_SINGLEANION: - // neutralMoleculePhase_->getActivityCoefficients(DATA_PTR(gammaNeutralMolecule_)); neutralMoleculePhase_->getLnActivityCoefficients(DATA_PTR(lnActCoeff_NeutralMolecule_)); fact2 = 2.0 * RT_ * log(2.0); @@ -553,8 +531,6 @@ void IonsFromNeutralVPSSTP::setState_TP(doublereal t, doublereal p) /* * Calculate the partial molar volumes, and then the density of the fluid */ - - //calcDensity(); double dd = neutralMoleculePhase_->density(); Phase::setDensity(dd); } @@ -617,8 +593,6 @@ void IonsFromNeutralVPSSTP::calcNeutralMoleculeMoleFractions() const } #endif - // bool fmSimple = true; - switch (ionSolnType_) { case cIonSolnType_PASSTHROUGH: @@ -718,9 +692,6 @@ void IonsFromNeutralVPSSTP::getNeutralMoleculeMoleGrads(const doublereal* const dy[k] = 0.0; } - - // bool fmSimple = true; - switch (ionSolnType_) { case cIonSolnType_PASSTHROUGH: @@ -1031,9 +1002,6 @@ void IonsFromNeutralVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string for (size_t k = 0; k < m_kk; k++) { fm_invert_ionForNeutral[k] = npos; } - /* for (int jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { - fm_invert_ionForNeutral[jNeut] = -1; - }*/ for (size_t jNeut = 0; jNeut < numNeutralMoleculeSpecies_; jNeut++) { for (size_t m = 0; m < nElementsN; m++) { elemVectorN[m] = neutralMoleculePhase_->nAtoms(jNeut, m); diff --git a/src/thermo/LatticePhase.cpp b/src/thermo/LatticePhase.cpp index 597d58bd6..981b3707a 100644 --- a/src/thermo/LatticePhase.cpp +++ b/src/thermo/LatticePhase.cpp @@ -96,17 +96,6 @@ doublereal LatticePhase::calcDensity() setMolarDensity(m_site_density); doublereal mw = meanMolecularWeight(); doublereal dens = mw * m_site_density; - /* - * Calculate the molarVolume of the solution (m**3 kmol-1) - */ - // const doublereal * const dtmp = moleFractdivMMW(); - // doublereal invDens = dot(m_speciesMolarVolume.begin(), m_speciesMolarVolume.end(), dtmp); - /* - * Set the density in the parent State object directly, - * by calling the Phase::setDensity() function. - */ - // doublereal dens = 1.0/invDens; - // Phase::setDensity(dens); return dens; } diff --git a/src/thermo/LatticeSolidPhase.cpp b/src/thermo/LatticeSolidPhase.cpp index 593b515d4..be71d8e32 100644 --- a/src/thermo/LatticeSolidPhase.cpp +++ b/src/thermo/LatticeSolidPhase.cpp @@ -441,10 +441,6 @@ void LatticeSolidPhase::initLengths() void LatticeSolidPhase::_updateThermo() const { doublereal tnow = temperature(); - // if (fabs(molarDensity() - m_molar_density)/m_molar_density > 0.0001) { - // throw CanteraError("_updateThermo","molar density changed from " - // +fp2str(m_molar_density)+" to "+fp2str(molarDensity())); - //} if (m_tlast != tnow) { getMoleFractions(DATA_PTR(m_x)); size_t strt = 0; diff --git a/src/thermo/MargulesVPSSTP.cpp b/src/thermo/MargulesVPSSTP.cpp index 8ac208855..755815f92 100644 --- a/src/thermo/MargulesVPSSTP.cpp +++ b/src/thermo/MargulesVPSSTP.cpp @@ -342,7 +342,6 @@ void MargulesVPSSTP::getPartialMolarVolumes(doublereal* vbar) const for (size_t iK = 0; iK < m_kk; iK++) { vbar[iK] += all; - // vbar[iK] += XA*XB*temp1+((delAK-XA)*XB+XA*(delBK-XB))*temp1+XB*XA*(delBK-XB)*g1; } vbar[iA] += XB * temp1; vbar[iB] += XA * temp1 + XA*XB*g1; @@ -483,7 +482,6 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dT() const const doublereal mult = 2.0 * invT; const doublereal dT2all = mult * all; for (iK = 0; iK < m_kk; iK++) { - // double temp = (delAK * XB + XA * delBK - XA * XB) * (g0 + g1 * XB) + XA * XB * (delBK - XB) * g1; dlnActCoeffdT_Scaled_[iK] += all; d2lnActCoeffdT2_Scaled_[iK] -= dT2all; } @@ -517,8 +515,6 @@ void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* double XA, XB, g0 , g1, dXA, dXB; double T = temperature(); double RT = GasConstant*T; - - //fvo_zero_dbl_1(dlnActCoeff, m_kk); s_update_dlnActCoeff_dT(); for (iK = 0; iK < m_kk; iK++) { @@ -538,8 +534,6 @@ void MargulesVPSSTP::getdlnActCoeffds(const doublereal dTds, const doublereal* const doublereal g2XAdXB = 2*g1*XA*dXB; const doublereal all = (-XB * dXA - XA *dXB) * g02g1XB - XB *g2XAdXB; for (iK = 0; iK < m_kk; iK++) { - // dlnActCoeffds[iK] += ((delBK-XB)*dXA + (delAK-XA)*dXB)*(g0+2*g1*XB) + (delBK-XB)*2*g1*XA*dXB - // + dlnActCoeffdT_Scaled_[iK]*dTds; dlnActCoeffds[iK] += all + dlnActCoeffdT_Scaled_[iK]*dTds; } dlnActCoeffds[iA] += dXB * g02g1XB; @@ -581,24 +575,8 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN_diag() const g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB))); - - // double gfac = g0 + g1 * XB; - // double gggg = (delBK - XB) * g1; - - - // dlnActCoeffdlnN_diag_[iK] += gfac * delAK * ( - XB + delBK); - - // dlnActCoeffdlnN_diag_[iK] += gfac * delBK * ( - XA + delAK); - - // dlnActCoeffdlnN_diag_[iK] += gfac * (2.0 * XA * XB - delAK * XB - XA * delBK); - - // dlnActCoeffdlnN_diag_[iK] += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBK); - - // dlnActCoeffdlnN_diag_[iK] += gggg * ( - 2.0 * XA * XB + delAK * XB + XA * delBK); - - // dlnActCoeffdlnN_diag_[iK] += - g1 * XA * XB * (- XB + delBK); } - dlnActCoeffdlnN_diag_[iK] = XK*dlnActCoeffdlnN_diag_[iK];//-XK; + dlnActCoeffdlnN_diag_[iK] = XK*dlnActCoeffdlnN_diag_[iK]; } } @@ -648,22 +626,6 @@ void MargulesVPSSTP::s_update_dlnActCoeff_dlnN() const dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB)); dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA); - - // double gfac = g0 + g1 * XB; - // double gggg = (delBK - XB) * g1; - - - // dlnActCoeffdlnN_(iK, iM) += gfac * delAK * ( - XB + delBM); - - // dlnActCoeffdlnN_(iK, iM) += gfac * delBK * ( - XA + delAM); - - // dlnActCoeffdlnN_(iK, iM) += gfac * (2.0 * XA * XB - delAM * XB - XA * delBM); - - // dlnActCoeffdlnN_(iK, iM) += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBM); - - // dlnActCoeffdlnN_(iK, iM) += gggg * ( - 2.0 * XA * XB + delAM * XB + XA * delBM); - - // dlnActCoeffdlnN_(iK, iM) += - g1 * XA * XB * (- XB + delBM); } dlnActCoeffdlnN_(iK,iM) = XM*dlnActCoeffdlnN_(iK,iM); } diff --git a/src/thermo/MineralEQ3.cpp b/src/thermo/MineralEQ3.cpp index f12799dcc..ea9a4af86 100644 --- a/src/thermo/MineralEQ3.cpp +++ b/src/thermo/MineralEQ3.cpp @@ -345,12 +345,6 @@ void MineralEQ3::convertDGFormation() totalSum += na * ge; } } - // Add in the charge - // if (m_charge_j != 0.0) { - // ename = "H"; - // ge = LookupGe(ename); - // totalSum -= m_charge_j * ge; - //} // Ok, now do the calculation. Convert to joules kmol-1 doublereal dg = m_deltaG_formation_pr_tr * 4.184 * 1.0E3; //! Store the result into an internal variable. diff --git a/src/thermo/MixedSolventElectrolyte.cpp b/src/thermo/MixedSolventElectrolyte.cpp index 09020ebb4..5020b1481 100644 --- a/src/thermo/MixedSolventElectrolyte.cpp +++ b/src/thermo/MixedSolventElectrolyte.cpp @@ -516,9 +516,7 @@ void MixedSolventElectrolyte::getdlnActCoeffds(const doublereal dTds, const dou double XA, XB, g0, g1, dXA, dXB; double T = temperature(); double RT = GasConstant*T; - - //fvo_zero_dbl_1(dlnActCoeff, m_kk); - s_update_dlnActCoeff_dT(); + s_update_dlnActCoeff_dT(); for (size_t iK = 0; iK < m_kk; iK++) { dlnActCoeffds[iK] = 0.0; @@ -585,22 +583,6 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN_diag() const g1 = (m_HE_c_ij[i] - T * m_SE_c_ij[i]) / RT; dlnActCoeffdlnN_diag_[iK] += 2*(delBK-XB)*(g0*(delAK-XA)+g1*(2*(delAK-XA)*XB+XA*(delBK-XB))); - - // double gfac = g0 + g1 * XB; - // double gggg = (delBK - XB) * g1; - - - // dlnActCoeffdlnN_diag_[iK] += gfac * delAK * ( - XB + delBK); - - // dlnActCoeffdlnN_diag_[iK] += gfac * delBK * ( - XA + delAK); - - // dlnActCoeffdlnN_diag_[iK] += gfac * (2.0 * XA * XB - delAK * XB - XA * delBK); - - // dlnActCoeffdlnN_diag_[iK] += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBK); - - // dlnActCoeffdlnN_diag_[iK] += gggg * ( - 2.0 * XA * XB + delAK * XB + XA * delBK); - - // dlnActCoeffdlnN_diag_[iK] += - g1 * XA * XB * (- XB + delBK); } dlnActCoeffdlnN_diag_[iK] = XK*dlnActCoeffdlnN_diag_[iK];//-XK; } @@ -650,22 +632,6 @@ void MixedSolventElectrolyte::s_update_dlnActCoeff_dlnN() const dlnActCoeffdlnN_(iK,iM) += g0*((delAM-XA)*(delBK-XB)+(delAK-XA)*(delBM-XB)); dlnActCoeffdlnN_(iK,iM) += 2*g1*((delAM-XA)*(delBK-XB)*XB+(delAK-XA)*(delBM-XB)*XB+(delBM-XB)*(delBK-XB)*XA); - - // double gfac = g0 + g1 * XB; - // double gggg = (delBK - XB) * g1; - - - // dlnActCoeffdlnN_(iK, iM) += gfac * delAK * ( - XB + delBM); - - // dlnActCoeffdlnN_(iK, iM) += gfac * delBK * ( - XA + delAM); - - // dlnActCoeffdlnN_(iK, iM) += gfac * (2.0 * XA * XB - delAM * XB - XA * delBM); - - // dlnActCoeffdlnN_(iK, iM) += (delAK * XB + XA * delBK - XA * XB) * g1 * (-XB + delBM); - - // dlnActCoeffdlnN_(iK, iM) += gggg * ( - 2.0 * XA * XB + delAM * XB + XA * delBM); - - // dlnActCoeffdlnN_(iK, iM) += - g1 * XA * XB * (- XB + delBM); } dlnActCoeffdlnN_(iK,iM) = XM*dlnActCoeffdlnN_(iK,iM); } diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp index a744ccc15..3ef243f16 100644 --- a/src/thermo/MixtureFugacityTP.cpp +++ b/src/thermo/MixtureFugacityTP.cpp @@ -311,13 +311,7 @@ void MixtureFugacityTP::setTemperature(const doublereal temp) void MixtureFugacityTP::setPressure(doublereal p) { setState_TP(temperature(), p); - // double chemPot[5], mf[5]; - // getMoleFractions(mf); - // getChemPotentials(chemPot); - // for (int i = 0; i < m_kk; i++) { - // printf(" MixFug:setPres: mu(%d = %g) = %18.8g\n", i, mf[i], chemPot[i]); - // } -} + } void MixtureFugacityTP::setMassFractions(const doublereal* const y) { @@ -380,9 +374,7 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) _updateReferenceStateThermo(); // Depends on the mole fractions and the temperature updateMixingExpressions(); - // setPressure(pres); m_Pcurrent = pres; - // double mmw = meanMolecularWeight(); if (forcedState_ == FLUID_UNDEFINED) { double rhoNow = Phase::density(); @@ -444,12 +436,6 @@ void MixtureFugacityTP::setState_TP(doublereal t, doublereal pres) } } - - - - //setTemperature(t); - //setPressure(pres); - //calcDensity(); } void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho) @@ -464,8 +450,6 @@ void MixtureFugacityTP::setState_TR(doublereal T, doublereal rho) m_Pcurrent = pressureCalc(T, mv); iState_ = phaseState(true); - - // printf("setState_TR: state at T = %g, rho = %g, mv = %g, P = %20.13g, iState = %d\n", T, rho, mv, m_Pcurrent, iState_); } void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublereal* x) @@ -477,10 +461,6 @@ void MixtureFugacityTP::setState_TPX(doublereal t, doublereal p, const doublerea void MixtureFugacityTP::initThermoXML(XML_Node& phaseNode, const std::string& id_) { MixtureFugacityTP::initLengths(); - - //m_VPSS_ptr->initThermo(); - - // m_VPSS_ptr->initThermoXML(phaseNode, id); ThermoPhase::initThermoXML(phaseNode, id_); } @@ -529,8 +509,6 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa, { double tcrit = critTemperature(); doublereal mmw = meanMolecularWeight(); - // double pcrit = critPressure(); - // doublereal deltaGuess = 0.0; if (rhoguess == -1.0) { if (phase != -1) { if (TKelvin > tcrit) { @@ -639,7 +617,6 @@ doublereal MixtureFugacityTP::densityCalc(doublereal TKelvin, doublereal presPa, if (n < 10) { dpdV = dpdVBase * 1.5; } - // if (dpdV > -0.001) dpdV = -0.001; /* * Formulate the update to the molar volume by @@ -734,9 +711,6 @@ int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& de int retn = 0; doublereal densLiq = densityCalc(TKelvin, pres, FLUID_LIQUID_0, densLiqGuess); if (densLiq <= 0.0) { - // throw Cantera::CanteraError("MixtureFugacityTP::corr0", - // "Error occurred trying to find liquid density at (T,P) = " - // + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres)); retn = -1; } else { densLiqGuess = densLiq; @@ -746,9 +720,6 @@ int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& de doublereal densGas = densityCalc(TKelvin, pres, FLUID_GAS, densGasGuess); if (densGas <= 0.0) { - //throw Cantera::CanteraError("MixtureFugacityTP::corr0", - // "Error occurred trying to find gas density at (T,P) = " - // + Cantera::fp2str(TKelvin) + " " + Cantera::fp2str(pres)); if (retn == -1) { throw Cantera::CanteraError("MixtureFugacityTP::corr0", "Error occurred trying to find gas density at (T,P) = " @@ -760,7 +731,6 @@ int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& de setState_TR(TKelvin, densGas); gasGRT = gibbs_mole() / _RT(); } - // delGRT = gibbsLiqRT - gibbsGasRT; return retn; } diff --git a/src/thermo/Nasa9Poly1.cpp b/src/thermo/Nasa9Poly1.cpp index 3d7c593ac..bad612aad 100644 --- a/src/thermo/Nasa9Poly1.cpp +++ b/src/thermo/Nasa9Poly1.cpp @@ -80,8 +80,6 @@ void Nasa9Poly1::updateProperties(const doublereal* tt, cp_R[m_index] = cpdivR; h_RT[m_index] = hdivRT; s_R[m_index] = sdivR; - //writelog("NASA9poly1: for species "+int2str(m_index)+", h_RT = "+ - // fp2str(h)+"\n"); } void Nasa9Poly1::updatePropertiesTemp(const doublereal temp, diff --git a/src/thermo/PDSS_ConstVol.cpp b/src/thermo/PDSS_ConstVol.cpp index 4ea147d67..a94c5b35d 100644 --- a/src/thermo/PDSS_ConstVol.cpp +++ b/src/thermo/PDSS_ConstVol.cpp @@ -97,7 +97,6 @@ void PDSS_ConstVol::constructPDSSXML(VPStandardStateTP* tp, size_t spindex, m_constMolarVolume = ctml::getFloat(*ss, "molarVolume", "toSI"); std::string id = ""; - // initThermoXML(phaseNode, id); } void PDSS_ConstVol::constructPDSSFile(VPStandardStateTP* tp, size_t spindex, diff --git a/src/thermo/PDSS_HKFT.cpp b/src/thermo/PDSS_HKFT.cpp index 5a77628f8..c3c1018b3 100644 --- a/src/thermo/PDSS_HKFT.cpp +++ b/src/thermo/PDSS_HKFT.cpp @@ -608,24 +608,18 @@ void PDSS_HKFT::constructPDSSXML(VPStandardStateTP* tp, size_t spindex, doublereal val = getFloat(*hh, "DG0_f_Pr_Tr"); m_deltaG_formation_tr_pr = val; hasDGO = 1; - } else { - // throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DG0_f_Pr_Tr field"); } if (hh->hasChild("DH0_f_Pr_Tr")) { doublereal val = getFloat(*hh, "DH0_f_Pr_Tr"); m_deltaH_formation_tr_pr = val; hasDHO = 1; - } else { - // throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing DH0_f_Pr_Tr field"); } if (hh->hasChild("S0_Pr_Tr")) { doublereal val = getFloat(*hh, "S0_Pr_Tr"); m_Entrop_tr_pr= val; hasSO = 1; - } else { - // throw CanteraError("PDSS_HKFT::constructPDSSXML", " missing S0_Pr_Tr field"); } const XML_Node* ss = speciesNode.findByName("standardState"); @@ -822,8 +816,6 @@ doublereal PDSS_HKFT::deltaH() const doublereal PDSS_HKFT::deltaG() const { doublereal pbar = m_pres * 1.0E-5; - //doublereal m_presR_bar = OneAtm * 1.0E-5; - doublereal sterm = - m_Entrop_tr_pr * (m_temp - 298.15); doublereal c1term = -m_c1 * (m_temp * log(m_temp/298.15) - (m_temp - 298.15)); @@ -1074,7 +1066,6 @@ doublereal PDSS_HKFT::gstar(const doublereal temp, const doublereal pres, const if (fabs(fvalT - fval) > 1.0E-9) { printf("we are here\n"); } - // return gvalT - fvalT; } #endif return res; diff --git a/src/thermo/PDSS_IdealGas.cpp b/src/thermo/PDSS_IdealGas.cpp index 8ce14167e..f602aba4b 100644 --- a/src/thermo/PDSS_IdealGas.cpp +++ b/src/thermo/PDSS_IdealGas.cpp @@ -74,8 +74,6 @@ PDSS* PDSS_IdealGas::duplMyselfAsPDSS() const void PDSS_IdealGas::constructPDSSXML(VPStandardStateTP* tp, size_t spindex, const XML_Node& phaseNode, const std::string& id) { - //initThermo(); - //initThermoXML(phaseNode, id); } void PDSS_IdealGas::constructPDSSFile(VPStandardStateTP* tp, size_t spindex, diff --git a/src/thermo/PDSS_Water.cpp b/src/thermo/PDSS_Water.cpp index 547c884ec..1e68fec23 100644 --- a/src/thermo/PDSS_Water.cpp +++ b/src/thermo/PDSS_Water.cpp @@ -202,14 +202,12 @@ void PDSS_Water::constructSet() } s = entropy_mole(); s -= GasConstant * log(oneBar/presLow); - //printf("s = %g\n", s); doublereal h = enthalpy_mole(); if (h != -241.826E6) { EW_Offset = -241.826E6 - h; } h = enthalpy_mole(); - //printf("h = %g\n", h); /* * Set the initial state of the system to 298.15 K and @@ -322,10 +320,6 @@ void PDSS_Water::setPressure(doublereal p) waterState = WATER_SUPERCRIT; } -#ifdef DEBUG_HKM - //printf("waterPDSS: set pres = %g t = %g, waterState = %d\n", - // p, T, waterState); -#endif doublereal dd = m_sub.density(T, p, waterState, dens); if (dd <= 0.0) { std::string stateString = "T = " + diff --git a/src/thermo/Phase.cpp b/src/thermo/Phase.cpp index 88f095c84..6dd777c61 100644 --- a/src/thermo/Phase.cpp +++ b/src/thermo/Phase.cpp @@ -349,7 +349,6 @@ void Phase::setMoleFractions(const doublereal* const x) * Set m_ym_ to the normalized mole fractions divided by the normalized mean molecular weight: * m_ym_k = X_k / (sum_k X_k M_k) */ - // transform(m_y.begin(), m_y.end(), m_ym.begin(), timesConstant(1.0/sum)); const doublereal invSum = 1.0/sum; for (size_t k=0; k < m_kk; k++) { m_ym[k] = m_y[k]*invSum; @@ -358,7 +357,6 @@ void Phase::setMoleFractions(const doublereal* const x) * Now set m_y to the normalized mass fractions * m_y = X_k M_k / (sum_k X_k M_k) */ - // transform(m_ym.begin(), m_ym.begin() + m_kk, m_molwts.begin(), m_y.begin(), multiplies()); for (size_t k=0; k < m_kk; k++) { m_y[k] = m_ym[k] * m_molwts[k]; } @@ -799,9 +797,7 @@ void Phase::addElementsFromXML(const XML_Node& phase) for (size_t i = 0; i < enames.size(); i++) { XML_Node* e = 0; if (local_db) { - //writelog("looking in local database."); e = local_db->findByAttr("name",enames[i]); - //if (!e) writelog(enames[i]+" not found."); } if (!e) { e = dbe->findByAttr("name",enames[i]); diff --git a/src/thermo/PhaseCombo_Interaction.cpp b/src/thermo/PhaseCombo_Interaction.cpp index 5fd1bbb73..6831897a2 100644 --- a/src/thermo/PhaseCombo_Interaction.cpp +++ b/src/thermo/PhaseCombo_Interaction.cpp @@ -202,9 +202,7 @@ void PhaseCombo_Interaction::getChemPotentials(doublereal* mu) const * Update the activity coefficients */ s_update_lnActCoeff(); - /* - * - */ + doublereal RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { xx = std::max(moleFractions_[k], SmallNumber); @@ -541,8 +539,6 @@ void PhaseCombo_Interaction::getdlnActCoeffds(const doublereal dTds, const doub doublereal T = temperature(); doublereal RT = GasConstant*T; doublereal xx; - - //fvo_zero_dbl_1(dlnActCoeff, m_kk); s_update_dlnActCoeff_dT(); for (size_t iK = 0; iK < m_kk; iK++) { @@ -605,7 +601,6 @@ void PhaseCombo_Interaction::s_update_dlnActCoeff_dlnN_diag() const /* * First wipe out the ideal solution mixing term */ - // lnActCoeff_Scaled_[iK] = - log(xx); if (xx > SmallNumber) { dlnActCoeffdlnN_diag_[iK] = - 1.0 + xx; } diff --git a/src/thermo/PureFluidPhase.cpp b/src/thermo/PureFluidPhase.cpp index dbecc0d84..2464d2723 100644 --- a/src/thermo/PureFluidPhase.cpp +++ b/src/thermo/PureFluidPhase.cpp @@ -248,7 +248,6 @@ void PureFluidPhase::getEnthalpy_RT_ref(doublereal* hrt) const { double psave = pressure(); double t = temperature(); - //double pref = m_spthermo->refPressure(); double plow = 1.0E-8; Set(tpx::PropertyPair::TP, t, plow); getEnthalpy_RT(hrt); diff --git a/src/thermo/RedlichKisterVPSSTP.cpp b/src/thermo/RedlichKisterVPSSTP.cpp index f5edcee58..a0192af16 100644 --- a/src/thermo/RedlichKisterVPSSTP.cpp +++ b/src/thermo/RedlichKisterVPSSTP.cpp @@ -203,9 +203,7 @@ void RedlichKisterVPSSTP::getChemPotentials(doublereal* mu) const * Update the activity coefficients */ s_update_lnActCoeff(); - /* - * - */ + doublereal RT = GasConstant * temperature(); for (size_t k = 0; k < m_kk; k++) { xx = std::max(moleFractions_[k], SmallNumber); @@ -480,9 +478,6 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const polyk *= fac; } // This gives the same result as above - // printf("RT lnActCoeff_Scaled_[iA] = %15.8E , lnA = %15.8E\n", lnActCoeff_Scaled_[iA], lnA); - // printf("RT lnActCoeff_Scaled_[iB] = %15.8E , lnB = %15.8E\n", lnActCoeff_Scaled_[iB], lnB); - #endif } @@ -492,8 +487,6 @@ void RedlichKisterVPSSTP::s_update_lnActCoeff() const void RedlichKisterVPSSTP::s_update_dlnActCoeff_dT() const { doublereal XA, XB; - // doublereal T = temperature(); - dlnActCoeffdT_Scaled_.assign(m_kk, 0.0); d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0); diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp index b5f194dba..abfb3206f 100644 --- a/src/thermo/RedlichKwongMFTP.cpp +++ b/src/thermo/RedlichKwongMFTP.cpp @@ -599,9 +599,6 @@ void RedlichKwongMFTP::getPartialMolarCp(doublereal* cpbar) const void RedlichKwongMFTP::getPartialMolarVolumes(doublereal* vbar) const { - // getStandardVolumes(vbar); - - for (size_t k = 0; k < m_kk; k++) { m_pp[k] = 0.0; for (size_t i = 0; i < m_kk; i++) { @@ -1041,9 +1038,6 @@ doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGu int nsol = NicholsSolve(TKelvin, pres, atmp, btmp, Vroot); - // printf("nsol = %d\n", nsol); - // printf("liquidVolEst start: T = %g , p = %g, a = %g, b = %g\n", TKelvin, pres, m_a_current, m_b_current); - if (nsol == 1 || nsol == 2) { double pc = critPressure(); if (pres > pc) { @@ -1062,7 +1056,6 @@ doublereal RedlichKwongMFTP::liquidVolEst(doublereal TKelvin, doublereal& presGu } else { v = -1.0; } - //printf (" RedlichKwongMFTP::liquidVolEst %g %g converged in %d its\n", TKelvin, pres, i); return v; } @@ -1117,8 +1110,6 @@ doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, if (phaseRequested == FLUID_GAS || phaseRequested == FLUID_SUPERCRIT || phaseRequested == FLUID_UNDEFINED) { molarVolLast = Vroot_[0]; } else { - //molarVolLast = Vroot_[0]; - //printf("DensityCalc(): Possible problem encountered\n"); return -2.0; } } else if (NSolns_ == -1) { @@ -1127,13 +1118,10 @@ doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa, } else if (TKelvin > tcrit) { molarVolLast = Vroot_[0]; } else { - // molarVolLast = Vroot_[0]; - // printf("DensityCalc(): Possible problem encountered\n"); return -2.0; } } else { molarVolLast = Vroot_[0]; - //printf("DensityCalc(): Possible problem encountered\n"); return -1.0; } return mmw / molarVolLast; @@ -1412,13 +1400,6 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do } else if (desc == 0.0) { nSolnValues = 2; // We are here as p goes to zero. - // double hleft = 3.0 * an * cn / (bn * bn); - //double ynleft = 9.0 * an * cn / (2.0 * bn * bn) - 27.0 * an * an * dn / (2.0 * bn * bn * bn); - //printf("hleft = %g , ynleft = %g\n", -3. / 2. * hleft, -ynleft); - //double h2left = - 3 * hleft + 3 * hleft * hleft - hleft * hleft * hleft; - //double y2left = - 2.0 * ynleft + ynleft * ynleft; - //printf("h2left = %g , yn2left = %g\n", h2left, y2left); - } else if (desc > 0.0) { nSolnValues = 1; } @@ -1562,8 +1543,6 @@ int RedlichKwongMFTP::NicholsSolve(double TKelvin, double pres, doublereal a, do } } } - // writelog("RedlichKwongMFTP::NicholsSolve(T = " + fp2str(TKelvin) + ", p = " + fp2str(pres) + "): finished"); - // writelogendl(); return nSolnValues; } diff --git a/src/thermo/SemiconductorPhase.cpp b/src/thermo/SemiconductorPhase.cpp index 458669fb1..88d6c1a80 100644 --- a/src/thermo/SemiconductorPhase.cpp +++ b/src/thermo/SemiconductorPhase.cpp @@ -18,20 +18,6 @@ static doublereal JoyceDixon(doublereal r) SemiconductorPhase::SemiconductorPhase(std::string infile, std::string id_) {} - -// doublereal SemiconductorPhase::ionizedDonorConcentration() { -// return 1.0/(1.0 + 2.0*exp( fermiLevel() - m_edonor)); -//} - -//doublereal SemiconductorPhase::ionizedAcceptorConcentration() { -// return 1.0/(1.0 + 2.0*exp( m_eacceptor - fermiLevel())); -//} - -//doublereal SemiconductorPhase::_dn(doublereal efermi) { -// m_fermi_level = efermi; -// return electronConcentration() - holeConcentration() + -// ionizedAcceptorConcentration() - ionizedDonorConcentration(); -//} void SemiconductorPhase::getChemPotentials(doublereal* mu) const { getActivityConcentrations(DATA_PTR(m_work)); diff --git a/src/thermo/SingleSpeciesTP.cpp b/src/thermo/SingleSpeciesTP.cpp index 4f243286c..8da09153f 100644 --- a/src/thermo/SingleSpeciesTP.cpp +++ b/src/thermo/SingleSpeciesTP.cpp @@ -100,7 +100,6 @@ doublereal SingleSpeciesTP::cp_mole() const * function in ThermoPhase. However, the standard * state heat capacity will do fine here for now. */ - //getPartialMolarCp(&cpbar); getCp_R(&cpbar); cpbar *= GasConstant; return cpbar; diff --git a/src/thermo/SpeciesThermoFactory.cpp b/src/thermo/SpeciesThermoFactory.cpp index 1eff0c530..e94c40e77 100644 --- a/src/thermo/SpeciesThermoFactory.cpp +++ b/src/thermo/SpeciesThermoFactory.cpp @@ -90,8 +90,6 @@ static void getSpeciesThermoTypes(std::vector & spDataNodeList, has_other = 1; } else { has_other = 1; - //throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:", - // spNode->attrib("name"), "missing"); } } else { throw CanteraError("getSpeciesThermoTypes:", @@ -135,7 +133,6 @@ SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(std::vector & s popError(); } if (iother) { - //writelog("returning new GeneralSpeciesThermo"); return new GeneralSpeciesThermo(); } return newSpeciesThermo(NASA*inasa diff --git a/src/thermo/SpeciesThermoInterpType.cpp b/src/thermo/SpeciesThermoInterpType.cpp index e3092e5c5..193f81f45 100644 --- a/src/thermo/SpeciesThermoInterpType.cpp +++ b/src/thermo/SpeciesThermoInterpType.cpp @@ -140,7 +140,6 @@ void STITbyPDSS::updatePropertiesTemp(const doublereal temp, doublereal* h_RT, doublereal* s_R) const { - //m_vpssmgr_ptr->setState_T(temp); m_PDSS_ptr->setTemperature(temp); AssertThrowMsg(m_index != npos, "STITbyPDSS::updatePropertiesTemp", "object was probably not installed correctly"); diff --git a/src/thermo/StatMech.cpp b/src/thermo/StatMech.cpp index 94213592e..32ead3a64 100644 --- a/src/thermo/StatMech.cpp +++ b/src/thermo/StatMech.cpp @@ -563,7 +563,6 @@ void StatMech::updateProperties(const doublereal* tt, if (name_map.find(sp_name) != name_map.end()) { s = name_map.find(sp_name)->second; } else { - //std::cout << sp_name << std::endl; throw CanteraError("StatMech.cpp", "species properties not found!. \n\n"); } @@ -628,7 +627,6 @@ void StatMech::reportParameters(size_t& n, int& type, if (name_map.find(sp_name) != name_map.end()) { s = name_map.find(sp_name)->second; } else { - //std::cout << sp_name << std::endl; throw CanteraError("StatMech.cpp", "species properties not found!. \n\n"); } diff --git a/src/thermo/VPSSMgr.cpp b/src/thermo/VPSSMgr.cpp index ddae62e46..36dd03a18 100644 --- a/src/thermo/VPSSMgr.cpp +++ b/src/thermo/VPSSMgr.cpp @@ -50,8 +50,6 @@ VPSSMgr::VPSSMgr(const VPSSMgr& right) : m_kk(0), m_vptp_ptr(0), m_spthermo(0), - // m_Tnow(300.), - // m_Pnow(OneAtm), m_tlast(-1.0), m_plast(-1.0), m_p0(-1.0), @@ -404,36 +402,6 @@ void VPSSMgr::initThermoXML(XML_Node& phaseNode, const std::string& id) m_maxTemp = mint; } } -#ifdef DEBUG_MODE - // Add a check to see that all references pressures are the same - double m_p0_k; - if (m_spthermo) { - for (size_t k = 0; k < m_kk; k++) { - m_p0_k = m_spthermo->refPressure(k); - if (m_p0 != m_p0_k) { - //throw CanteraError("VPSSMgr::initThermoXML", - // "inconsistent ref pressures" + fp2str(m_p0) + " " - // + fp2str(m_p0_k)); - // writelog("VPSSMgr::initThermoXML:" - // "inconsistent ref pressures: " + fp2str(m_p0) + " " - // + fp2str(m_p0_k) + " for SpeciesThermo k = " + int2str(k) + "\n"); - } - } - } - - for (size_t k = 0; k < m_kk; k++) { - const PDSS* kPDSS = m_vptp_ptr->providePDSS(k); - m_p0_k = kPDSS->refPressure(); - if (m_p0 != m_p0_k) { - //throw CanteraError("VPSSMgr::initThermoXML", - // "inconsistent ref pressures" + fp2str(m_p0) + " " - // + fp2str(m_p0_k)); - //writelog("VPSSMgr::initThermoXML" - // "inconsistent ref pressures: " + fp2str(m_p0) + " " - // + fp2str(m_p0_k) + " for PDSS k = " + int2str(k) + "\n"); - } - } -#endif } void VPSSMgr::installSTSpecies(size_t k, const XML_Node& s, diff --git a/src/thermo/VPSSMgr_ConstVol.cpp b/src/thermo/VPSSMgr_ConstVol.cpp index 1f3e09c29..cc27beef1 100644 --- a/src/thermo/VPSSMgr_ConstVol.cpp +++ b/src/thermo/VPSSMgr_ConstVol.cpp @@ -131,7 +131,6 @@ PDSS* VPSSMgr_ConstVol::createInstallPDSS(size_t k, const XML_Node& speciesNode, const XML_Node* const phaseNode_ptr) { - //VPSSMgr::installSpecies(k, speciesNode, phaseNode_ptr); const XML_Node* ss = speciesNode.findByName("standardState"); if (!ss) { throw CanteraError("VPSSMgr_ConstVol::installSpecies", diff --git a/src/thermo/VPSSMgr_General.cpp b/src/thermo/VPSSMgr_General.cpp index 0c09558ac..85a1eef63 100644 --- a/src/thermo/VPSSMgr_General.cpp +++ b/src/thermo/VPSSMgr_General.cpp @@ -164,7 +164,6 @@ VPSSMgr_General::returnPDSS_ptr(size_t k, const XML_Node& speciesNode, throw CanteraError("VPSSMgr_General::returnPDSS_ptr", "new PDSS_ConstVol failed"); } } else if (model == "waterIAPWS" || model == "waterPDSS") { - // VPSSMgr::installSTSpecies(k, speciesNode, phaseNode_ptr); kPDSS = new PDSS_Water(m_vptp_ptr, 0); if (!genSpthermo) { throw CanteraError("VPSSMgr_General::returnPDSS_ptr", @@ -213,7 +212,6 @@ VPSSMgr_General::createInstallPDSS(size_t k, const XML_Node& speciesNode, { bool doST; PDSS* kPDSS = returnPDSS_ptr(k, speciesNode, phaseNode_ptr, doST); - // VPSSMgr::installSTSpecies(k, speciesNode, phaseNode_ptr); if (m_PDSS_ptrs.size() < k+1) { m_PDSS_ptrs.resize(k+1, 0); } diff --git a/src/thermo/VPSSMgr_IdealGas.cpp b/src/thermo/VPSSMgr_IdealGas.cpp index d0b044302..39084c860 100644 --- a/src/thermo/VPSSMgr_IdealGas.cpp +++ b/src/thermo/VPSSMgr_IdealGas.cpp @@ -92,7 +92,6 @@ PDSS* VPSSMgr_IdealGas::createInstallPDSS(size_t k, const XML_Node& speciesNode, const XML_Node* const phaseNode_ptr) { - //VPSSMgr::installSpecies(k, speciesNode, phaseNode_ptr); const XML_Node* ss = speciesNode.findByName("standardState"); if (ss) { std::string model = (*ss)["model"]; diff --git a/src/thermo/VPStandardStateTP.cpp b/src/thermo/VPStandardStateTP.cpp index 37e12d349..070dc2aa5 100644 --- a/src/thermo/VPStandardStateTP.cpp +++ b/src/thermo/VPStandardStateTP.cpp @@ -356,7 +356,6 @@ void VPStandardStateTP::initThermoXML(XML_Node& phaseNode, const std::string& id // routine. VPStandardStateTP::initLengths(); - //m_VPSS_ptr->initThermo(); for (size_t k = 0; k < m_kk; k++) { PDSS* kPDSS = m_PDSS_storage[k]; AssertTrace(kPDSS != 0); diff --git a/src/thermo/WaterProps.cpp b/src/thermo/WaterProps.cpp index eb0d77e11..d72ae1fe9 100644 --- a/src/thermo/WaterProps.cpp +++ b/src/thermo/WaterProps.cpp @@ -202,16 +202,11 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc) doublereal psat = satPressure(T); doublereal P; if (psat > P_input) { - //printf("ADebye WARNING: p_input < psat: %g %g\n", - // P_input, psat); P = psat; } else { P = P_input; } doublereal epsRelWater = relEpsilon(T, P, 0); - //printf("releps calc = %g, compare to 78.38\n", epsRelWater); - //doublereal B_Debye = 3.28640E9; - doublereal epsilon = epsilon_0 * epsRelWater; doublereal dw = density_IAPWS(T, P); doublereal tmp = sqrt(2.0 * Avogadro * dw / 1000.); @@ -232,24 +227,13 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc) doublereal depsRelWaterdT = relEpsilon(T, P, 1); dAdT -= A_Debye * (1.5 * depsRelWaterdT / epsRelWater); - //int methodD = 1; - //doublereal ddwdT = density_T_new(T, P, 1); - // doublereal contrib1 = A_Debye * (0.5 * ddwdT / dw); - /* * calculate d(lnV)/dT _constantP, i.e., the cte */ doublereal cte = coeffThermalExp_IAPWS(T, P); doublereal contrib2 = - A_Debye * (0.5 * cte); - - //dAdT += A_Debye * (0.5 * ddwdT / dw); dAdT += contrib2; -#ifdef DEBUG_HKM - //printf("dAdT = %g, contrib1 = %g, contrib2 = %g\n", - // dAdT, contrib1, contrib2); -#endif - if (ifunc == 1) { return dAdT; } @@ -264,13 +248,6 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc) doublereal d2epsRelWaterdT2 = relEpsilon(T, P, 2); - //doublereal dT = -0.01; - //doublereal TT = T + dT; - //doublereal depsRelWaterdTdel = relEpsilon(TT, P, 1); - //doublereal d2alt = (depsRelWaterdTdel- depsRelWaterdT ) / dT; - //printf("diff %g %g\n",d2epsRelWaterdT2, d2alt); - // HKM -> checks out, i.e., they are the same. - d2AdT2 += 1.5 * (- dAdT * depsRelWaterdT / epsRelWater - A_Debye / epsRelWater * (d2epsRelWaterdT2 - depsRelWaterdT * depsRelWaterdT / epsRelWater)); @@ -279,10 +256,6 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc) doublereal Tdel = T + deltaT; doublereal cte_del = coeffThermalExp_IAPWS(Tdel, P); doublereal dctedT = (cte_del - cte) / Tdel; - - - //doublereal d2dwdT2 = density_T_new(T, P, 2); - doublereal contrib3 = 0.5 * (-(dAdT * cte) -(A_Debye * dctedT)); d2AdT2 += contrib3; @@ -309,8 +282,6 @@ doublereal WaterProps::ADebye(doublereal T, doublereal P_input, int ifunc) dAdP -= A_Debye * (1.5 * depsRelWaterdP / epsRelWater); doublereal kappa = isothermalCompressibility_IAPWS(T,P); - - //doublereal ddwdP = density_T_new(T, P, 3); dAdP += A_Debye * (0.5 * kappa); return dAdP; @@ -375,25 +346,13 @@ doublereal WaterProps::viscosityWater() const doublereal temp = m_waterIAPWS->temperature(); doublereal dens = m_waterIAPWS->density(); - //WaterPropsIAPWS *waterP = new WaterPropsIAPWS(); - //m_waterIAPWS->setState_TR(temp, dens); - //doublereal pressure = m_waterIAPWS->pressure(); - //printf("pressure = %g\n", pressure); - //dens = 18.02 * pressure / (GasConstant * temp); - //printf ("mod dens = %g\n", dens); - doublereal rhobar = dens/rhoStar; doublereal tbar = temp / TStar; - // doublereal pbar = pressure / presStar; - doublereal tbar2 = tbar * tbar; doublereal tbar3 = tbar2 * tbar; doublereal mu0bar = std::sqrt(tbar) / (H[0] + H[1]/tbar + H[2]/tbar2 + H[3]/tbar3); - //printf("mu0bar = %g\n", mu0bar); - //printf("mu0 = %g\n", mu0bar * muStar); - doublereal tfac1 = 1.0 / tbar - 1.0; doublereal tfac2 = tfac1 * tfac1; doublereal tfac3 = tfac2 * tfac1; @@ -464,8 +423,6 @@ doublereal WaterProps::thermalConductivityWater() const doublereal tbar3 = tbar2 * tbar; doublereal lambda0bar = sqrt(tbar) / (L[0] + L[1]/tbar + L[2]/tbar2 + L[3]/tbar3); - //doublereal lambdagas = lambda0bar * lambdastar * 1.0E3; - doublereal tfac1 = 1.0 / tbar - 1.0; doublereal tfac2 = tfac1 * tfac1; doublereal tfac3 = tfac2 * tfac1; diff --git a/src/thermo/WaterPropsIAPWS.cpp b/src/thermo/WaterPropsIAPWS.cpp index 57c9d224f..d2a2de2ff 100644 --- a/src/thermo/WaterPropsIAPWS.cpp +++ b/src/thermo/WaterPropsIAPWS.cpp @@ -411,7 +411,6 @@ int WaterPropsIAPWS::phaseState(bool checkState) const } else { doublereal T = T_c / tau; doublereal rho = delta * Rho_c; - //doublereal psatTable = psat_est(T); doublereal rhoMidAtm = 0.5 * (OneAtm * M_water / (Rgas * 373.15) + 1.0E3); doublereal rhoMid = Rho_c + (T - T_c) * (Rho_c - rhoMidAtm) / (T_c - 373.15); int iStateGuess = WATER_LIQUID; @@ -424,8 +423,6 @@ int WaterPropsIAPWS::phaseState(bool checkState) const } else { // When we are here we are between the spinodal curves doublereal rhoDel = rho * 1.000001; - - //setState_TR(T, rhoDel); doublereal deltaSave = delta; doublereal deltaDel = rhoDel / Rho_c; delta = deltaDel; @@ -438,7 +435,6 @@ int WaterPropsIAPWS::phaseState(bool checkState) const } else { iState = WATER_UNSTABLEGAS; } - //setState_TR(T, rho); delta = deltaSave; m_phi->tdpolycalc(tau, delta); @@ -489,7 +485,6 @@ doublereal WaterPropsIAPWS::densSpinodalWater() const slope = std::max(slope, dpdrho_new *5.0/ dens_new); } else { slope = -dpdrho_new; - //slope = MIN(slope, dpdrho_new *5.0 / dens_new); // shouldn't be here for liquid spinodal } doublereal delta_rho = - dpdrho_new / slope; @@ -579,10 +574,8 @@ doublereal WaterPropsIAPWS::densSpinodalSteam() const doublereal slope = (dpdrho_new - dpdrho_old)/(dens_new - dens_old); if (slope >= 0.0) { slope = dpdrho_new; - //slope = MAX(slope, dpdrho_new *5.0/ dens_new); // shouldn't be here for gas spinodal } else { - //slope = -dpdrho_new; slope = std::min(slope, dpdrho_new *5.0 / dens_new); } diff --git a/src/thermo/WaterPropsIAPWSphi.cpp b/src/thermo/WaterPropsIAPWSphi.cpp index ffa379461..5e996d99f 100644 --- a/src/thermo/WaterPropsIAPWSphi.cpp +++ b/src/thermo/WaterPropsIAPWSphi.cpp @@ -1151,7 +1151,6 @@ doublereal WaterPropsIAPWSphi::cp_R() const doublereal tau = TAUsave; doublereal delta = DELTAsave; doublereal cvR = cv_R(); - //doublereal nd = phi0_d(); doublereal rd = phiR_d(); doublereal rdd = phiR_dd(); doublereal rdt = phiR_dt(); diff --git a/src/thermo/WaterSSTP.cpp b/src/thermo/WaterSSTP.cpp index 8e4d31e97..173717d03 100644 --- a/src/thermo/WaterSSTP.cpp +++ b/src/thermo/WaterSSTP.cpp @@ -170,7 +170,6 @@ void WaterSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id) } s = entropy_mole(); s -= GasConstant * log(oneBar/presLow); - //printf("s = %g\n", s); doublereal h = enthalpy_mole(); if (h != -241.826E6) { @@ -178,9 +177,6 @@ void WaterSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id) } h = enthalpy_mole(); - //printf("h = %g\n", h); - - /* * Set the initial state of the system to 298.15 K and * 1 bar. diff --git a/src/tpx/CarbonDioxide.cpp b/src/tpx/CarbonDioxide.cpp index 45d453bbe..6abb5b03d 100644 --- a/src/tpx/CarbonDioxide.cpp +++ b/src/tpx/CarbonDioxide.cpp @@ -212,7 +212,6 @@ double CarbonDioxide::up() double CarbonDioxide::sp() { - //double Tinverse = 1.0/T; double T2inverse = pow(T, -2); double T3inverse = pow(T, -3); double T4inverse = pow(T, -4); @@ -267,8 +266,6 @@ double CarbonDioxide::Psat() log = ((Tc/T)-1)*sum; P=exp(log)*Pc; - - //cout << "Psat is returning " << P << " at T " << T << " and Pc " << Pc << " and Tp " << Tp << endl; return P; } diff --git a/src/tpx/HFC134a.cpp b/src/tpx/HFC134a.cpp index 3a9b7a6c5..fcdae9f8a 100644 --- a/src/tpx/HFC134a.cpp +++ b/src/tpx/HFC134a.cpp @@ -165,19 +165,6 @@ double HFC134a::Psat() return Pc*exp(f/x1); } -/* - double HFC134a::dPsatdT(){ - if ((T < Tmn) || (T > Tc)) set_Err(TempError); - double x1 = T/Tc; - double x2 = 1.0 - x1; - double f = -7.686556*x2 + 2.311791*pow(x2,1.5) - - 2.039554*x2*x2 - 3.583758*pow(x2,4); - double fp = -7.686556 + 1.5*2.311791*pow(x2,0.5) - - 2.0*2.039554*x2 - 4.0*3.583758*pow(x2,3); - return -Pc*exp(f/x1)*(fp/T + f/(x1*T)); - } -*/ - double HFC134a::ldens() { if ((T < Tmn) || (T > Tc)) { diff --git a/src/tpx/Hydrogen.cpp b/src/tpx/Hydrogen.cpp index 40ea9d696..992bf6730 100644 --- a/src/tpx/Hydrogen.cpp +++ b/src/tpx/Hydrogen.cpp @@ -206,7 +206,6 @@ double hydrogen::Pp() { double rt = 1.0/T; double rt2 = rt*rt; - // double rt3 = rt*rt2; double egrho = exp(-Gamma*Rho*Rho); double P = Rho*R*T; diff --git a/src/tpx/Oxygen.cpp b/src/tpx/Oxygen.cpp index f2b38e1fd..4d05892e6 100644 --- a/src/tpx/Oxygen.cpp +++ b/src/tpx/Oxygen.cpp @@ -180,7 +180,6 @@ double oxygen::Pp() { double rt = 1.0/T; double rt2 = rt*rt; - //double rt3 = rt*rt2; double egrho = exp(-Gamma*Rho*Rho); double P = Rho*R*T; diff --git a/src/tpx/RedlichKwong.cpp b/src/tpx/RedlichKwong.cpp index 421f344ac..9860f2f05 100644 --- a/src/tpx/RedlichKwong.cpp +++ b/src/tpx/RedlichKwong.cpp @@ -28,7 +28,6 @@ double RedlichKwong::sp() { const double Pref = 101325.0; double rgas = 8314.3/m_mw; - //double ss = rgas*(log(Pref/(Rho*rgas*T))); double sr = sresid(); double p = Pp(); return rgas*(log(Pref/p)) + sr + m_entropy_offset; diff --git a/src/tpx/Water.cpp b/src/tpx/Water.cpp index 5369c455c..379835e6b 100644 --- a/src/tpx/Water.cpp +++ b/src/tpx/Water.cpp @@ -168,21 +168,6 @@ double water::Psat() return exp(log)*Pc; } -/* -double water::dPsatdT(){ - double log, sum1=0, sum2=0; - int i; - if ((T < Tmn) || (T > Tc)) - set_Err(TempError); // Error("water::dPsatdT",TempError,T); - for (i=1;i<=8;i++) - sum1 += F[i-1]*pow(a*(T-Tp),double(i-1)); - for (i=2;i<=8;i++) - sum2 += F[i-1]*a*(i-1)*pow(a*(T-Tp),double(i-2)); - log = (Tc/T-1)*sum2 - Tc*sum1/(T*T); - return log*Psat(); -} -*/ - double water::ldens() { double sum=0; diff --git a/src/transport/AqueousTransport.cpp b/src/transport/AqueousTransport.cpp index 814b986b3..af8b3a107 100644 --- a/src/transport/AqueousTransport.cpp +++ b/src/transport/AqueousTransport.cpp @@ -56,9 +56,6 @@ bool AqueousTransport::initLiquid(LiquidTransportParams& tr) m_thermo->molecularWeights().end(), m_mw.begin()); // copy polynomials and parameters into local storage - //m_visccoeffs = tr.visccoeffs; - //m_condcoeffs = tr.condcoeffs; - //m_diffcoeffs = tr.diffcoeffs; cout << "In AqueousTransport::initLiquid we need to replace" << endl << "LiquidTransportParams polynomial coefficients with" << endl << "those in LiquidTransportData as in SimpleTransport." << endl; @@ -336,15 +333,6 @@ void AqueousTransport::update_C() { doublereal pres = m_thermo->pressure(); - // Check for changes in the mole fraction vector. - //int iStateNew = m_thermo->getIStateMF(); - //if (iStateNew == m_iStateMF) { - // if (pres == m_press) { - // return; - // } - // } else { - // m_iStateMF = iStateNew; - //} m_press = pres; // signal that concentration-dependent quantities will need to @@ -463,16 +451,6 @@ void AqueousTransport::stefan_maxwell_solve() // grab a local copy of the molecular weights const vector_fp& M = m_thermo->molecularWeights(); - - // get the mean molecular weight of the mixture - //double M_mix = m_thermo->meanMolecularWeight(); - - - // get the concentration of the mixture - //double rho = m_thermo->density(); - //double c = rho/M_mix; - - m_thermo->getMoleFractions(DATA_PTR(m_molefracs)); double T = m_thermo->temperature(); @@ -534,9 +512,6 @@ void AqueousTransport::stefan_maxwell_solve() } } - //! invert and solve the system Ax = b. Answer is in m_B - //solve(m_A, m_B); - m_flux = m_B; @@ -563,9 +538,6 @@ void AqueousTransport::stefan_maxwell_solve() } } - //! invert and solve the system Ax = b. Answer is in m_B - //solve(m_A, m_B); - m_flux = m_B; diff --git a/src/transport/LiquidTranInteraction.cpp b/src/transport/LiquidTranInteraction.cpp index bf1a64e47..7f6ca445b 100644 --- a/src/transport/LiquidTranInteraction.cpp +++ b/src/transport/LiquidTranInteraction.cpp @@ -74,18 +74,6 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, size_t nsp = thermo->nSpecies(); m_Dij.resize(nsp, nsp, 0.0); m_Eij.resize(nsp, nsp, 0.0); - /* - m_Aij.resize(nsp); - m_Bij.resize(nsp); - m_Hij.resize(nsp); - m_Sij.resize(nsp); - for (int k = 0; k < nsp; k++ ){ - (*m_Aij[k]).resize(nsp, nsp, 0.0); - (*m_Bij[k]).resize(nsp, nsp, 0.0); - (*m_Hij[k]).resize(nsp, nsp, 0.0); - (*m_Sij[k]).resize(nsp, nsp, 0.0); - } - */ std::string speciesA; std::string speciesB; @@ -110,10 +98,6 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, throw CanteraError("TransportFactory::getLiquidInteractionsTransportData", "Unknown species " + speciesB); } - /* if (xmlChild.hasChild("Aij" ) ) { - m_Aij(iSpecies,jSpecies) = getFloat(xmlChild, "Aij", "toSI" ); - m_Aij(jSpecies,iSpecies) = m_Aij(iSpecies,jSpecies) ; - }*/ if (xmlChild.hasChild("Eij")) { m_Eij(iSpecies,jSpecies) = getFloat(xmlChild, "Eij", "actEnergy"); @@ -123,9 +107,7 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, if (xmlChild.hasChild("Aij")) { vector_fp poly; - // poly0 = getFloat(poly, xmlChild, "Aij", "toSI" ); getFloatArray(xmlChild, poly, true, "toSI", "Aij"); - // if (!poly.size() ) poly.push_back(poly0); while (m_Aij.size()resize(nsp, nsp, 0.0); @@ -133,14 +115,12 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, } for (int i = 0; i < (int)poly.size(); i++) { (*m_Aij[i])(iSpecies,jSpecies) = poly[i]; - //(*m_Aij[i])(jSpecies,iSpecies) = (*m_Aij[i])(iSpecies,jSpecies) ; } } if (xmlChild.hasChild("Bij")) { vector_fp poly; getFloatArray(xmlChild, poly, true, "toSI", "Bij"); - //if (!poly.size() ) poly.push_back(poly0); while (m_Bij.size() < poly.size()) { DenseMatrix* bTemp = new DenseMatrix(); bTemp->resize(nsp, nsp, 0.0); @@ -148,15 +128,12 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, } for (size_t i=0; iresize(nsp, nsp, 0.0); @@ -165,15 +142,12 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, for (size_t i=0; iresize(nsp, nsp, 0.0); @@ -182,16 +156,9 @@ void LiquidTranInteraction::init(const XML_Node& compModelNode, for (size_t i=0; i LTPptrs) doublereal value = 0.0; for (size_t k = 0; k < nsp; k++) { - //molefracs[k] = molefracs[k]; // should be: molefracs[k] = molefracs[k]*LTPptrs[k]->getMixWeight(); for consistency, but weight(solvent)=1? } @@ -455,11 +413,9 @@ doublereal LTI_Log_MoleFracs::getMixTransProp(doublereal* speciesValues, doubler for (size_t j = 0; j < nsp; j++) { for (size_t k = 0; k < m_Hij.size(); k++) { value += molefracs[i]*molefracs[j]*(*m_Hij[k])(i,j)/temp*pow(molefracs[i], (int) k); - //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; } for (size_t k = 0; k < m_Sij.size(); k++) { value -= molefracs[i]*molefracs[j]*(*m_Sij[k])(i,j)*pow(molefracs[i], (int) k); - //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; } } } @@ -488,21 +444,14 @@ doublereal LTI_Log_MoleFracs::getMixTransProp(std::vector LTPptrs) for (size_t j = 0; j < nsp; j++) { for (size_t k = 0; k < m_Hij.size(); k++) { value += molefracs[i]*molefracs[j]*(*m_Hij[k])(i,j)/temp*pow(molefracs[i], (int) k); - //cout << "1 = " << molefracs[i]+molefracs[j] << endl; - //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; } for (size_t k = 0; k < m_Sij.size(); k++) { value -= molefracs[i]*molefracs[j]*(*m_Sij[k])(i,j)*pow(molefracs[i], (int) k); - //cout << "1 = " << molefracs[i]+molefracs[j] << endl; - //cout << "value = " << value << ", m_Sij = " << (*m_Sij[k])(i,j) << ", m_Hij = " << (*m_Hij[k])(i,j) << endl; } } } value = exp(value); - // cout << ", viscSpeciesA = " << LTPptrs[0]->getSpeciesTransProp() << endl; - //cout << ", viscSpeciesB = " << LTPptrs[1]->getSpeciesTransProp() << endl; - //cout << "value = " << value << " FINAL" << endl; return value; } @@ -568,11 +517,9 @@ void LTI_StefanMaxwell_PPN::setParameters(LiquidTransportParams& trParam) { size_t nsp = m_thermo->nSpecies(); size_t nsp2 = nsp*nsp; - //vectorgetMixTransProp(m_ionCondTempDep_Ns); return m_ionCondmix; - - /* - // update m_ionCondSpecies[] if necessary - if (!m_ionCond_temp_ok) { - updateIonConductivity_T(); - } - - if (!m_ionCond_conc_ok) { - updateIonConductivity_C(); - } - */ } void LiquidTransport::getSpeciesIonConductivity(doublereal* ionCond) @@ -570,9 +557,6 @@ void LiquidTransport::getBinaryDiffCoeffs(size_t ld, doublereal* d) for (size_t i = 0; i < m_nsp; i++) { for (size_t j = 0; j < m_nsp; j++) { - //if (!( ( m_bdiff(i,j) > 0.0 ) | ( m_bdiff(i,j) < 0.0 ))){ - // throw CanteraError("LiquidTransport::getBinaryDiffCoeffs ", - // "m_bdiff has zero entry in non-diagonal.");} d[ld*j + i] = 1.0 / m_bdiff(i,j); } @@ -955,12 +939,6 @@ void LiquidTransport::updateHydrodynamicRadius_T() void LiquidTransport::update_Grad_lnAC() { doublereal grad_T; - // static vector_fp grad_lnAC(m_nsp), grad_X(m_nsp); - // IonsFromNeutralVPSSTP * tempIons = dynamic_cast m_thermo; - //MargulesVPSSTP * tempMarg = dynamic_cast (tempIons->neutralMoleculePhase_); - - - //m_thermo->getdlnActCoeffdlnX( DATA_PTR(grad_lnAC) ); for (size_t k = 0; k < m_nDim; k++) { grad_T = m_Grad_T[k]; const int start = m_nsp*k; @@ -971,7 +949,6 @@ void LiquidTransport::update_Grad_lnAC() } else { m_Grad_lnAC[start+i] += m_Grad_X[start+i]/m_molefracs[i]; } - // std::cout << k << " m_Grad_lnAC = " << m_Grad_lnAC[k] << std::endl; } return; @@ -1001,7 +978,6 @@ void LiquidTransport::stefan_maxwell_solve() update_Grad_lnAC() ; - //m_thermo->getStandardVolumes(DATA_PTR(m_volume_spec)); m_thermo->getActivityCoefficients(DATA_PTR(m_actCoeff)); /* @@ -1033,7 +1009,6 @@ void LiquidTransport::stefan_maxwell_solve() for (size_t i = 0; i < m_nsp; i++) { m_Grad_mu[a*m_nsp + i] = m_chargeSpecies[i] * Faraday * m_Grad_V[a] - //+ (m_volume_spec[i] - M[i]/dens_) * m_Grad_P[a] + GasConstant * T * m_Grad_lnAC[a*m_nsp+i]; } } @@ -1085,9 +1060,6 @@ void LiquidTransport::stefan_maxwell_solve() m_A(i,i) = 0.0; for (size_t j = 0; j < m_nsp; j++) { if (j != i) { - //if ( !( m_bdiff(i,j) > 0.0 ) ) - //throw CanteraError("LiquidTransport::stefan_maxwell_solve", - // "m_bdiff has zero entry in non-diagonal."); tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = tmp; @@ -1098,29 +1070,10 @@ void LiquidTransport::stefan_maxwell_solve() //! invert and solve the system Ax = b. Answer is in m_B solve(m_A, m_B); - /* - condSum2 = m_chargeSpecies[1]*m_chargeSpecies[1]*m_molefracs_tran[1]*m_bdiff(2,3) + - m_chargeSpecies[2]*m_chargeSpecies[2]*m_molefracs_tran[2]*m_bdiff(1,3) + - m_chargeSpecies[3]*m_chargeSpecies[3]*m_molefracs_tran[3]*m_bdiff(1,2); - condSum1 = m_molefracs_tran[1]*m_bdiff(1,2)*m_bdiff(1,3) + - m_molefracs_tran[2]*m_bdiff(2,3)*m_bdiff(1,2) + - m_molefracs_tran[3]*m_bdiff(1,3)*m_bdiff(2,3); - condSum2 = condSum2/condSum1*Faraday*Faraday/GasConstant/T/vol; - */ - condSum1 = 0; for (size_t i = 0; i < m_nsp; i++) { condSum1 -= Faraday*m_chargeSpecies[i]*m_B(i,0)*m_molefracs_tran[i]/vol; } - - /* - Check Mobility Ratio of Cations - cout << "mobility ratio = " << m_chargeSpecies[1]*(m_B(1,0)-m_B(2,0))/m_chargeSpecies[0]/(m_B(0,0)-m_B(2,0)) << endl; - */ - - // cout << condSum1 << " = " << condSum2 << endl; - - break; case 2: /* 2-D approximation */ m_B(0,0) = 0.0; @@ -1149,9 +1102,6 @@ void LiquidTransport::stefan_maxwell_solve() m_A(i,i) = 0.0; for (size_t j = 0; j < m_nsp; j++) { if (j != i) { - //if ( !( m_bdiff(i,j) > 0.0 ) ) - //throw CanteraError("LiquidTransport::stefan_maxwell_solve", - // "m_bdiff has zero entry in non-diagonal."); tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = tmp; @@ -1194,9 +1144,6 @@ void LiquidTransport::stefan_maxwell_solve() m_A(i,i) = 0.0; for (size_t j = 0; j < m_nsp; j++) { if (j != i) { - //if ( !( m_bdiff(i,j) > 0.0 ) ) - //throw CanteraError("LiquidTransport::stefan_maxwell_solve", - // "m_bdiff has zero entry in non-diagonal."); tmp = m_molefracs_tran[j] * m_bdiff(i,j); m_A(i,i) -= tmp; m_A(i,j) = tmp; diff --git a/src/transport/MultiTransport.cpp b/src/transport/MultiTransport.cpp index 0c1bd0dee..656732963 100644 --- a/src/transport/MultiTransport.cpp +++ b/src/transport/MultiTransport.cpp @@ -188,19 +188,11 @@ void MultiTransport::solveLMatrixEquation() // in m_a should provide a good starting guess, so convergence // should be fast. - //if (m_gmres) { - // gmres(m_mgmres, 3*m_nsp, m_Lmatrix, m_b.begin(), - // m_a.begin(), m_eps_gmres); - // m_lmatrix_soln_ok = true; - // m_l0000_ok = true; // L matrix not modified by GMRES - //} - //else { copy(m_b.begin(), m_b.end(), m_a.begin()); try { solve(m_Lmatrix, DATA_PTR(m_a)); } catch (CanteraError& err) { err.save(); - //if (info != 0) { throw CanteraError("MultiTransport::solveLMatrixEquation", "error in solving L matrix."); } @@ -263,11 +255,6 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ for (size_t n = 0; n < ldx*ndim; n++) { grx[n] = grad_X[n]; } - //for (n = 0; n < ndim; n++) { - // gsave[n] = grad_X[jmax + n*ldx]; // save the input mole frac gradient - //grad_X[jmax + n*ldx] = 0.0; - // grx[jmax + n*ldx] = 0.0; - // } // copy grad_X to fluxes const doublereal* gx; @@ -300,7 +287,6 @@ void MultiTransport::getSpeciesFluxes(size_t ndim, const doublereal* const grad_ for (size_t i = 0; i < m_nsp; i++) { fluxes[i + offset] *= rho * y[i] / pp; } - //grad_X[jmax + n*ldx] = gsave[n]; } // thermal diffusion @@ -453,7 +439,6 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d) m_l0000_ok = false; // matrix is overwritten by inverse m_lmatrix_soln_ok = false; - //doublereal pres = m_thermo->pressure(); doublereal prefactor = 16.0 * m_temp * m_thermo->meanMolecularWeight()/(25.0 * p); doublereal c; @@ -563,7 +548,6 @@ void MultiTransport::updateThermal_T() m_cinternal[k] = cp[k] - 2.5; } - // m_thermo->update_T(m_update_thermal_T); m_thermal_tlast = m_thermo->temperature(); } @@ -604,7 +588,6 @@ void MultiTransport::eval_L0010(const doublereal* const x) doublereal sum, wj, xj; for (size_t j = 0; j < m_nsp; j++) { - //constant = prefactor * x[j]; xj = x[j]; wj = m_mw[j]; sum = 0.0; diff --git a/src/transport/PecosTransport.cpp b/src/transport/PecosTransport.cpp index 5387a11e9..37a2df9c5 100755 --- a/src/transport/PecosTransport.cpp +++ b/src/transport/PecosTransport.cpp @@ -427,10 +427,6 @@ void PecosTransport::updateDiff_T() void PecosTransport::updateSpeciesViscosities() { - - // blottner - // return 0.10*std::exp(_a*(logT*logT) + _b*logT + _c); - int k; // iterate over species, update pure-species viscosity for (k = 0; k < m_nsp; k++) { @@ -444,28 +440,6 @@ void PecosTransport::updateSpeciesViscosities() void PecosTransport::read_blottner_transport_table() { - // istringstream blot - // ("Air 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" - // "CPAir 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" - // "N 1.15572000000e-02 6.03167900000e-01 -1.24327495000e+01\n" - // "N2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" - // "CPN2 2.68142000000e-02 3.17783800000e-01 -1.13155513000e+01\n" - // "NO 4.36378000000e-02 -3.35511000000e-02 -9.57674300000e+00\n" - // "O 2.03144000000e-02 4.29440400000e-01 -1.16031403000e+01\n" - // "O2 4.49290000000e-02 -8.26158000000e-02 -9.20194750000e+00\n" - // "C -8.3285e-3 0.7703240 -12.7378000\n" - // "C2 -8.4311e-3 0.7876060 -13.0268000\n" - // "C3 -8.4312e-3 0.7876090 -12.8240000\n" - // "C2H -2.4241e-2 1.0946550 -14.5835500\n" - // "CN -8.3811e-3 0.7860330 -12.9406000\n" - // "CO -0.019527394 1.013295 -13.97873\n" - // "CO2 -0.019527387 1.047818 -14.32212\n" - // "HCN -2.4241e-2 1.0946550 -14.5835500\n" - // "H -8.3912e-3 0.7743270 -13.6653000\n" - // "H2 -8.3346e-3 0.7815380 -13.5351000\n" - // "e 0.00000000000e+00 0.00000000000e+00 -1.16031403000e+01\n"); - - // // from: AIAA-1997-2474 and Sandia Report SC-RR-70-754 // // # Air -- Identical to N2 fit @@ -544,23 +518,6 @@ void PecosTransport::read_blottner_transport_table() } // done with for loop } - - - // for (k = 0; k < m_nsp; k++) - // { - // string sss = m_thermo->speciesName(k); - // cout << sss << endl; - // cout << a[k] << endl; - // cout << b[k] << endl; - // cout << c[k] << endl; - // } - - // simple sanity check - // if(i != m_nsp-1) - // { - // std::cout << "error\n" << i << std::endl; - // } - } void PecosTransport::updateViscosity_T() diff --git a/src/transport/SimpleTransport.cpp b/src/transport/SimpleTransport.cpp index dca6dca18..ff69ae0f7 100644 --- a/src/transport/SimpleTransport.cpp +++ b/src/transport/SimpleTransport.cpp @@ -192,39 +192,11 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr) m_coeffVisc_Ns.clear(); m_coeffVisc_Ns.resize(m_nsp); - //Cantera::LiquidTransportData <d0 = tr.LTData[0]; std::string spName = m_thermo->speciesName(0); - /* - LiquidTR_Model vm0 = ltd0.model_viscosity; - std::string spName0 = m_thermo->speciesName(0); - if (vm0 == LTR_MODEL_CONSTANT) { - tempDepType_ = 0; - } else if (vm0 == LTR_MODEL_ARRHENIUS) { - tempDepType_ = 1; - } else if (vm0 == LTR_MODEL_NOTSET) { - throw CanteraError("SimpleTransport::initLiquid", - "Viscosity Model is not set for species " + spName0 + " in the input file"); - } else { - throw CanteraError("SimpleTransport::initLiquid", - "Viscosity Model for species " + spName0 + " is not handled by this object"); - } - */ for (size_t k = 0; k < m_nsp; k++) { spName = m_thermo->speciesName(k); Cantera::LiquidTransportData& ltd = tr.LTData[k]; - //LiquidTR_Model vm = ltd.model_viscosity; - //vector_fp &kentry = m_coeffVisc_Ns[k]; - /* - if (vm != vm0) { - if (compositionDepType_ != 0) { - throw CanteraError(" SimpleTransport::initLiquid", - "different viscosity models for species " + spName + " and " + spName0 ); - } else { - kentry = m_coeffVisc_Ns[0]; - } - } - */ m_coeffVisc_Ns[k] = ltd.viscosity; ltd.viscosity = 0; } @@ -235,27 +207,10 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr) m_condSpecies.resize(m_nsp); m_coeffLambda_Ns.clear(); m_coeffLambda_Ns.resize(m_nsp); - //LiquidTR_Model cm0 = ltd0.model_thermalCond; - //if (cm0 != vm0) { - // throw CanteraError("SimpleTransport::initLiquid", - // "Conductivity model is not the same as the viscosity model for species " + spName0); - // } for (size_t k = 0; k < m_nsp; k++) { spName = m_thermo->speciesName(k); Cantera::LiquidTransportData& ltd = tr.LTData[k]; - //LiquidTR_Model cm = ltd.model_thermalCond; - //vector_fp &kentry = m_coeffLambda_Ns[k]; - /* - if (cm != cm0) { - if (compositionDepType_ != 0) { - throw CanteraError(" SimpleTransport::initLiquid", - "different thermal conductivity models for species " + spName + " and " + spName0); - } else { - kentry = m_coeffLambda_Ns[0]; - } - } - */ m_coeffLambda_Ns[k] = ltd.thermalCond; ltd.thermalCond = 0; } @@ -268,52 +223,10 @@ bool SimpleTransport::initLiquid(LiquidTransportParams& tr) m_diffSpecies.resize(m_nsp); m_coeffDiff_Ns.clear(); m_coeffDiff_Ns.resize(m_nsp); - //LiquidTR_Model dm0 = ltd0.model_speciesDiffusivity; - /* - if (dm0 != vm0) { - if (dm0 == LTR_MODEL_NOTSET) { - LiquidTR_Model rm0 = ltd0.model_hydroradius; - if (rm0 != vm0) { - throw CanteraError("SimpleTransport::initLiquid", - "hydroradius model is not the same as the viscosity model for species " + spName0); - } else { - useHydroRadius_ = true; - } - } - } - */ for (size_t k = 0; k < m_nsp; k++) { spName = m_thermo->speciesName(k); Cantera::LiquidTransportData& ltd = tr.LTData[k]; - /* - LiquidTR_Model dm = ltd.model_speciesDiffusivity; - if (dm == LTR_MODEL_NOTSET) { - LiquidTR_Model rm = ltd.model_hydroradius; - if (rm == LTR_MODEL_NOTSET) { - throw CanteraError("SimpleTransport::initLiquid", - "Neither diffusivity nor hydroradius is set for species " + spName); - } - if (rm != vm0) { - throw CanteraError("SimpleTransport::initLiquid", - "hydroradius model is not the same as the viscosity model for species " + spName); - } - if (rm != LTR_MODEL_CONSTANT) { - throw CanteraError("SimpleTransport::initLiquid", - "hydroradius model is not constant for species " + spName0); - } - vector_fp &kentry = m_coeffHydroRadius_Ns[k]; - kentry = ltd.hydroradius; - } else { - if (dm != dm0) { - throw CanteraError(" SimpleTransport::initLiquid", - "different diffusivity models for species " + spName + " and " + spName0 ); - } - vector_fp &kentry = m_coeffDiff_Ns[k]; - kentry = ltd.speciesDiffusivity; - } - */ - m_coeffDiff_Ns[k] = ltd.speciesDiffusivity; ltd.speciesDiffusivity = 0; diff --git a/src/transport/SolidTransport.cpp b/src/transport/SolidTransport.cpp index 4505a8890..6a0842138 100644 --- a/src/transport/SolidTransport.cpp +++ b/src/transport/SolidTransport.cpp @@ -76,15 +76,6 @@ bool SolidTransport::initSolid(SolidTransportData& tr) { m_thermo = tr.thermo; tr.thermo = 0; - //m_nsp = m_thermo->nSpecies(); - //m_tmin = m_thermo->minTemp(); - //m_tmax = m_thermo->maxTemp(); - - // make a local copy of the molecular weights - //m_mw.resize(m_nsp, 0.0); - //copy(m_thermo->molecularWeights().begin(), - // m_thermo->molecularWeights().end(), m_mw.begin()); - m_ionConductivity = tr.ionConductivity; tr.ionConductivity = 0; m_electConductivity = tr.electConductivity; diff --git a/src/transport/TransportFactory.cpp b/src/transport/TransportFactory.cpp index 446f9fb7b..b133d96ba 100644 --- a/src/transport/TransportFactory.cpp +++ b/src/transport/TransportFactory.cpp @@ -178,7 +178,6 @@ TransportFactory::TransportFactory() : m_models["User"] = cUserTransport; m_models["Pecos"] = cPecosTransport; m_models["None"] = None; - //m_models["Radiative"] = cRadiative; for (map::iterator iter = m_models.begin(); iter != m_models.end(); iter++) { @@ -314,8 +313,6 @@ LiquidTranInteraction* TransportFactory::newLTI(const XML_Node& trNode, // lti = new LiquidTranInteraction(tp_ind); lti->init(trNode, thermo); - // throw CanteraError("TransportFactory::newLTI()", - // "unknown Liquid Transport Interaction submodel: " + model ); } return lti; } @@ -563,8 +560,6 @@ void TransportFactory::setupLiquidTransport(std::ostream& flog, thermo_t* thermo // Need to identify a method to obtain interaction matrices. // This will fill LiquidTransportParams members visc_Eij, visc_Sij - // trParam.visc_Eij.resize(nsp,nsp); - // trParam.visc_Sij.resize(nsp,nsp); trParam.thermalCond_Aij.resize(nsp,nsp); trParam.diff_Dij.resize(nsp,nsp); trParam.radius_Aij.resize(nsp,nsp); @@ -602,15 +597,8 @@ void TransportFactory::setupSolidTransport(std::ostream& flog, thermo_t* thermo, copy(trParam.thermo->molecularWeights().begin(), trParam.thermo->molecularWeights().end(), trParam.mw.begin()); - // Resize all other vectors in trParam - //trParam.LTData.resize(nsp); - XML_Node root, log; - // Note that getSolidSpeciesTransportData just populates the pure species transport data. - // const std::vector & species_database = thermo->speciesData(); - // getSolidSpeciesTransportData(species_database, log, trParam.thermo->speciesNames(), trParam); - // getSolidTransportData() populates the // phase transport models like electronic conductivity // thermal conductivity, interstitial diffusion @@ -1148,7 +1136,6 @@ void TransportFactory::getSolidTransportData(const XML_Node& transportNode, ThermoPhase* temp_thermo = trParam.thermo; //tranTypeNode contains the interaction model - // XML_Node &compDepNode = tranTypeNode.child("compositionDependence"); switch (m_tranPropMap[nodeName]) { case TP_IONCONDUCTIVITY: trParam.ionConductivity = newLTP(tranTypeNode, phaseName, @@ -1183,9 +1170,6 @@ void TransportFactory::getSolidTransportData(const XML_Node& transportNode, } catch (CanteraError) { showErrors(std::cout); } - //catch(CanteraError) { - // ; - //} return; } @@ -1440,8 +1424,6 @@ void TransportFactory::fitProperties(GasTransportParams& tr, // NOTE: THIS CORRECTION IS NOT APPLIED doublereal fkj, fjk; getBinDiffCorrection(t, tr, integrals, k, j, 1.0, 1.0, fkj, fjk); - //diffcoeff *= fkj; - if (mode == CK_Mode) { diff[n] = log(diffcoeff); diff --git a/src/zeroD/FlowReactor.cpp b/src/zeroD/FlowReactor.cpp index 297541997..9e5ccdacc 100644 --- a/src/zeroD/FlowReactor.cpp +++ b/src/zeroD/FlowReactor.cpp @@ -50,7 +50,6 @@ void FlowReactor::updateState(doublereal* y) m_dist = y[0]; m_speed = y[1]; doublereal* mss = y + 2; - // doublereal mass = accumulate(y+2, y+2+m_nsp, 0.0); m_thermo->setMassFractions(mss); doublereal rho = m_rho0 * m_speed0/m_speed; diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index 585360480..2a9a8b062 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -21,7 +21,7 @@ ReactorNet::ReactorNet() : Cantera::FuncEval(), m_nr(0), m_nreactors(0), #ifdef DEBUG_MODE m_verbose = true; #endif - m_integ = newIntegrator("CVODE");// CVodeInt; + m_integ = newIntegrator("CVODE"); // use backward differencing, with a full Jacobian computed // numerically, and use a Newton linear iterator