diff --git a/Cantera/clib/src/ctreactor.cpp b/Cantera/clib/src/ctreactor.cpp index b23ee8ad2..3dc4eb8d3 100755 --- a/Cantera/clib/src/ctreactor.cpp +++ b/Cantera/clib/src/ctreactor.cpp @@ -16,7 +16,6 @@ #include "Cabinet.h" #include "Storage.h" -using namespace CanteraZeroD; using namespace Cantera; using namespace std; diff --git a/Cantera/clib/src/ctstagn.cpp b/Cantera/clib/src/ctstagn.cpp deleted file mode 100755 index f32076cc1..000000000 --- a/Cantera/clib/src/ctstagn.cpp +++ /dev/null @@ -1,456 +0,0 @@ - -// Cantera includes -#include "oneD/OneDim.h" -#include "oneD/StFlow.h" -#include "oneD/Inlet1D.h" -#include "oneD/MultiNewton.h" -#include "DenseMatrix.h" -#include "Cabinet.h" -#include "Storage.h" - -// Build as a DLL under Windows -#ifdef WIN32 -#ifdef NO_DLL_BUILD -#define DLL_EXPORT -#else -#define DLL_EXPORT __declspec(dllexport) -#endif -#pragma warning(disable:4786) -#pragma warning(disable:4503) -#else -#define DLL_EXPORT -#endif - -// Values returned for error conditions -#define ERR -999 -#define DERR -999.999 - -using namespace FlowBdry; - -Cabinet* Cabinet::__storage = 0; -Cabinet* Cabinet::__storage = 0; -Cabinet* Cabinet::__storage = 0; -//Cabinet* Cabinet::__storage = 0; - -inline OneDim* _onedim(int i) { - return Cabinet::cabinet()->item(i); -} - -inline StFlow* _flow(int i) { - return Cabinet::cabinet()->item(i); -} - -inline Boundary* _boundary(int i) { - return Cabinet::cabinet()->item(i); -} - -inline Bdry1D* _bndry(int i) { - return Cabinet::cabinet()->item(i); -} - -//inline SurfKinetics* _surfkin(int i) { -// return Cabinet::cabinet()->item(i); -//} - -//inline Surf1D* _surface(int i) { -// return Cabinet::cabinet()->item(i); -//} - -inline DenseMatrix* _matrix(int i) { - return Cabinet::cabinet()->item(i); -} - -inline ThermoPhase* _phase(int n) { - return Storage::__storage->__thtable[n]; -} - -inline Kinetics* _kinetics(int n) { - return Storage::__storage->__ktable[n]; -} - -inline ThermoPhase* _thermo(int n) { - return Storage::__storage->__thtable[n]; -} - -inline Transport* _transport(int n) { - return Storage::__storage->__trtable[n]; -} - - -extern "C" { - - int DLL_EXPORT flow_new(int type, int iph, int np) { - IdealGasPhase* ph = (IdealGasPhase*)_thermo(iph); - StFlow* x; - try { - switch (type) { - case 0: - x = new AxiStagnFlow(ph, ph->nSpecies(), np); break; - case 1: - x = new OneDFlow(ph, ph->nSpecies(), np); break; - default: - return -2; - } - return Cabinet::cabinet()->add(x); - } - catch (CanteraError) { return -1; } - } - - - int DLL_EXPORT flow_del(int i) { - Cabinet::cabinet()->del(i); - return 0; - } - - int DLL_EXPORT flow_copy(int i) { - return Cabinet::cabinet()->newCopy(i); - } - - int DLL_EXPORT flow_assign(int i, int j) { - return Cabinet::cabinet()->assign(i,j); - } - -// int DLL_EXPORT flow_readinputs(int i, char* infile) { -// try { -// ifstream f(infile); -// if (!f) throw CanteraError("flow_readinputs", -// "error opening input file"); -// // _flow(i)->readInputs(f); -// f.close(); -// return 0; -// } -// catch (CanteraError) { return -1; } -// catch (...) { return ERR; } -// } - - int DLL_EXPORT flow_setupgrid(int i, int npts, double* grid) { - try { - _flow(i)->setupGrid(npts, grid); - return 0; - } - catch (CanteraError) { return -1; } - //catch (...) { return ERR; } - } - - int DLL_EXPORT flow_setthermo(int i, int k) { - IdealGasPhase* th = (IdealGasPhase*)_thermo(k); - _flow(i)->setThermo(*th); - return 0; - } - - int DLL_EXPORT flow_setkinetics(int i, int k) { - Kinetics* kin = _kinetics(k); - _flow(i)->setKinetics(*kin); - return 0; - } - - int DLL_EXPORT flow_settransport(int i, int k, int soret) { - try { - Transport* tr = _transport(k); - bool withSoret = (soret == 1); - _flow(i)->setTransport(*tr, withSoret); - return 0; - } - catch (CanteraError) { return -1; } - } - - int DLL_EXPORT flow_settemperature(int i, int j, double t) { - _flow(i)->setTemperature(j, t); - return 0; - } - - int DLL_EXPORT flow_setmassfraction(int i, int j, int k, double t) { - _flow(i)->setMassFraction(j, k, t); - return 0; - } - - int DLL_EXPORT flow_setpressure(int i, double p) { - _flow(i)->setPressure(p); - return 0; - } - - int DLL_EXPORT flow_showsolution(int i, char* fname, double* soln) { - string fn = string(fname); - if (fn == "-") - _flow(i)->showSolution(cout, soln); - else { - ofstream fout(fname); - _flow(i)->showSolution(fout, soln); - fout.close(); - } - return 0; - } - - int DLL_EXPORT flow_outputtec(int i, doublereal* x, - char* fname, char* title, int zone) { - ofstream f(fname); - //DenseMatrix* mat = _matrix(m); - _flow(i)->outputTEC(f, x, string(title), zone); - return 0; - } - - - // solve / fix - - int DLL_EXPORT flow_solveenergyeqn(int i, int j) { - _flow(i)->solveEnergyEqn(j); - return 0; - } - - int DLL_EXPORT flow_fixtemperature(int i, int j) { - _flow(i)->fixTemperature(j); - return 0; - } - - int DLL_EXPORT flow_setenergyfactor(int i, double e) { - _flow(i)->setEnergyFactor(e); - return 0; - } - - int DLL_EXPORT flow_fixspecies(int i, int j) { - _flow(i)->fixSpecies(j); - return 0; - } - - int DLL_EXPORT flow_solvespecies(int i, int j) { - _flow(i)->solveSpecies(j); - return 0; - } - - int DLL_EXPORT flow_resize(int i, int points) { - _flow(i)->resize(points); - return 0; - } - -// int DLL_EXPORT flow_integratechem(int i, doublereal* x, double dt) { -// try{ -// _flow(i)->integrateChem(x, dt); -// return 0; -// } -// catch (CanteraError) { return -1; } -// } - - int DLL_EXPORT flow_settolerances(int i, int nr, - doublereal* rtol, int na, doublereal* atol) { - try { - _flow(i)->setTolerances(nr, rtol, na, atol); - return 0; - } - catch (CanteraError) { return -1; } - //catch (...) { return ERR; } - } - - int DLL_EXPORT flow_eval(int i, int j, doublereal* x, doublereal* r, integer* m) { - try { - _flow(i)->eval(j, x, r, m); - return 0; - } - catch (CanteraError) { return -1; } - } - - int DLL_EXPORT flow_restore(int i, int job, char* fname, char* id, - int& size_z, doublereal* z, int& size_soln, doublereal* soln) { - try { - _flow(i)->restore(job, fname, string(id), size_z, z, - size_soln, soln); - return 0; - } - catch (CanteraError) { return -1; } - catch (...) { return ERR; } - } - - int DLL_EXPORT flow_setfixedpoint(int i, int j0, doublereal t0) { - _flow(i)->setFixedPoint(j0, t0); - return 0; - } - - - int DLL_EXPORT flow_setboundaries(int i, int nleft, int nright) { - Boundary *left=0, *right=0; - if (nleft > 0) left = _boundary(nleft); - if (nright > 0) right = _boundary(nright); - _flow(i)->setBoundaries(left, right); - return 0; - } - - - //========================================================== - - int DLL_EXPORT bdry_new(int type, int iph, int kin) { - Boundary* x=0; - //const doublereal* wt = _phase(iph)->molecularWeights().begin(); - int nsp = _phase(iph)->nSpecies(); - switch (type) { - case 0: - x = new Inlet(nsp); break; - case 1: - x = new Outlet(nsp); break; - //case 2: - //if (kin > 0) - // x = new Surface(nsp, _surfkin(kin)); - //else - // x = new Surface(nsp, 0); - //break; - case 3: - x = new SymmPlane(nsp); break; - default: - return -2; - } - return Cabinet::cabinet()->add(x); - } - - - int DLL_EXPORT bdry_del(int i) { - Cabinet::cabinet()->del(i); - return 0; - } - - int DLL_EXPORT bdry_copy(int i) { - return Cabinet::cabinet()->newCopy(i); - } - - int DLL_EXPORT bdry_assign(int i, int j) { - return Cabinet::cabinet()->assign(i,j); - } - - int DLL_EXPORT bdry_set(int i, int n, doublereal* v) { - switch (n) { - case 1: - _boundary(i)->set_mdot(*v); break; - case 2: - _boundary(i)->set_V(*v); break; - case 3: - _boundary(i)->set_T(*v); break; - case 4: - _boundary(i)->set_Y(v); break; - default: - throw CanteraError("bdry_set","unknown option"); - } - return 0; - } - - //========================================================= - - - int DLL_EXPORT onedim_new(int nd, int* domains, int* types) { - int i; - vector doms; - for (i = 0; i < nd; i++) { - switch (types[i]) { - case 0: - doms.push_back(_flow(domains[i])); break; - //case 1: - //doms.push_back(_surface(domains[i])); break; - case 2: - doms.push_back(_bndry(domains[i])); break; - default: - throw CanteraError("onedim_new", "unknown domain type"); - } - } - try { - OneDim* x = new OneDim(doms); - return Cabinet::cabinet()->add(x); - } - catch (CanteraError) { return -1; } - } - - - int DLL_EXPORT onedim_del(int i) { - Cabinet::cabinet()->del(i); - return 0; - } - - int DLL_EXPORT onedim_addFlow(int i, int n) { - try { - _onedim(i)->addDomain(_flow(n)); - return 0; - } - catch (CanteraError) { return -1; } - // catch (...) { return ERR; } - } - -// int DLL_EXPORT onedim_addSurf(int i, int n) { -// try { -// _onedim(i)->addDomain(_surface(n)); -// return 0; -// } -// catch (CanteraError) { return -1; } -// } - - int DLL_EXPORT onedim_eval(int i, doublereal* x0, doublereal* r) { - try { - _onedim(i)->eval(-1, x0, r, 0.0); - return 0; - } - catch (CanteraError) { return -1; } - // catch (...) { return ERR; } - } - - int DLL_EXPORT onedim_solve(int i, doublereal* x0, doublereal* x1, - int loglevel) { - try { - int m = _onedim(i)->solve(x0, x1, loglevel); - return m; - } - catch (CanteraError) { return -1; } - //catch (...) { return ERR; } - } - - double DLL_EXPORT onedim_ssnorm(int i, doublereal* x0, doublereal* x1) { - return _onedim(i)->ssnorm(x0, x1); - } - - int DLL_EXPORT onedim_setsteadymode(int i) { - if (_onedim(i)->transient()) { - _onedim(i)->setSteadyMode(); - //_onedim(i)->jacobian().setAge(10000); - return 1; - } - return 0; - } - - int DLL_EXPORT onedim_settransientmode(int i, doublereal dt, doublereal* x) { - _onedim(i)->initTimeInteg(dt, x); - double rr = fabs(_onedim(i)->rdt()*dt - 1.0); - if ((rr > 1.e-5) || _onedim(i)->steady()) { - //_onedim(i)->jacobian().setAge(10000); - return 1; - } - return 0; - } - - int DLL_EXPORT onedim_setnewtonoptions(int i, int maxage) { - _onedim(i)->newton().setOptions(maxage); - return 0; - } - - int DLL_EXPORT onedim_resize(int i) { - _onedim(i)->resize(); - return 0; - } - - int DLL_EXPORT onedim_writeStats(int i, int printTime) { - _onedim(i)->writeStats(printTime); - return 0; - } - - double DLL_EXPORT onedim_timestep(int i, int nsteps, doublereal dt, - doublereal* x, doublereal* xnew, int loglevel) { - try { - return _onedim(i)->timeStep(nsteps, dt, x, xnew, loglevel); - } - catch (CanteraError) { return -1.0; } - } - - int DLL_EXPORT onedim_save(int i, char* fname, char* id, - char* desc, doublereal* soln) { - try { - _onedim(i)->save(string(fname), string(id), string(desc), soln); - return 0; - } - catch (CanteraError) { return -1; } - //catch (...) { return ERR; } - } - - -} diff --git a/Cantera/clib/src/ctstagn.h b/Cantera/clib/src/ctstagn.h deleted file mode 100755 index 86742055f..000000000 --- a/Cantera/clib/src/ctstagn.h +++ /dev/null @@ -1,78 +0,0 @@ -/** - * @file ctstagn.h - */ -/* - * $Id$ - */ - -#ifndef CTC_STAGN_H -#define CTC_STAGN_H - -// Cantera includes -//#include "stagn.h" - -//#include "Cabinet.h" -//#include "Storage.h" -#include "clib_defs.h" - -//inline StFlow* _flow(int i) { -// return Cabinet::cabinet()->item(i); -//} - -extern "C" { - - int DLL_IMPORT flow_new(int type, int iph, int np); - int DLL_IMPORT flow_del(int i); - int DLL_IMPORT flow_copy(int i); - int DLL_IMPORT flow_assign(int i, int j); - int DLL_IMPORT flow_setupgrid(int i, int npts, double* grid); - int DLL_EXPORT flow_setthermo(int i, int k); - int DLL_IMPORT flow_setkinetics(int i, int k); - int DLL_IMPORT flow_settransport(int i, int k, int soret); - int DLL_IMPORT flow_solveenergyeqn(int i, int j); - int DLL_IMPORT flow_fixtemperature(int i, int j); - int DLL_IMPORT flow_setenergyfactor(int i, double e); - int DLL_IMPORT flow_fixspecies(int i, int j); - int DLL_IMPORT flow_solvespecies(int i, int j); - // int DLL_IMPORT flow_integratechem(int i, double* x, double dt); - int DLL_IMPORT flow_settemperature(int i, int j, double t); - int DLL_IMPORT flow_setpressure(int i, double p); - int DLL_IMPORT flow_setmassfraction(int i, int j, int k, double t); - int DLL_IMPORT flow_outputtec(int i, double* x, char* fname, - char* title, int zone); - int DLL_IMPORT flow_showsolution(int i, char* fname, double* x); - int DLL_IMPORT flow_settolerances(int i, int nr, - double* rtol, int na, double* atol); - int DLL_IMPORT flow_resize(int i, int points); - int DLL_IMPORT flow_setsteadymode(int i); - int DLL_IMPORT flow_settransientmode(int i, double dt, double* x); - - int DLL_IMPORT flow_restore(int i, int job, char* fname, char* id, - int& size_z, double* z, int& size_soln, double* soln); - int DLL_IMPORT flow_setfixedpoint(int i, int j0, double t0); - int DLL_IMPORT flow_setboundaries(int i, int nleft, int nright); - int DLL_IMPORT bdry_new(int type, int iph, int kin); - int DLL_IMPORT bdry_del(int i); - int DLL_IMPORT bdry_copy(int i); - int DLL_IMPORT bdry_assign(int i, int j); - int DLL_IMPORT bdry_set(int i, int n, double* v); - - int DLL_IMPORT onedim_new(int nd, int* domains, int* types); - int DLL_IMPORT onedim_del(int i); - int DLL_IMPORT onedim_addFlow(int i, int n); - //int DLL_IMPORT onedim_addSurf(int i, int n); - int DLL_EXPORT onedim_eval(int i, double* x0, double* r); - int DLL_IMPORT onedim_solve(int i, double* x0, double* x1, int loglevel); - double DLL_IMPORT onedim_ssnorm(int i, double* x0, double* x1); - int DLL_IMPORT onedim_setsteadymode(int i); - int DLL_IMPORT onedim_settransientmode(int i, double dt, double* x); - int DLL_IMPORT onedim_setnewtonoptions(int i, int maxage); - int DLL_IMPORT onedim_resize(int i); - int DLL_IMPORT onedim_writeStats(int i, int printTime = 1); - double DLL_IMPORT onedim_timestep(int i, int nsteps, double dt, - double* x, double* xnew, int loglevel); - int DLL_IMPORT onedim_save(int i, char* fname, char* id, char* desc, double* soln); - -} - -#endif diff --git a/Cantera/src/equil/BasisOptimize.cpp b/Cantera/src/equil/BasisOptimize.cpp index 79e681d61..9f1406961 100644 --- a/Cantera/src/equil/BasisOptimize.cpp +++ b/Cantera/src/equil/BasisOptimize.cpp @@ -48,7 +48,7 @@ static size_t amax(double *x, size_t j, size_t n); * @param jr first position * @param kspec second species */ -static void switch_pos(vector_int &orderVector, size_t jr, size_t kspec); +static void switch_pos(std::vector &orderVector, size_t jr, size_t kspec); //! Invert an nxn matrix and solve m rhs's diff --git a/Cantera/src/equil/MultiPhase.cpp b/Cantera/src/equil/MultiPhase.cpp index 903cb5f28..38fec09e5 100644 --- a/Cantera/src/equil/MultiPhase.cpp +++ b/Cantera/src/equil/MultiPhase.cpp @@ -286,7 +286,7 @@ namespace Cantera { return sum; } //==================================================================================================================== - int MultiPhase::speciesIndex(std::string speciesName, std::string phaseName) { + size_t MultiPhase::speciesIndex(std::string speciesName, std::string phaseName) { if (!m_init) { init(); } diff --git a/Cantera/src/equil/vcs_VolPhase.cpp b/Cantera/src/equil/vcs_VolPhase.cpp index f09d08743..0a050b8dd 100644 --- a/Cantera/src/equil/vcs_VolPhase.cpp +++ b/Cantera/src/equil/vcs_VolPhase.cpp @@ -502,7 +502,7 @@ namespace VCSnonideal { return Xmol_; } - double vcs_VolPhase::moleFraction(int k) const { + double vcs_VolPhase::moleFraction(size_t k) const { return Xmol_[k]; } /***************************************************************************/ diff --git a/Cantera/src/equil/vcs_phaseStability.cpp b/Cantera/src/equil/vcs_phaseStability.cpp index db557b4af..78a1fde43 100644 --- a/Cantera/src/equil/vcs_phaseStability.cpp +++ b/Cantera/src/equil/vcs_phaseStability.cpp @@ -61,7 +61,7 @@ namespace VCSnonideal { } #endif size_t irxn = kspec - m_numComponents; - if (kspec >= m_numCoimponents) { + if (kspec >= m_numComponents) { bool iPopPossible = true; for (size_t j = 0; j < m_numComponents; ++j) { if (m_elType[j] == VCS_ELEM_TYPE_ABSPOS) { diff --git a/Cantera/src/equil/vcs_rxnadj.cpp b/Cantera/src/equil/vcs_rxnadj.cpp index f98fdc4f0..3885bdfb8 100644 --- a/Cantera/src/equil/vcs_rxnadj.cpp +++ b/Cantera/src/equil/vcs_rxnadj.cpp @@ -40,7 +40,7 @@ namespace VCSnonideal { * in this routine. The species is a noncomponent * - 2 : Same as one but, the zeroed species is a component. */ - int VCS_SOLVE::vcs_RxnStepSizes(int & forceComponentCalc, int &kSpecial) { + int VCS_SOLVE::vcs_RxnStepSizes(int & forceComponentCalc, size_t &kSpecial) { int j, irxn, kspec, iph; int iphDel = -1; double s, xx, dss; diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index b3634bd60..aa1765eb1 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -584,10 +584,12 @@ namespace VCSnonideal { } } else { if (m_doEstimateEquil == 0) { + double sum; for (size_t j = 0; j < nelements; j++) { m_elemAbundancesGoal[j] = 0.0; for (size_t kspec = 0; kspec < nspecies; kspec++) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + sum += m_molNumSpecies_old[kspec]; m_elemAbundancesGoal[j] += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec]; } } diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 593daf3d8..4367114fc 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -490,7 +490,7 @@ public: * the same T and P as the solution. * tg : Total Number of moles in the phase. */ - void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const int ltop); + void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop); //! Print out a table of chemical potentials /*! @@ -577,7 +577,7 @@ public: * * @return Returns an int representing which phase may need to be zeroed */ - int vcs_RxnStepSizes(int & forceComponentCalc, int & kSpecial); + int vcs_RxnStepSizes(int & forceComponentCalc, size_t & kSpecial); //! Calculates the total number of moles of species in all phases. /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index c3ecd40b8..bd1889a53 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -106,6 +106,7 @@ namespace VCSnonideal { int finalElemAbundAttempts = 0; bool uptodate_minors = true; bool justDeletedMultiPhase = false; + bool MajorSpeciesHaveConverged; bool usedZeroedSpecies; /* return flag from basopt indicating that one of the components had a zero concentration */ size_t doPhaseDeleteIph = npos; @@ -224,7 +225,7 @@ namespace VCSnonideal { plogf("% -7.3g ", m_formulaMatrix[j][i]); } plogf(" %3d ", m_phaseID[i]); - print_space(55-m_numElemConstraints*8); + print_space(std::max(55-int(m_numElemConstraints)*8, 0)); plogf("%12.5E %12.5E", RT * m_SSfeSpecies[i], m_molNumSpecies_old[i]); if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { plogf(" Mol_Num"); @@ -563,11 +564,13 @@ namespace VCSnonideal { /********************************************************************/ /************************ VOLTAGE SPECIES ***************************/ /********************************************************************/ + bool soldel_ret; #ifdef DEBUG_MODE - dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); + dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE); #else - dx = vcs_minor_alt_calc(kspec, irxn, &soldel); + dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret); #endif + soldel = soldel_ret; m_deltaMolNumSpecies[kspec] = dx; } else if (m_speciesStatus[kspec] < VCS_SPECIES_MINOR) { @@ -705,11 +708,13 @@ namespace VCSnonideal { * If soldel is true on return, then we branch to the section * that deletes a species from the current set of active species. */ + bool soldel_ret; #ifdef DEBUG_MODE - dx = vcs_minor_alt_calc(kspec, irxn, &soldel, ANOTE); + dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret, ANOTE); #else - dx = vcs_minor_alt_calc(kspec, irxn, &soldel); + dx = vcs_minor_alt_calc(kspec, irxn, &soldel_ret); #endif + soldel = soldel_ret; m_deltaMolNumSpecies[kspec] = dx; m_molNumSpecies_new[kspec] = m_molNumSpecies_old[kspec] + dx; diff --git a/Cantera/src/equil/vcs_solve_phaseStability.cpp b/Cantera/src/equil/vcs_solve_phaseStability.cpp index 73e851416..170e0b994 100644 --- a/Cantera/src/equil/vcs_solve_phaseStability.cpp +++ b/Cantera/src/equil/vcs_solve_phaseStability.cpp @@ -206,7 +206,7 @@ namespace VCSnonideal { int printLvl) { int retn = 0; double test = -1.0E-10; - int usedZeroedSpecies; + bool usedZeroedSpecies; std::vector phasePopPhaseIDs(0); int iphasePop; int iStab = 0; @@ -219,7 +219,7 @@ namespace VCSnonideal { std::vector wx(m_numElemConstraints, 0.0); - retn = vcs_basopt(FALSE, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), + retn = vcs_basopt(false, VCS_DATA_PTR(aw), VCS_DATA_PTR(sa), VCS_DATA_PTR(sm), VCS_DATA_PTR(ss), test, &usedZeroedSpecies); vcs_evaluate_speciesType(); diff --git a/Cantera/src/equil/vcs_util.cpp b/Cantera/src/equil/vcs_util.cpp index df2eb2f49..0070f0c92 100644 --- a/Cantera/src/equil/vcs_util.cpp +++ b/Cantera/src/equil/vcs_util.cpp @@ -274,6 +274,12 @@ namespace VCSnonideal { x[i2] = t; } + void vcsUtil_ssw(size_t x[], size_t i1, size_t i2) { + size_t t = x[i1]; + x[i1] = x[i2]; + x[i2] = t; + } + //==================================================================================================================== #ifdef DEBUG_HKM static void mlequ_matrixDump(double *c, int idem, int n) { @@ -446,7 +452,7 @@ namespace VCSnonideal { * (each column is a new rhs) * @param m number of rhs's */ - int vcsUtil_mlequ(double *c, int idem, int n, double *b, int m) { + int vcsUtil_mlequ(double *c, size_t idem, size_t n, double *b, size_t m) { #ifdef DEBUG_HKM // mlequ_matrixDump(c, idem, n); #endif diff --git a/Cantera/src/kinetics/InterfaceKinetics.h b/Cantera/src/kinetics/InterfaceKinetics.h index 472ac41af..ba41b9198 100644 --- a/Cantera/src/kinetics/InterfaceKinetics.h +++ b/Cantera/src/kinetics/InterfaceKinetics.h @@ -592,7 +592,7 @@ namespace Cantera { * @param iphase Index of the phase. This is the order within the internal thermo vector object * @param exists Boolean indicating whether the phase exists or not */ - void setPhaseExistence(const int iphase, const int exists); + void setPhaseExistence(const size_t iphase, const bool exists); //! Set the stability of a phase in the reaction object /*! diff --git a/Cantera/src/numerics/BandMatrix.cpp b/Cantera/src/numerics/BandMatrix.cpp index 919512281..2d020355d 100644 --- a/Cantera/src/numerics/BandMatrix.cpp +++ b/Cantera/src/numerics/BandMatrix.cpp @@ -31,7 +31,7 @@ namespace Cantera { ludata.clear(); } //==================================================================================================================== - BandMatrix::BandMatrix(size_t n, size_ kl, size_t ku, doublereal v) : + BandMatrix::BandMatrix(size_t n, size_t kl, size_t ku, doublereal v) : GeneralMatrix(1), m_factored(false), m_n(n), diff --git a/Cantera/src/numerics/BandMatrix.h b/Cantera/src/numerics/BandMatrix.h index 975d44de9..a27076d35 100644 --- a/Cantera/src/numerics/BandMatrix.h +++ b/Cantera/src/numerics/BandMatrix.h @@ -148,7 +148,7 @@ namespace Cantera { * * @return Returns the value of the matrix entry */ - doublereal _value(size_t i, size_t j) const; + doublereal _value(int i, int j) const; //! Returns the number of rows virtual size_t nRows() const; diff --git a/Cantera/src/numerics/DAE_Solver.h b/Cantera/src/numerics/DAE_Solver.h new file mode 100644 index 000000000..cbc4d4b03 --- /dev/null +++ b/Cantera/src/numerics/DAE_Solver.h @@ -0,0 +1,250 @@ +/** + * + * @file DAE_Solver.h + * + * Header file for class DAE_Solver + */ + +/* + * $Date$ + * $Revision$ + * + * Copyright 2006 California Institute of Technology + * + */ +#ifndef CT_DAE_Solver_H +#define CT_DAE_Solver_H + +#include + +#include "ct_defs.h" +#include "ResidJacEval.h" +#include "global.h" + +namespace Cantera { + +#define DAE_DEVEL +#ifdef DAE_DEVEL + + class Jacobian { + public: + Jacobian(){} + virtual ~Jacobian(){} + virtual bool supplied() { return false; } + virtual bool isBanded() { return false; } + virtual int lowerBandWidth() { return 0; } + virtual int upperBandWidth() { return 0; } + }; + + class BandedJacobian : public Jacobian { + public: + BandedJacobian(int ml, int mu) { + m_ml = ml; m_mu = mu; + } + virtual bool supplied() { return false; } + virtual bool isBanded() { return true; } + virtual int lowerBandWidth() { return m_ml; } + virtual int upperBandWidth() { return m_mu; } + protected: + int m_ml, m_mu; + }; + + const int cDirect = 0; + const int cKrylov = 1; + + + + /** + * Wrapper for DAE solvers + */ + class DAE_Solver { + public: + + DAE_Solver(ResidJacEval& f) : + m_resid(f), + m_neq(f.nEquations()), + m_time(0.0) + { + } + + virtual ~DAE_Solver(){} + + /** + * Set error tolerances. This version specifies a scalar + * relative tolerance, and a vector absolute tolerance. + */ + virtual void setTolerances(doublereal reltol, + doublereal* abstol) { + warn("setTolerances"); + } + + /** + * Set error tolerances. This version specifies a scalar + * relative tolerance, and a scalar absolute tolerance. + */ + virtual void setTolerances(doublereal reltol, doublereal abstol) { + warn("setTolerances"); + } + + /** + * Specify a Jacobian evaluator. If this method is not called, + * the Jacobian will be computed by finite difference. + */ + void setJacobian(Jacobian& jac) { + warn("setJacobian"); + } + + virtual void setLinearSolverType(int solverType) { + warn("setLinearSolverType"); + } + + virtual void setDenseLinearSolver() { + warn("setDenseLinearSolver"); + } + + virtual void setBandedLinearSolver(int m_upper, int m_lower) { + warn("setBandedLinearSolver"); + } + virtual void setMaxStepSize(doublereal dtmax) { + warn("setMaxStepSize"); + } + virtual void setMaxOrder(int n) { + warn("setMaxOrder"); + } + virtual void setMaxNumSteps(int n) { + warn("setMaxNumSteps"); + } + virtual void setInitialStepSize(doublereal h0) { + warn("setInitialStepSize"); + } + virtual void setStopTime(doublereal tstop) { + warn("setStopTime"); + } + virtual void setMaxErrTestFailures(int n) { + warn("setMaxErrTestFailures"); + } + virtual void setMaxNonlinIterations(int n) { + warn("setMaxNonlinIterations"); + } + virtual void setMaxNonlinConvFailures(int n) { + warn("setMaxNonlinConvFailures"); + } + virtual void inclAlgebraicInErrorTest(bool yesno) { + warn("inclAlgebraicInErrorTest"); + } + + /** + * This method may be called if the initial conditions do not + * satisfy the residual equation F = 0. Given the derivatives + * of all variables, this method computes the initial y + * values. + */ + virtual void correctInitial_Y_given_Yp(doublereal* y, doublereal* yp, + doublereal tout) { + warn("correctInitial_Y_given_Yp"); + } + + /** + * This method may be called if the initial conditions do not + * satisfy the residual equation F = 0. Given the initial + * values of all differential variables, it computes the + * initial values of all algebraic variables and the initial + * derivatives of all differential variables. + */ + virtual void correctInitial_YaYp_given_Yd(doublereal* y, doublereal* yp, + doublereal tout) + { + warn("correctInitial_YaYp_given_Yd"); + } + + /** + * Solve the system of equations up to time tout. + */ + virtual int solve(doublereal tout) { + warn("solve"); return 0; + } + + /** + * Take one internal step. + */ + virtual doublereal step(doublereal tout) { + warn("step"); return 0; + } + + /// Number of equations. + int nEquations() const { + return m_resid.nEquations(); + } + + /** + * initialize. Base class method does nothing. + */ + virtual void init(doublereal t0) {} + + /** + * Set a solver-specific input parameter. + */ + virtual void setInputParameter(int flag, doublereal value) { + warn("setInputParameter"); + } + + /** + * Get the value of a solver-specific output parameter. + */ + virtual doublereal getOutputParameter(int flag) const { + warn("getOutputParameter"); return 0.0; + } + + /// the current value of solution component k. + virtual doublereal solution(int k) const { + warn("solution"); return 0.0; + } + + virtual const doublereal* solutionVector() const { + warn("solutionVector"); return &m_dummy; + } + + /// the current value of the derivative of solution component k. + virtual doublereal derivative(int k) const { + warn("derivative"); return 0.0; + } + + virtual const doublereal* derivativeVector() const { + warn("derivativeVector"); return &m_dummy; + } + + protected: + + doublereal m_dummy; + + ResidJacEval& m_resid; + + //! Number of total equations in the system + integer m_neq; + doublereal m_time; + + + private: + void warn(std::string msg) const { + writelog(">>>> Warning: method "+msg+" of base class " + +"DAE_Solver called. Nothing done.\n"); + } + }; + + + //! Factor method for choosing a DAE solver + /*! + * + * @param itype String identifying the type + * (IDA is the only option) + * @param f Residual function to be solved by the DAE algorithm + * + * @return Returns a point to the instantiated DAE_Solver object + */ + DAE_Solver* newDAE_Solver(std::string itype, ResidJacEval& f); + +#endif + +} + +#endif diff --git a/Cantera/src/numerics/FuncEval.h b/Cantera/src/numerics/FuncEval.h index 385a6a7e9..a7a0132a8 100644 --- a/Cantera/src/numerics/FuncEval.h +++ b/Cantera/src/numerics/FuncEval.h @@ -51,7 +51,7 @@ namespace Cantera { virtual size_t neq()=0; //! Number of parameters. - virtual int nparams() { return 0; } + virtual size_t nparams() { return 0; } protected: diff --git a/Cantera/src/numerics/RootFind.cpp b/Cantera/src/numerics/RootFind.cpp index 4d5760fca..a84153f55 100644 --- a/Cantera/src/numerics/RootFind.cpp +++ b/Cantera/src/numerics/RootFind.cpp @@ -14,6 +14,11 @@ #include "ct_defs.h" #include "RootFind.h" +// turn on debugging for now +#ifndef DEBUG_MODE +#define DEBUG_MODE +#endif + #include "global.h" #ifdef DEBUG_MODE #include "mdp_allo.h" @@ -53,11 +58,6 @@ namespace Cantera { #endif #ifndef SWAP #define SWAP(x1, x2, tmp) ((tmp) = (x2), (x2) = (x1), (x1) = (tmp)) -#endif - -// turn on debugging for now -#ifndef DEBUG_MODE -#define DEBUG_MODE #endif /*****************************************************************************/ diff --git a/Cantera/src/oneD/Solid1D.cpp b/Cantera/src/oneD/Solid1D.cpp deleted file mode 100644 index 4f6be4b61..000000000 --- a/Cantera/src/oneD/Solid1D.cpp +++ /dev/null @@ -1,580 +0,0 @@ -/** - * @file Solid1D.cpp - */ - -/* - * $Author$ - * $Revision$ - * $Date$ - */ - -// Copyright 2003 California Institute of Technology - - -// turn off warnings under Windows -#ifdef WIN32 -#pragma warning(disable:4786) -#pragma warning(disable:4503) -#endif - -#include -#include - -#include "Solid1D.h" -#include "../ArrayViewer.h" -#include "../ctml.h" -#include "MultiJac.h" - -using namespace ctml; - -namespace Cantera { - - - int Solid1D::c_T_loc = 0; - int Solid1D::c_C_loc = 1; - - Solid1D::Solid1D(ThermoPhase* ph, int points) : - Domain1D(1, points), - m_kin(0), - m_trans(0), - m_jac(0), - m_ok(false) - { - m_type = cSolidType; - m_points = points; - m_thermo = ph; - - if (ph == 0) { m_nsp = 1; return; }// used to create a dummy object - - m_nsp = m_thermo->nSpecies(); - Domain1D::resize(m_nsp+1, points); - - // make a local copy of the species molecular weight vector - m_wt = m_thermo->molecularWeights(); - - m_nv = m_nsp + 1; - - // turn off the energy equation at all points - m_do_energy.resize(m_points,false); - m_do_species.resize(m_nsp,false); - - m_diff.resize(m_nsp*m_points); - m_flux.resize(m_nsp,m_points); - m_wdot.resize(m_nsp,m_points, 0.0); - m_cbar.resize(m_nsp); - - - //-------------- default solution bounds -------------------- - - vector_fp vmin(m_nv), vmax(m_nv); - - // temperature bounds - vmin[c_T_loc] = 200.0; - vmax[c_T_loc]= 1.e9; - - // concentration bounds - int k; - for (k = 0; k < m_nsp; k++) { - vmin[c_C_loc + k] = -1.0e-5; - vmax[c_C_loc + k] = 1.0e5; - } - setBounds(vmin.size(), vmin.begin(), vmax.size(), vmax.begin()); - - - //-------------------- default error tolerances ---------------- - vector_fp rtol(m_nv, 1.0e-8); - vector_fp atol(m_nv, 1.0e-15); - setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin(),false); - setTolerances(rtol.size(), rtol.begin(), atol.size(), atol.begin(),true); - - //-------------------- grid refinement ------------------------- - m_refiner->setActive(c_T_loc, false); - - vector_fp gr; - for (int ng = 0; ng < m_points; ng++) gr.push_back(1.0*ng/m_points); - setupGrid(m_points, gr.begin()); - setID("solid"); - } - - - /** - * Change the grid size. Called after grid refinement. - */ - void Solid1D::resize(int points) { - Domain1D::resize(m_nv, points); - - m_rho.resize(m_points, 0.0); - m_wtm.resize(m_points, 0.0); - m_cp.resize(m_points, 0.0); - m_tcon.resize(m_points, 0.0); - m_diff.resize(m_nsp*m_points); - - m_flux.resize(m_nsp,m_points); - m_wdot.resize(m_nsp,m_points, 0.0); - - m_do_energy.resize(m_points,false); - m_fixedtemp.resize(m_points); - - m_dz.resize(m_points-1); - m_z.resize(m_points); - } - - - - void Solid1D::setupGrid(int n, const doublereal* z) { - resize(n); - int j; - m_z[0] = z[0]; - for (j = 1; j < m_points; j++) { - m_z[j] = z[j]; - m_dz[j-1] = m_z[j] - m_z[j-1]; - } - } - - - /** - * Install a transport manager. - */ - void Solid1D::setTransport(Transport& trans) { - m_trans = &trans; - - if (m_trans->model() != cSolidTransport) { - throw CanteraError("setTransport","unknown transport model."); - } - - - /** - * Set the solid object state to be consistent with the solution at - * point j. - */ - void Solid1D::setThermoState(const doublereal* x,int j) { - m_thermo->setTemperature(T(x,j)); - const doublereal* yy = x + m_nv*j + 1; - m_thermo->setConcentrations(yy); - } - - - /** - * Set the state to be consistent with the solution at the - * midpoint between j and j + 1. - */ - void Solid1D::setStateAtMidpoint(const doublereal* x,int j) { - m_thermo->setTemperature(0.5*(T(x,j)+T(x,j+1))); - const doublereal* ccj = x + m_nv*j + 1; - const doublereal* ccjp = x + m_nv*(j+1) + 1; - for (int k = 0; k < m_nsp; k++) - m_ybar[k] = 0.5*(ccj[k] + ccjp[k]); - m_thermo->setConcentrations(m_cbar.begin()); - } - - - - void Solid1D::eval(int jg, doublereal* xg, - doublereal* rg, integer* diagg, doublereal rdt) { - - // if evaluating a Jacobian, and the global point is outside - // the domain of influence for this domain, then skip - // evaluating the residual - if (jg >=0 && (jg < firstPoint() - 1 || jg > lastPoint() + 1)) return; - - // if evaluating a Jacobian, compute the steady-state residual - if (jg >= 0) rdt = 0.0; - - // start of local part of global arrays - doublereal* x = xg + loc(); - doublereal* rsd = rg + loc(); - integer* diag = diagg + loc(); - - int jmin, jmax, jpt; - jpt = jg - firstPoint(); - - if (jg < 0) { // evaluate all points - jmin = 0; - jmax = m_points - 1; - } - else { // evaluate points for Jacobian - jmin = max(jpt-1, 0); - jmax = min(jpt+1,m_points-1); - } - - // properties are computed for grid points from j0 to j1 - int j0 = max(jmin-1,0); - int j1 = min(jmax+1,m_points-1); - - - int j, k; - - //----------------------------------------------------- - // update properties - //----------------------------------------------------- - - // thermodynamic properties only if a Jacobian is - // not being evaluated - if (jpt < 0) - updateThermo(x, j0, j1); - - // update transport properties only if a Jacobian is - // not being evaluated - if (jpt < 0) - updateTransport(x, j0, j1); - - // update the species diffusive mass fluxes whether or not a - // Jacobian is being evaluated - updateDiffFluxes(x, j0, j1); - - for (j = j0; j <= j1; j++) { - setThermoState(j); - } - - //---------------------------------------------------- - // evaluate the residual equations at all required - // grid points - //---------------------------------------------------- - - for (j = jmin; j <= jmax; j++) { - - - //---------------------------------------------- - // left boundary - //---------------------------------------------- - - if (j == 0) { - rsd[index(c_T_loc,0)] = T(x,0); - - // The default boundary condition for species is zero - // flux. However, the boundary object may modify - // this. - for (k = 0; k < m_nsp; k++) { - rsd[index(c_C_loc + k, 0)] = - m_flux(k,0); - } - } - - - //---------------------------------------------- - // - // right boundary - // - //---------------------------------------------- - - else if (j == m_points - 1) { - rsd[index(c_T_loc,j)] = T(x,j); - for (k = 0; k < m_nsp; k++) { - rsd[index(k+c_C_loc,j)] = m_flux(k,j-1); - } - } - - - //------------------------------------------ - // interior points - //------------------------------------------ - - else { - - //------------------------------------------------- - // Species equations - // - // \rho u dY_k/dz + dJ_k/dz + M_k\omega_k - // - //------------------------------------------------- - getWdot(x,j); - - doublereal diffus; - for (k = 0; k < m_nsp; k++) { - diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1)) - /(z(j+1) - z(j-1)); - rsd[index(c_C_loc + k, j)] - = wdot(k,j) - diffus - - rdt*(C(x,k,j) - C_prev(k,j)); - diag[index(c_C_loc + k, j)] = 1; - } - - //----------------------------------------------- - // energy equation - //----------------------------------------------- - - if (m_do_energy[j]) { - - rsd[index(c_T_loc, j)] = - divHeatFlux(x,j); - rsd[index(c_T_loc, j)] /= (m_rho[j]*m_cp[j]); - - rsd[index(c_T_loc, j)] -= rdt*(T(x,j) - T_prev(j)); - diag[index(c_T_loc, j)] = 1; - } - } - - // residual equations if the energy or species equations - // are disabled - - if (!m_do_energy[j]) { - rsd[index(c_T_loc, j)] = T(x,j) - T_fixed(j); - diag[index(c_T_loc, j)] = 0; - } - } - } - - - - /** - * Update the transport properties at grid points in the range - * from j0 to j1, based on solution x. - */ - void Surf1D::updateTransport(doublereal* x,int j0, int j1) { - int j; - for (j = j0; j < j1; j++) { - setStateAtMidpoint(x,j); - m_trans->getMixDiffCoeffs(m_diff.begin() + j*m_nsp); - m_tcon[j] = m_trans->thermalConductivity(); - } - } - - - /** - * Print the solution. - */ - void Solid1D::showSolution(const doublereal* x) { - int nn = m_nv/5; - int i, j, n; - char* buf = new char[100]; - - // The mean molecular weight is needed to convert - updateThermo(x, 0, m_points-1); - - for (i = 0; i < nn; i++) { - drawline(); - sprintf(buf, "\n z "); - writelog(buf); - for (n = 0; n < 5; n++) { - sprintf(buf, " %10s ",componentName(i*5 + n).c_str()); - writelog(buf); - } - drawline(); - for (j = 0; j < m_points; j++) { - sprintf(buf, "\n %10.4g ",m_z[j]); - writelog(buf); - for (n = 0; n < 5; n++) { - sprintf(buf, " %10.4g ",component(x, i*5+n,j)); - writelog(buf); - } - } - writelog("\n"); - } - int nrem = m_nv - 5*nn; - drawline(); - sprintf(buf, "\n z "); - writelog(buf); - for (n = 0; n < nrem; n++) { - sprintf(buf, " %10s ", componentName(nn*5 + n).c_str()); - writelog(buf); - } - drawline(); - for (j = 0; j < m_points; j++) { - sprintf(buf, "\n %10.4g ",m_z[j]); - writelog(buf); - for (n = 0; n < nrem; n++) { - sprintf(buf, " %10.4g ",component(x, nn*5+n,j)); - writelog(buf); - } - } - writelog("\n"); - } - - - /** - * Update the diffusive mass fluxes. - */ - void Solid1D::updateDiffFluxes(const doublereal* x, int j0, int j1) { - int j, k, m; - doublereal sum, wtm, rho, dz, gradlogT, s; - doublereal dphidz, a1; - for (j = j0; j < j1; j++) { - sum = 0.0; - rho = density(j); - dz = z(j+1) - z(j); - for (k = 0; k < m_nsp; k++) { - m_flux(k,j) = m_diff[k+m_nsp*j] * - (C(x,k,j) - C(x,k,j+1))/dz; - sum -= m_flux(k,j); - } - for (k = 0; k < m_nsp; k++) m_flux(k,j) += C(x,k,j)*sum; - } - break; - } - - - void Solid1D::outputTEC(ostream &s, const doublereal* x, - string title, int zone) { - int j,k; - s << "TITLE = \"" + title + "\"" << endl; - s << "VARIABLES = \"Z (m)\"" << endl; - s << "\"T (K)\"" << endl; - - for (k = 0; k < m_nsp; k++) { - s << "\"" << m_thermo->speciesName(k) << "\"" << endl; - } - s << "ZONE T=\"c" << zone << "\"" << endl; - s << " I=" << m_points << ",J=1,K=1,F=POINT" << endl; - s << "DT=(SINGLE"; - for (k = 0; k < m_nsp; k++) s << " SINGLE"; - s << " )" << endl; - for (j = 0; j < m_points; j++) { - s << z(j) << " "; - for (k = 0; k < m_nv; k++) { - s << component(x, k, j) << " "; - } - s << endl; - } - } - - - string Solid1D::componentName(int n) const { - switch(n) { - case c_T_loc: return "T"; - default: - if (n >= (int) 1 && n < (int) (c_C_loc + m_nsp)) { - return m_thermo->speciesName(n - 1); - } - else - return ""; - } - } - - - void Solid1D::restore(XML_Node& dom, doublereal* soln) { - - vector ignored; - int nsp = m_thermo->nSpecies(); - vector_int did_species(nsp, 0); - - vector str; - dom.getChildren("string",str); - int nstr = str.size(); - for (int istr = 0; istr < nstr; istr++) { - XML_Node& nd = *str[istr]; - writelog(nd["title"]+": "+nd.value()+"\n"); - } - - map params; - getFloats(dom, params); - - vector d; - dom.child("grid_data").getChildren("floatArray",d); - int nd = d.size(); - - vector_fp x; - int n, np, j, ks, k; - string nm; - bool readgrid = false, wrote_header = false; - for (n = 0; n < nd; n++) { - XML_Node& fa = *d[n]; - nm = fa["title"]; - if (nm == "z") { - getFloatArray(fa,x,false); - np = x.size(); - writelog("Grid contains "+int2str(np)+ - " points.\n"); - readgrid = true; - - // note that setupGrid also resizes the domain. - setupGrid(np, x.begin()); - } - } - if (!readgrid) { - throw CanteraError("Solid1D::restore", - "domain contains no grid points."); - } - - writelog("Importing datasets:\n"); - for (n = 0; n < nd; n++) { - XML_Node& fa = *d[n]; - nm = fa["title"]; - getFloatArray(fa,x,false); - if (nm == "z") { - ; // already read grid - } - else if (nm == "T") { - writelog("temperature "); - if ((int) x.size() == np) { - for (j = 0; j < np; j++) - soln[index(c_T_loc,j)] = x[j]; - - // For fixed-temperature simulations, use the imported temperature profile by default. - // If this is not desired, call setFixedTempProfile *after* restoring the solution. - vector_fp zz(np); - for (int jj = 0; jj < np; jj++) zz[jj] = (grid(jj) - zmin())/(zmax() - zmin()); - setFixedTempProfile(zz, x); - } - else goto error; - } - else if (m_thermo->speciesIndex(nm) >= 0) { - writelog(nm+" "); - if ((int) x.size() == np) { - k = m_thermo->speciesIndex(nm); - did_species[k] = 1; - for (j = 0; j < np; j++) - soln[index(k+c_C_loc,j)] = x[j]; - } - } - else - ignored.push_back(nm); - } - - if (ignored.size() != 0) { - writelog("\n\n"); - writelog("Ignoring datasets:\n"); - int nn = ignored.size(); - for (int n = 0; n < nn; n++) { - writelog(ignored[n]+" "); - } - } - - for (ks = 0; ks < nsp; ks++) { - if (did_species[ks] == 0) { - if (!wrote_header) { - writelog("Missing data for species:\n"); - wrote_header = true; - } - writelog(m_thermo->speciesName(ks)+" "); - } - } - - return; - error: - throw CanteraError("Solid1D::restore","Data size error"); - } - - - - void Solid1D::save(XML_Node& o, const doublereal * const sol) { - int k; - - ArrayViewer soln(m_nv, m_points, sol + loc()); - - XML_Node& flow = (XML_Node&)o.addChild("domain"); - flow.addAttribute("type",flowType()); - flow.addAttribute("id",m_id); - flow.addAttribute("points",m_points); - flow.addAttribute("components",m_nv); - - if (m_desc != "") addString(flow,"description",m_desc); - XML_Node& gv = flow.addChild("grid_data"); - addFloatArray(gv,"z",m_z.size(),m_z.begin(), - "m","length"); - vector_fp x(soln.nColumns()); - - soln.getRow(c_T_loc,x.begin()); - addFloatArray(gv,"T",x.size(),x.begin(),"K","temperature",0.0); - - for (k = 0; k < m_nsp; k++) { - soln.getRow(c_C_loc+k,x.begin()); - addFloatArray(gv,m_thermo->speciesName(k), - x.size(),x.begin(),"","concentration",0.0,1.0); - } - } - - - void Solid1D::setJac(MultiJac* jac) { - m_jac = jac; - } - - -} diff --git a/Cantera/src/oneD/Solid1D.h b/Cantera/src/oneD/Solid1D.h deleted file mode 100644 index c0f32964f..000000000 --- a/Cantera/src/oneD/Solid1D.h +++ /dev/null @@ -1,410 +0,0 @@ -/** - * @file Solid1D.h - * - */ - -/* - * $Author$ - * $Revision$ - * $Date$ - */ - -// Copyright 2001 California Institute of Technology - -#ifndef CT_SOLID1D_H -#define CT_SOLID1D_H - -#include "../transport/TransportBase.h" -#include "Domain1D.h" -#include "../Array.h" -#include "../sort.h" -#include "../ThermoPhase.h" -#include "../Kinetics.h" -#include "../funcs.h" - - -namespace Cantera { - - class MultiJac; - - - //----------------------------------------------------------- - // Class Solid1D - //----------------------------------------------------------- - - - /** - * A class for one-dimensional reacting solids with current - * transport. This class implements the one-dimensional - * similarity solution for a chemically-reacting, axisymmetric, - * stagnation-point flow. - */ - class Solid1D : public Domain1D { - - public: - - - //------------------------------------------ - // constants - //------------------------------------------ - - /** - * Offsets of solution components in the solution array. - */ - const unsigned int c_phi_loc; // electric potential - const unsigned int c_T_loc; // temperature - const unsigned int c_C_loc; // concentrations - - - //-------------------------------- - // construction and destruction - //-------------------------------- - - // Constructor. - Solid1D(ThermoPhase* ph = 0, int nsp = 1, int points = 1); - - /// Destructor. - virtual ~Solid1D(){} - - - /** - * @name Problem Specification - */ - //@{ - - virtual void setupGrid(int n, const doublereal* z); - - thermo_t& phase() { return *m_thermo; } - kinetics_t& kinetics() { return *m_kin; } - - /** - * Set the thermo manager. - */ - void setThermo(thermo_t& th) { - m_thermo = &th; - } - - /// set the kinetics manager - void setKinetics(kinetics_t& kin) { m_kin = &kin; } - - /// set the transport manager - void setTransport(Transport& trans); - - virtual void setState(int point, const doublereal* state) { - setTemperature(point, state[c_T_loc]); - setElectricPotential(point, state[c_phi_loc]); - int k; - for (k = 0; k < m_nsp; k++) { - setConcentration(point, k, state[c_C_loc+k]); - } - } - - - virtual void _getInitialSoln(doublereal* x) { - int k, j; - for (j = 0; j < m_points; j++) { - x[index(c_T_loc,j)] = T_fixed(j); - x[index(c_phi_loc,j)] = phi_fixed(j); - for (k = 0; k < m_nsp; k++) { - x[index(c_C_loc+k,j)] = C_fixed(k,j); - } - } - } - - virtual void _finalize(const doublereal* x) { - int k, j; - doublereal zz, tt; - int nz = m_zfix.size(); - bool e = m_do_energy[0]; - for (j = 0; j < m_points; j++) { - if (e || nz == 0) - setTemperature(j, T(x, j)); - else { - zz = (z(j) - z(0))/(z(m_points - 1) - z(0)); - tt = linearInterp(zz, m_zfix, m_tfix); - setTemperature(j, tt); - } - setElectricPotential(j, phi(x,j)); - for (k = 0; k < m_nsp; k++) { - setConcentration(j, k, C(x, k, j)); - } - } - if (e) solveEnergyEqn(); - } - - - void setFixedTempProfile(vector_fp& zfixed, vector_fp& tfixed) { - m_zfix = zfixed; - m_tfix = tfixed; - } - - /** - * Set the temperature fixed point at grid point j, and - * disable the energy equation so that the solution will be - * held to this value. - */ - void setTemperature(int j, doublereal t) { - m_fixedtemp[j] = t; - m_do_energy[j] = false; - } - - /** - * Set the electric potential fixed point at grid point j, and - * disable Gauss's equation so that the solution will be - * held to this value. - */ - void setElectricPotential(int j, doublereal phi) { - m_fixedphi[j] = phi; - m_do_gauss[j] = false; - } - - /** - * Set the mass fraction fixed point for species k at grid - * point j, and disable the species equation so that the - * solution will be held to this value. - */ - void setConcentration(int j, int k, doublereal c) { - m_fixedc(k,j) = c; - m_do_species[k] = true; // false; - } - - /** - * The fixed temperature value at point j. - */ - doublereal T_fixed(int j) const {return m_fixedtemp[j];} - - /** - * The fixed potential value at point j. - */ - doublereal phi_fixed(int j) const {return m_fixedphi[j];} - - /** - * The fixed mass fraction value of species k at point j. - */ - doublereal C_fixed(int k, int j) const {return m_fixedc(k,j);} - - virtual std::string componentName(int n) const; - - void setDielectricConstant(doublereal e) { m_eps = e; } - doublereal dielectricConstant() { return m_eps; } - - /** - * Write a Tecplot zone corresponding to the current solution. - * May be called multiple times to generate animation. - */ - void outputTEC(ostream &s, const doublereal* x, - std::string title, int zone); - - virtual void showSolution(const doublereal* x); - - //! Save the current solution for this domain into an XML_Node - /*! - * - * @param o XML_Node to save the solution to. - * @param sol Current value of the solution vector. - * The object will pick out which part of the solution - * vector pertains to this object. - */ - virtual void save(XML_Node& o, const doublereal * const sol); - - virtual void restore(XML_Node& dom, doublereal* soln); - - // overloaded in subclasses - virtual std::string solidType() { return ""; } - - void solveEnergyEqn(int j=-1) { - if (j < 0) - for (int i = 0; i < m_points; i++) - m_do_energy[i] = true; - else - m_do_energy[j] = true; - m_refiner->setActive(c_T_loc, true); - needJacUpdate(); - } - - void fixTemperature(int j=-1) { - if (j < 0) - for (int i = 0; i < m_points; i++) { - m_do_energy[i] = false; - } - else m_do_energy[j] = false; - m_refiner->setActive(c_T_loc, false); - needJacUpdate(); - } - - void solveGaussEqn(int j=-1) { - if (j < 0) - for (int i = 0; i < m_points; i++) - m_do_gauss[i] = true; - else - m_do_gauss[j] = true; - m_refiner->setActive(c_phi_loc, true); - needJacUpdate(); - } - - void fixElectricPotential(int j=-1) { - if (j < 0) - for (int i = 0; i < m_points; i++) { - m_do_gauss[i] = false; - } - else m_do_gauss[j] = false; - m_refiner->setActive(c_phi_loc, false); - needJacUpdate(); - } - - bool doSpecies(int k) { return m_do_species[k]; } - bool doEnergy(int j) { return m_do_energy[j]; } - bool doGauss(int j) { return m_do_gauss[j]; } - - void solveSpecies(int k=-1) { - if (k == -1) { - for (int i = 0; i < m_nsp; i++) - m_do_species[i] = true; - } - else m_do_species[k] = true; - needJacUpdate(); - } - - void fixSpecies(int k=-1) { - if (k == -1) { - for (int i = 0; i < m_nsp; i++) - m_do_species[i] = false; - } - else m_do_species[k] = false; - needJacUpdate(); - } - - void resize(int points); - - void setJac(MultiJac* jac); - void setThermoState(const doublereal* x,int j); - void setStateAtMidpoint(const doublereal* x,int j); - - - protected: - - doublereal component(const doublereal* x, int i, int j) const { - doublereal xx = x[index(i,j)]; - return xx; - } - - doublereal wdot(int k, int j) const {return m_wdot(k,j);} - - /// write the net production rates at point j into array m_wdot - void getWdot(doublereal* x,int j) { - setThermoState(x,j); - m_kin->getNetProductionRates(&m_wdot(0,j)); - } - - /** - * update the thermodynamic properties from point - * j0 to point j1 (inclusive), based on solution x. - */ - void updateThermo(const doublereal* x, int j0, int j1) { - int j; - for (j = j0; j <= j1; j++) { - setThermoState(x,j); - m_cp[j] = m_thermo->cp_mass(); - } - } - - - //-------------------------------- - // solution components - //-------------------------------- - - - doublereal T(const doublereal* x,int j) const { - return x[index(c_T_loc, j)]; - } - - doublereal& T(doublereal* x,int j) {return x[index(c_T_loc, j)];} - - doublereal T_prev(int j) const {return prevSoln(c_T_loc, j);} - - doublereal C(const doublereal* x,int k, int j) const { - return x[index(c_C_loc + k, j)]; - } - - doublereal& C(doublereal* x,int k, int j) { - return x[index(c_C_loc + k, j)]; - } - - doublereal C_prev(int k, int j) const { - return prevSoln(c_C_loc + k, j); - } - - doublereal flux(int k, int j) const { - return m_flux(k, j); - } - - doublereal phi(const doublereal* x, int j) { - return x[index(c_phi_loc, j)]; - } - - doublereal divHeatFlux(const doublereal* x, int j) const { - doublereal c1 = m_tcon[j-1]*(T(x,j) - T(x,j-1)); - doublereal c2 = m_tcon[j]*(T(x,j+1) - T(x,j)); - return -2.0*(c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); - } - - doublereal divDisplCurr(const doublereal* x, int j) const { - doublereal c1 = (phi(x,j) - phi(x,j-1)); - doublereal c2 = (phi(x,j+1) - phi(x,j)); - return -2.0*m_eps*epsilon_0* - (c2/(z(j+1) - z(j)) - c1/(z(j) - z(j-1)))/(z(j+1) - z(j-1)); - } - - void updateDiffFluxes(const doublereal* x, int j0, int j1); - - //--------------------------------------------------------- - // - // member data - // - //--------------------------------------------------------- - - doublereal m_eps; // relative dielectric constant - - // grid parameters - vector_fp m_dz; - - // mixture thermo properties - vector_fp m_cdens; - - // transport properties - vector_fp m_tcon; - vector_fp m_diff; - Array2D m_flux; - - // production rates - Array2D m_wdot; - - int m_nsp; - - thermo_t* m_thermo; - kinetics_t* m_kin; - Transport* m_trans; - - MultiJac* m_jac; - - bool m_ok; - - // flags - std::vector m_do_energy; - std::vector m_do_species; - std::vector m_do_gauss; - - // fixed T and Y values - Array2D m_fixedy; - Array2D m_fixedphi; - vector_fp m_fixedtemp; - vector_fp m_zfix; - vector_fp m_tfix; - - private: - - vector_fp m_cbar; - }; -} - -#endif diff --git a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp index 1f7d5b4b6..93bf88fb6 100644 --- a/Cantera/src/thermo/GibbsExcessVPSSTP.cpp +++ b/Cantera/src/thermo/GibbsExcessVPSSTP.cpp @@ -212,11 +212,11 @@ namespace Cantera { } - doublereal GibbsExcessVPSSTP::standardConcentration(int k) const { + doublereal GibbsExcessVPSSTP::standardConcentration(size_t k) const { return 1.0; } - doublereal GibbsExcessVPSSTP::logStandardConc(int k) const { + doublereal GibbsExcessVPSSTP::logStandardConc(size_t k) const { return 0.0; } diff --git a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp index 526ad84f3..f5c32b3ea 100644 --- a/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp +++ b/Cantera/src/thermo/IonsFromNeutralVPSSTP.cpp @@ -324,7 +324,7 @@ namespace Cantera { */ //=========================================================================================================== void IonsFromNeutralVPSSTP::getDissociationCoeffs(vector_fp& coeffs, - vector_fp& charges, std::vector& neutMolIndex) const { + vector_fp& charges, std::vector& neutMolIndex) const { coeffs = fm_neutralMolec_ions_; charges = m_speciesCharge; neutMolIndex = fm_invert_ionForNeutral; @@ -1519,7 +1519,7 @@ namespace Cantera { */ GibbsExcessVPSSTP *geThermo = dynamic_cast(neutralMoleculePhase_); if (!geThermo) { - fvo_zero_dbl_1(dlnActCoeffdlnX_diag_, m_kk); + dlnActCoeffdlnX_diag_.assign(m_kk, 0.0); return; } @@ -1577,7 +1577,7 @@ namespace Cantera { */ GibbsExcessVPSSTP *geThermo = dynamic_cast(neutralMoleculePhase_); if (!geThermo) { - fvo_zero_dbl_1(dlnActCoeffdlnN_diag_, m_kk); + dlnActCoeffdlnN_diag_.assign(m_kk, 0.0); return; } diff --git a/Cantera/src/thermo/LatticeSolidPhase.h b/Cantera/src/thermo/LatticeSolidPhase.h index 952aeb53c..c6caf93b6 100644 --- a/Cantera/src/thermo/LatticeSolidPhase.h +++ b/Cantera/src/thermo/LatticeSolidPhase.h @@ -549,13 +549,13 @@ namespace Cantera { * Returns the standard concentration. The units are by definition * dependent on the ThermoPhase and kinetics manager representation. */ - virtual doublereal standardConcentration(int k=0) const; + virtual doublereal standardConcentration(size_t k=0) const; //! Natural logarithm of the standard concentration of the kth species. /*! * @param k index of the species (defaults to zero) */ - virtual doublereal logStandardConc(int k=0) const; + virtual doublereal logStandardConc(size_t k=0) const; //@} /// @name Thermodynamic Values for the Species Reference States -------------------- //@{ diff --git a/Cantera/src/thermo/MargulesVPSSTP.cpp b/Cantera/src/thermo/MargulesVPSSTP.cpp index 1c3053a4b..ca9dcd2ca 100644 --- a/Cantera/src/thermo/MargulesVPSSTP.cpp +++ b/Cantera/src/thermo/MargulesVPSSTP.cpp @@ -715,7 +715,6 @@ namespace Cantera { double RT = GasConstant*T; lnActCoeff_Scaled_.assign(m_kk, 0.0); for (iK = 0; iK < m_kk; iK++) { - XK = moleFractions_[iK]; for (size_t i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; iB = m_pSpecies_B_ij[i]; @@ -744,8 +743,8 @@ namespace Cantera { doublereal XA, XB, g0, g1; doublereal T = temperature(); doublereal RTT = GasConstant*T*T; - fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); - fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); + dlnActCoeffdT_Scaled_.assign(m_kk, 0.0); + d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0); for (iK = 0; iK < m_kk; iK++) { for (size_t i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; @@ -848,7 +847,7 @@ namespace Cantera { double T = temperature(); double RT = GasConstant*T; - dlnActCoeffdlnN_Scaled_.assign(m_kk, 0.0); + dlnActCoeffdlnN_diag_.assign(m_kk, 0.0); for ( iK = 0; iK < m_kk; iK++ ){ @@ -969,7 +968,7 @@ namespace Cantera { doublereal XA, XB, g0 , g1; doublereal T = temperature(); - fvo_zero_dbl_1(dlnActCoeffdlnX_diag_, m_kk); + dlnActCoeffdlnX_diag_.assign(m_kk, 0.0); doublereal RT = GasConstant * T; diff --git a/Cantera/src/thermo/MixedSolventElectrolyte.cpp b/Cantera/src/thermo/MixedSolventElectrolyte.cpp index 922a40bfa..5ee2d55db 100644 --- a/Cantera/src/thermo/MixedSolventElectrolyte.cpp +++ b/Cantera/src/thermo/MixedSolventElectrolyte.cpp @@ -721,7 +721,7 @@ namespace Cantera { double XA, XB, XK, g0 , g1; double T = temperature(); double RT = GasConstant*T; - fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + lnActCoeff_Scaled_.assign(m_kk, 0.0); for (iK = 0; iK < m_kk; iK++) { XK = moleFractions_[iK]; for (int i = 0; i < numBinaryInteractions_; i++) { @@ -752,8 +752,8 @@ namespace Cantera { doublereal XA, XB, g0, g1; doublereal T = temperature(); doublereal RTT = GasConstant*T*T; - fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); - fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); + dlnActCoeffdT_Scaled_.assign(m_kk, 0.0); + d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0); for (iK = 0; iK < m_kk; iK++) { for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; @@ -856,7 +856,7 @@ namespace Cantera { double T = temperature(); double RT = GasConstant*T; - fvo_zero_dbl_1(dlnActCoeffdlnN_diag_, m_kk); + dlnActCoeffdlnN_diag_.assign(m_kk, 0); for ( iK = 0; iK < m_kk; iK++ ){ @@ -977,7 +977,7 @@ namespace Cantera { doublereal XA, XB, g0 , g1; doublereal T = temperature(); - fvo_zero_dbl_1(dlnActCoeffdlnX_diag_, m_kk); + dlnActCoeffdlnX_diag_.assign(m_kk, 0); doublereal RT = GasConstant * T; diff --git a/Cantera/src/thermo/PhaseCombo_Interaction.cpp b/Cantera/src/thermo/PhaseCombo_Interaction.cpp index ef01eb6d7..f82212e8d 100644 --- a/Cantera/src/thermo/PhaseCombo_Interaction.cpp +++ b/Cantera/src/thermo/PhaseCombo_Interaction.cpp @@ -731,7 +731,7 @@ namespace Cantera { doublereal xx; doublereal T = temperature(); doublereal RT = GasConstant*T; - fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + lnActCoeff_Scaled_.assign(m_kk, 0.0); for (iK = 0; iK < m_kk; iK++) { /* @@ -776,8 +776,8 @@ namespace Cantera { doublereal XA, XB, g0, g1; doublereal T = temperature(); doublereal RTT = GasConstant*T*T; - fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); - fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); + dlnActCoeffdT_Scaled_.assign(m_kk, 0.0); + d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0); for (iK = 0; iK < m_kk; iK++) { for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; @@ -906,7 +906,7 @@ namespace Cantera { doublereal RT = GasConstant*T; doublereal xx; - fvo_zero_dbl_1(dlnActCoeffdlnN_diag_, m_kk); + dlnActCoeffdlnN_diag_.assign(m_kk, 0.0); for (iK = 0; iK < m_kk; iK++) { @@ -1023,7 +1023,7 @@ namespace Cantera { doublereal XA, XB, g0 , g1; doublereal T = temperature(); - fvo_zero_dbl_1(dlnActCoeffdlnX_diag_, m_kk); + dlnActCoeffdlnX_diag_.assign(m_kk, 0.0); doublereal RT = GasConstant * T; diff --git a/Cantera/src/thermo/PureFluidPhase.cpp b/Cantera/src/thermo/PureFluidPhase.cpp index 79f1a344b..b90995619 100644 --- a/Cantera/src/thermo/PureFluidPhase.cpp +++ b/Cantera/src/thermo/PureFluidPhase.cpp @@ -659,61 +659,61 @@ namespace Cantera { std::vector cpbar(kk, 0.0); std::vector vbar(kk, 0.0); std::vector pNames; - std::vector data; + std::vector > data; - getMoleFractions(x); + getMoleFractions(&x[0]); pNames.push_back("X"); data.push_back(x); try{ - getMassFractions(y); + getMassFractions(&y[0]); pNames.push_back("Y"); data.push_back(y); } catch (CanteraError) {;} try{ - getChemPotentials(mu); + getChemPotentials(&mu[0]); pNames.push_back("Chem. Pot (J/kmol)"); data.push_back(mu); } catch (CanteraError) {;} try{ - getActivities(a); + getActivities(&a[0]); pNames.push_back("Activity"); data.push_back(a); } catch (CanteraError) {;} try{ - getActivityCoefficients(ac); + getActivityCoefficients(&ac[0]); pNames.push_back("Act. Coeff."); data.push_back(ac); } catch (CanteraError) {;} try{ - getPartialMolarEnthalpies(hbar); + getPartialMolarEnthalpies(&hbar[0]); pNames.push_back("Part. Mol Enthalpy (J/kmol)"); data.push_back(hbar); } catch (CanteraError) {;} try{ - getPartialMolarEntropies(sbar); + getPartialMolarEntropies(&sbar[0]); pNames.push_back("Part. Mol. Entropy (J/K/kmol)"); data.push_back(sbar); } catch (CanteraError) {;} try{ - getPartialMolarIntEnergies(ubar); + getPartialMolarIntEnergies(&ubar[0]); pNames.push_back("Part. Mol. Energy (J/kmol)"); data.push_back(ubar); } catch (CanteraError) {;} try{ - getPartialMolarCp(cpbar); + getPartialMolarCp(&cpbar[0]); pNames.push_back("Part. Mol. Cp (J/K/kmol"); data.push_back(cpbar); } catch (CanteraError) {;} try{ - getPartialMolarVolumes(vbar); + getPartialMolarVolumes(&vbar[0]); pNames.push_back("Part. Mol. Cv (J/K/kmol)"); data.push_back(vbar); } diff --git a/Cantera/src/thermo/RedlichKisterVPSSTP.cpp b/Cantera/src/thermo/RedlichKisterVPSSTP.cpp index 1f3944cd3..5722f41a0 100644 --- a/Cantera/src/thermo/RedlichKisterVPSSTP.cpp +++ b/Cantera/src/thermo/RedlichKisterVPSSTP.cpp @@ -687,7 +687,7 @@ namespace Cantera { doublereal T = temperature(); doublereal RT = GasConstant * T; - fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + lnActCoeff_Scaled_.assign(m_kk, 0.0); /* * Scaling: I moved the division of RT higher so that we are always dealing with G/RT dimensionless terms @@ -764,8 +764,8 @@ namespace Cantera { doublereal XA, XB; // doublereal T = temperature(); - fvo_zero_dbl_1(dlnActCoeffdT_Scaled_, m_kk); - fvo_zero_dbl_1(d2lnActCoeffdT2_Scaled_, m_kk); + dlnActCoeffdT_Scaled_.assign(m_kk, 0.0); + d2lnActCoeffdT2_Scaled_.assign(m_kk, 0.0); for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; @@ -1056,7 +1056,7 @@ namespace Cantera { doublereal RT = GasConstant * T; double Volts = 0.0; - fvo_zero_dbl_1(lnActCoeff_Scaled_, m_kk); + lnActCoeff_Scaled_.assign(m_kk, 0.0); for (int i = 0; i < numBinaryInteractions_; i++) { iA = m_pSpecies_A_ij[i]; diff --git a/Cantera/src/thermo/ShomateThermo.h b/Cantera/src/thermo/ShomateThermo.h index a2b233952..4dfb47d56 100644 --- a/Cantera/src/thermo/ShomateThermo.h +++ b/Cantera/src/thermo/ShomateThermo.h @@ -159,7 +159,7 @@ namespace Cantera { * @see ShomatePoly * @see ShomatePoly2 */ - virtual void install(std::string name, size_t index, size_t type, + virtual void install(std::string name, size_t index, int type, const doublereal* c, doublereal minTemp, doublereal maxTemp, doublereal refPressure) { diff --git a/Cantera/src/thermo/SimpleThermo.h b/Cantera/src/thermo/SimpleThermo.h index a3000e149..74db5e4c2 100644 --- a/Cantera/src/thermo/SimpleThermo.h +++ b/Cantera/src/thermo/SimpleThermo.h @@ -147,7 +147,7 @@ namespace Cantera { * * @see ConstCpPoly */ - virtual void install(std::string name, size_t index, size_t type, const doublereal* c, + virtual void install(std::string name, size_t index, int type, const doublereal* c, doublereal minTemp, doublereal maxTemp, doublereal refPressure) { m_logt0.push_back(log(c[0])); diff --git a/Cantera/src/transport/DustyGasTransport.h b/Cantera/src/transport/DustyGasTransport.h index 3420d5f5f..43e7fc511 100644 --- a/Cantera/src/transport/DustyGasTransport.h +++ b/Cantera/src/transport/DustyGasTransport.h @@ -130,7 +130,7 @@ namespace Cantera { * d[ld*j+i] is the D_ij diffusion coefficient (the diffusion * coefficient for species i due to species j). */ - virtual void getMultiDiffCoeffs(const int ld, doublereal* const d); + virtual void getMultiDiffCoeffs(const size_t ld, doublereal* const d); //! Get the molar fluxes [kmol/m^2/s], given the thermodynamic state at two nearby points. /*! diff --git a/Cantera/src/transport/L_matrix.h b/Cantera/src/transport/L_matrix.h index 7bd409f52..4dbb3224d 100644 --- a/Cantera/src/transport/L_matrix.h +++ b/Cantera/src/transport/L_matrix.h @@ -133,7 +133,7 @@ namespace Cantera { doublereal constant, sum; size_t n2 = 2*m_nsp; int npoly = 0; - for (j = 0; j < m_nsp; j++) { + for (size_t j = 0; j < m_nsp; j++) { // collect terms that depend only on "j" if (hasInternalModes(j)) { constant = prefactor*m_mw[j]*x[j]*m_crot[j]/(m_cinternal[j]*m_rotrelax[j]); @@ -209,7 +209,7 @@ namespace Cantera { - constant1*sum; } else { - for (k = 0; k < m_nsp; k++) { + for (size_t k = 0; k < m_nsp; k++) { m_Lmatrix(i+n2,i+n2) = 1.0; } } diff --git a/Cantera/src/transport/LiquidTranInteraction.cpp b/Cantera/src/transport/LiquidTranInteraction.cpp index 7d52b2127..8971b63d3 100644 --- a/Cantera/src/transport/LiquidTranInteraction.cpp +++ b/Cantera/src/transport/LiquidTranInteraction.cpp @@ -676,15 +676,15 @@ namespace Cantera { m_thermo->getMoleFractions(molefracs ); vector_fp neut_molefracs; ions_thermo->getNeutralMolecMoleFractions(neut_molefracs); - vector cation; - vector anion; + vector cation; + vector anion; ions_thermo->getCationList(cation); ions_thermo->getAnionList(anion); // Reaction Coeffs and Charges std::vector viS(6); std::vector charges(3); - std::vector neutMolIndex(3); + std::vector neutMolIndex(3); ions_thermo->getDissociationCoeffs(viS,charges,neutMolIndex); if ((int)anion.size() != 1) { diff --git a/Cantera/src/transport/LiquidTransport.cpp b/Cantera/src/transport/LiquidTransport.cpp index 7403e5a08..435912ca0 100644 --- a/Cantera/src/transport/LiquidTransport.cpp +++ b/Cantera/src/transport/LiquidTransport.cpp @@ -1661,7 +1661,7 @@ namespace Cantera { */ condSum1 = 0; - for (i = 0; i < m_nsp; i++){ + for (size_t i = 0; i < m_nsp; i++){ condSum1 -= Faraday*m_chargeSpecies[i]*m_B(i,0)*m_molefracs_tran[i]/vol; } diff --git a/Cantera/src/transport/LiquidTransport.h b/Cantera/src/transport/LiquidTransport.h index 33fcef55f..a77eef9ac 100644 --- a/Cantera/src/transport/LiquidTransport.h +++ b/Cantera/src/transport/LiquidTransport.h @@ -501,7 +501,7 @@ namespace Cantera { * Flat vector with the m_nsp in the inner loop. * length = ldx * ndim */ - virtual void getSpeciesVdiff(int ndim, + virtual void getSpeciesVdiff(size_t ndim, const doublereal* grad_T, int ldx, const doublereal* grad_X, @@ -534,7 +534,7 @@ namespace Cantera { * Flat vector with the m_nsp in the inner loop. * length = ldx * ndim */ - virtual void getSpeciesVdiffES(int ndim, const doublereal* grad_T, + virtual void getSpeciesVdiffES(size_t ndim, const doublereal* grad_T, int ldx, const doublereal* grad_X, int ldf, const doublereal* grad_Phi, doublereal* Vdiff) ; @@ -658,7 +658,7 @@ namespace Cantera { * Flat vector with the m_nsp in the inner loop. * length = ldx * ndim */ - virtual void getSpeciesVdiffExt(int ldf, doublereal* Vdiff); + virtual void getSpeciesVdiffExt(size_t ldf, doublereal* Vdiff); //! Return the species diffusive fluxes relative to //! the averaged velocity. @@ -675,7 +675,7 @@ namespace Cantera { * Flat vector with the m_nsp in the inner loop. * length = ldx * ndim */ - virtual void getSpeciesFluxesExt(int ldf, doublereal* fluxes); + virtual void getSpeciesFluxesExt(size_t ldf, doublereal* fluxes); protected: //! Returns true if temperature has changed, diff --git a/Cantera/src/transport/SimpleTransport.cpp b/Cantera/src/transport/SimpleTransport.cpp index b028cdf7a..ae1efb3dc 100644 --- a/Cantera/src/transport/SimpleTransport.cpp +++ b/Cantera/src/transport/SimpleTransport.cpp @@ -818,8 +818,8 @@ namespace Cantera { rhoVc[n] += fluxes[n*ldf + k] / mw[k]; } } - for (n = 0; n < m_nDim; n++) { - for (k = 0; k < m_nsp; k++) { + for (size_t n = 0; n < m_nDim; n++) { + for (size_t k = 0; k < m_nsp; k++) { fluxes[n*ldf + k] -= m_molefracs[k] * rhoVc[n] * mw[k]; } fluxes[n*ldf + m_velocityBasis] = 0.0; diff --git a/Cantera/src/transport/TransportFactory.cpp b/Cantera/src/transport/TransportFactory.cpp index af14c419a..bc7cacdcf 100644 --- a/Cantera/src/transport/TransportFactory.cpp +++ b/Cantera/src/transport/TransportFactory.cpp @@ -1462,8 +1462,6 @@ namespace Cantera { doublereal eps, sigma; for (size_t k = 0; k < tr.nsp_; k++) { for (size_t j = k; j < tr.nsp_; j++) { - - ipoly = tr.poly[k][j]; for (size_t n = 0; n < np; n++) { t = tr.tmin + dt*n;