From 7c9bf2f41eeee51f94982cd408c0b5cc46173689 Mon Sep 17 00:00:00 2001 From: DetlevCM Date: Fri, 3 Jul 2015 18:55:53 +0200 Subject: [PATCH] updates to structure & output updates to structure & output --- headers/Get_Input.h | 3 +- headers/Global_Variables.hpp | 18 + headers/MyFunctions.h | 12 + headers/MyHeaders.h | 2 + headers/Namespaces.hpp | 25 + headers/Solver_Calculations.h | 6 +- headers/Structs.h | 28 +- headers/Write_Output.h | 2 +- source/Main.cpp | 917 +++--------------- source/Make_Irreversible_v2.cpp | 8 +- source/Run_Integrator.cpp | 704 ++++++++++++++ source/get_input/Get_Species.cpp | 28 +- source/get_input/Handle_Mechanism_Input.cpp | 70 +- source/get_input/Read_Input_Data_v2.cpp | 41 - source/get_input/Read_Reaction_Matrix.cpp | 4 +- .../Read_Thermodynamic_Data_New_Format.cpp | 66 +- .../Combine_Species-k-based.cpp | 18 +- .../Compare_Concentrations.cpp | 111 --- .../Get_Combine_Species_Mapping.cpp | 4 +- source/mechanism_reduction/Remove_Species.cpp | 4 +- .../Species_Class_Mapping.cpp | 11 +- .../Calculate_Rate_Constants.cpp | 7 +- .../Calculate_Thermodynamics.cpp | 19 +- source/write_output/WriteSpecies.cpp | 2 +- source/write_output/Write_Mechanism.cpp | 2 +- .../Write_Stoichiometry_Matrix_For_Opt.cpp | 2 +- 26 files changed, 1003 insertions(+), 1111 deletions(-) create mode 100644 headers/Global_Variables.hpp create mode 100644 headers/Namespaces.hpp create mode 100644 source/Run_Integrator.cpp diff --git a/headers/Get_Input.h b/headers/Get_Input.h index 64b3802..4b4da88 100644 --- a/headers/Get_Input.h +++ b/headers/Get_Input.h @@ -13,7 +13,8 @@ #define INPUT_FUNCTIONS_ -void Handle_Mechanism_Input( +//void Handle_Mechanism_Input( +bool Handle_Mechanism_Input( string , vector< string >& , vector< ThermodynamicData >& , diff --git a/headers/Global_Variables.hpp b/headers/Global_Variables.hpp new file mode 100644 index 0000000..a53f434 --- /dev/null +++ b/headers/Global_Variables.hpp @@ -0,0 +1,18 @@ +/* + * Global_Variables.hpp + * + * Created on: 02.07.2015 + * Author: DetlevCM + */ + +#ifndef SOURCE_GLOBAL_VARIABLES_HPP_ +#define SOURCE_GLOBAL_VARIABLES_HPP_ + + + +// Shoddy implementation of PetroOxy logic +extern PetroOxyCalculation PetroOxyData; +extern int OxyGasSpeciesID; +// end PetroOxy Logic + +#endif /* SOURCE_GLOBAL_VARIABLES_HPP_ */ diff --git a/headers/MyFunctions.h b/headers/MyFunctions.h index bc6e362..0cd4d78 100644 --- a/headers/MyFunctions.h +++ b/headers/MyFunctions.h @@ -103,4 +103,16 @@ vector< SingleReactionData > Make_Irreversible( ); +void Integrate( + string, + string, + string, + string, + vector< double >, + vector< string > , + vector< SingleReactionData > &, + InitParam, + vector< double >& + ); + #endif /* _MY_OTHER_FUNCTIONS_ */ diff --git a/headers/MyHeaders.h b/headers/MyHeaders.h index 24ad3bd..dead906 100644 --- a/headers/MyHeaders.h +++ b/headers/MyHeaders.h @@ -40,6 +40,8 @@ using std::istringstream; #include +#include +#include // Definitions of my functions #include diff --git a/headers/Namespaces.hpp b/headers/Namespaces.hpp new file mode 100644 index 0000000..55573d6 --- /dev/null +++ b/headers/Namespaces.hpp @@ -0,0 +1,25 @@ +/* + * Namespaces.hpp + * + * Created on: 02.07.2015 + * Author: DetlevCM + */ + +#ifndef SOURCE_NAMESPACES_HPP_ +#define SOURCE_NAMESPACES_HPP_ + + +namespace GlobalArrays{ +// Key parameters that define the whole reaction scheme - used globally via namespaces +// Not sure if this is a good place to put it... +extern vector< vector < str_RatesAnalysis > > RatesAnalysisData; +extern vector< TrackSpecies > ProductsForRatesAnalysis; + +} + + + + + + +#endif /* SOURCE_NAMESPACES_HPP_ */ diff --git a/headers/Solver_Calculations.h b/headers/Solver_Calculations.h index 50bde3e..e9d9fef 100644 --- a/headers/Solver_Calculations.h +++ b/headers/Solver_Calculations.h @@ -18,7 +18,8 @@ void Calculate_Rate_Constant( vector< double >& , const double , const vector< ReactionParameter >& , - const vector< vector< double > >& , + //const vector< vector< double > >& , + const vector< CalculatedThermodynamics >& , const vector< TrackSpecies >& , const vector< double >& ); @@ -37,7 +38,8 @@ void CalculateReactionRates( void Calculate_Thermodynamics( - vector< vector< double > >& , + //vector< vector< double > >& , + vector< CalculatedThermodynamics >& , const double& , const vector< ThermodynamicData >& ); diff --git a/headers/Structs.h b/headers/Structs.h index f81cf43..dc9c143 100644 --- a/headers/Structs.h +++ b/headers/Structs.h @@ -23,9 +23,9 @@ struct JacobianSpecies{ }; struct JacobianData { - int ColumnWiseArrayPosition; bool IsForward; bool IsProduction; + int ColumnWiseArrayPosition; int ReactionID; double coefficient; vector< JacobianSpecies > Species; @@ -34,22 +34,22 @@ struct JacobianData { struct SingleReactionData { - vector Reactants; - vector Products; + bool Reversible; + bool IsDuplicate; double paramA; double paramN; double paramEa; - bool Reversible; - bool IsDuplicate; + vector Reactants; + vector Products; }; //*/ //* struct ReactionParameter { + bool Reversible; double paramA; double paramN; double paramEa; - bool Reversible; }; //*/ @@ -71,6 +71,20 @@ struct InitSpecies { }; +struct ConstantInitRHSODE { + bool EnforceStability; + bool PetroOxy; + double temperature; + double PetroOxyTemperatureRise; +}; + +struct CalculatedThermodynamics { + double Hf; + double Cp; + double Cv; + double S; +}; + struct InitParam { bool irrev; // shoud be bool bool PrintReacRates; @@ -118,13 +132,13 @@ struct InitParam { struct PetroOxyCalculation { + bool HenryLawDiffusionLimitSet; double SampleSize ; // Old 0 double HeadSpaceGas ; // Old 4 double HeadSpaceGasMol; // Old 6 double HenryConstantk; // Old 10 double HeadSpaceGasPressure; // Old 7 double HeadSpaceSolventComponentPressure; // Old 8 - bool HenryLawDiffusionLimitSet; double HenryLawDiffusionLimit; }; diff --git a/headers/Write_Output.h b/headers/Write_Output.h index 7bb201d..e057abe 100644 --- a/headers/Write_Output.h +++ b/headers/Write_Output.h @@ -73,7 +73,7 @@ void WriteReactions( void Write_Stoichiometric_Matrix_For_Opt( string , - const vector< string >& , +// const vector< string >& , const vector< SingleReactionData >& ); diff --git a/source/Main.cpp b/source/Main.cpp index 8645682..eeefa67 100644 --- a/source/Main.cpp +++ b/source/Main.cpp @@ -1,118 +1,54 @@ #include -//struct SpeciesLoss; - -// The actual ODEs, only required locally, therefore defined in here only -extern void ODE_RHS(int*, double*, double*, double*); -extern void Jacobian_Matrix(int*, double*, double*, double*); - -/* It seems a limited number of global variables cannot be avoided. - * The reaction scheme as well as species concentrations will need - * to be provided as a global variable as the ODEs need access to - * them without passing any variables to the relevant void */ - - namespace GlobalArrays{ // Key parameters that define the whole reaction scheme - used globally via namespaces -InitParam InitialParameters; // Not sure if this is a good place to put it... vector< vector < str_RatesAnalysis > > RatesAnalysisData; vector< TrackSpecies > ProductsForRatesAnalysis; - -//vector< double > JacobeanColumnWise; -vector< JacobianData > JacobianMatrix; } -// Needed for right hand side equation for the solve to function - -namespace ODESolverConstant -{ -// Below Used in ODE RHS -int Number_Species; -int Number_Reactions; -vector< double > Delta_N; -vector< ThermodynamicData > Thermodynamics; -vector< ReactionParameter > ReactionParameters; // tidier than reactions vector -vector< TrackSpecies > ReactantsForReactions; -vector< TrackSpecies > ProductsForReactions; -vector< TrackSpecies > SpeciesLossAll; // vector for recording species loss -} - -namespace ODESolverVariable -{ -vector< vector< double > > CalculatedThermo; -vector< double > Kf; -vector< double > Kr; -vector< double > Rates; -vector< double > Concentration; -vector< double > SpeciesConcentrationChange; - - -//for limited Oxy -double time_previous = 0; -} - -// End solver function - - - -// By making the following global, I avoid vector re-declarations -// This has a positive impact on performance -// by sticking them in a namespace I make them more manageable -// It also allows me to size the array beforehand - // Shoddy implementation of PetroOxy logic PetroOxyCalculation PetroOxyData; int OxyGasSpeciesID; // end PetroOxy Logic - - -void Integrate( - string, - string, - string, - string, - vector< double >, - vector< string > , - vector< SingleReactionData >& -); - - -// Global for now, clean up !!!!!!!!!!! -vector< double > KeyRates; - - - int main(int argc, char* argv[]) { - // http://stackoverflow.com/questions/10150468/how-to-redirect-cin-and-cout-to-files - std::ofstream out("log.txt"); - //std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf - std::cout.rdbuf(out.rdbuf()); - - vector< string > Species; - vector< SingleReactionData > Reactions; + // Add a switch to have regular ouput to a log file or debug to command line + // http://stackoverflow.com/questions/12774207/fastest-way-to-check-if-a-file-exist-using-standard-c-c11-c + ifstream file_exists("debug"); + // File does not exist, so redirect to log file + // http://stackoverflow.com/questions/10150468/how-to-redirect-cin-and-cout-to-files + std::ofstream out("log.txt"); + if (!file_exists.good()){ + //std::streambuf *coutbuf = std::cout.rdbuf(); //save old buf + std::cout.rdbuf(out.rdbuf()); + } + file_exists.close(); - using namespace GlobalArrays; - using namespace ODESolverConstant; // what is needed? - PetroOxyCalculation PetroOxyDataInitial; - // Globally useful counter - int i; + // Definition of the main variables to store basic information to run the sover + vector< string > Species; // stores the species list + vector< ThermodynamicData > Thermodynamics; // stores the entire Thermodynamic data + vector< SingleReactionData > Reactions; // stores the entire mechanism + InitParam InitialParameters; // Initial Conditions/Parameters vector< double > InitialSpeciesConcentration; + PetroOxyCalculation PetroOxyDataInitial; // Petro Oxy Specific + + + // Handle All the Data Input - The Arrays Contain the required information - Handle_Mechanism_Input( + bool Mechanism_Read_In = Handle_Mechanism_Input( "chem.inp", Species, Thermodynamics, @@ -122,767 +58,198 @@ int main(int argc, char* argv[]) PetroOxyDataInitial ); - if(InitialParameters.StoichimoetryMatrixForOpt) - { - Write_Stoichiometric_Matrix_For_Opt - ( - "stoichiometry_matrix.txt" , - Species, - Reactions - ); - } - - Number_Species = (int) Species.size(); - Number_Reactions = (int) Reactions.size(); - - // We have now pre-processed all information, time to set up the ODEs and the solver - // Let us set up the reactions first for the ODE solver - - - printf("On to the solver. \n"); - - - if(InitialParameters.RatesMaxAnalysis) - {// Initialise array - vector < str_RatesAnalysis > temp((int) Reactions.size()); - for(i=0;i ReducedReactions; - ReducedReactions = ReduceReactionsNew(Species, Reactions, KeyRates); - - // start a second run only if reduced scheme is not empty and has size different - // to original scheme - if(Number_Reactions > 0 && Number_Reactions != (int) ReducedReactions.size()){ - Number_Reactions = (int) ReducedReactions.size(); - - WriteReactions("reduced_scheme.txt", Species, ReducedReactions); - - InitialParameters.ReduceReactions = 0; // switch off reduction... - - if(InitialParameters.StoichimoetryMatrixForOpt) - { - Write_Stoichiometric_Matrix_For_Opt - ( - "reduced_stoichiometry_matrix.txt" , - Species, - ReducedReactions - ); - } - - - // Option considered experimental, cannot see why it won't work... - if(InitialParameters.RatesMaxAnalysis) - {// Initialise array - RatesAnalysisData.clear(); // empty for new run - vector < str_RatesAnalysis > temp(Number_Reactions); - for(i=0;i SpeciesConcentration, - vector< string > Species, - vector< SingleReactionData >& Reactions -) -{ - using namespace GlobalArrays; - - using namespace ODESolverConstant; - using namespace ODESolverVariable; - - - ofstream ReactionRatesOutput; - ofstream ConcentrationOutput (species_filename.c_str(),ios::app); - - if(InitialParameters.PrintReacRates) - { - ReactionRatesOutput.open(rates_filename.c_str(),ios::app); - } - - - - int n, ierr, ipar[128], dpar_size; - double h, hm, ep, tr; - int i; - - /* Because we do not know the size of kd and dpar in advance we need to be clever - or at least act as we are. - * By defining kd and dpar as vectors we can assign a size "on the fly" as required. We can further pass a - * pointer to the vector to the solver, thus this is what we do. Rather than arrays, kd and dpar are pointers - * to vectors of doubles in our case. - */ - - - vector kdstep(Number_Species + 1); - - - int* kd = &kdstep[0]; - - n = Number_Species + 1; - if (13 * n > (7 + 2 * n) * n) { - dpar_size = 13 * n; - } else { - dpar_size = (7 + 2 * n) * n; - } - vector dparstep(dpar_size); - double* dpar = &dparstep[0]; - - // For performance assessment, use a clock: - clock_t cpu_time_begin, cpu_time_end, cpu_time_current; - - // Some tolerances for the solver: - hm = InitialParameters.hm; // minimal step size for the methods, 1.0e-12 recommended for normalised problems - ep = InitialParameters.rtol ; // relative tolerance. The code cannot guarantee the requested accuracy for ep<1.d-9 - tr = InitialParameters.threshold; // Threshold, absolute tolerance is ep*tr - h = InitialParameters.h; - - - Delta_N = Get_Delta_N(Reactions); // just make sure the Delta_N is current - // Reduce the matrix from a sparse matrix to something more manageable and quicker to use - - - ReactionParameters = Process_Reaction_Parameters(Reactions); - - ReactantsForReactions = Reactants_ForReactionRate(Reactions); - - ProductsForReactions = Products_ForReactionRate(Reactions,false); - - if(InitialParameters.RatesMaxAnalysis || InitialParameters.StreamRatesAnalysis || InitialParameters.RatesAnalysisAtTime) - { - ProductsForRatesAnalysis = Products_ForReactionRate(Reactions,true); - } - - SpeciesLossAll = PrepareSpecies_ForSpeciesLoss(Reactions); // New method of listing species - - // And now it is time to call the solver again... with the right information... - // According to Intel initialise ipar array with zeros before the first call to dodesol - for (i = 0; i < 128; i++) { - ipar[i] = 0; - } - - double* y = &SpeciesConcentration[0]; - double time_current, time_step, time_step1, time_end; - - time_current = 0;// -> Solver is designed for t_0 = 0 - time_step1 = InitialParameters.TimeStep[0]; - time_end = InitialParameters.TimeEnd[0]; - int TimeChanges = (int) InitialParameters.TimeStep.size(); - int tracker = 0; - - cout << "\nEnd Time: " << time_end << " Time Step: " << time_step1 << "\n"; - /* -- Initial values at t = 0 -- */ - - Number_Reactions = (int) ReactionParameters.size(); - - CalculatedThermo.resize(Number_Species); - for(i=0;i SpeciesConcentrationChange = SpeciesLossRate(SpeciesLossAll,Number_Species, Rates); - - - /* -- Got values at t = 0 -- */ - - // Resize the vector for the concentrations in the ODE void - Concentration.clear(); // don't think it makes a difference if I clear or not - Concentration.resize(Number_Species + 1); - - // enables reset of Rates Analysis - int RatesAnalysisTimepoint = 0; - - - // start the clock: - cpu_time_begin = cpu_time_current = clock(); - - - do - { - time_step = time_current + time_step1; - - - if(!InitialParameters.UseStiffSolver && InitialParameters.Jacobian) - { - dodesol_rkm9mka(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, (void*) Jacobian_Matrix, &h, &hm, &ep, &tr, dpar, kd, &ierr); + Write_Stoichiometric_Matrix_For_Opt + ( + "stoichiometry_matrix.txt" , + Reactions + ); } - if(InitialParameters.UseStiffSolver && InitialParameters.Jacobian) - { - // stiff solver with automatics numerical Jacobi matric computations - dodesol_mk52lfa(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, (void*) Jacobian_Matrix, &h, &hm, &ep, &tr, dpar, kd, &ierr); - } + using namespace GlobalArrays; - if(InitialParameters.UseStiffSolver && !InitialParameters.Jacobian) - { - // stiff solver with automatics numerical Jacobi matric computations - dodesol_mk52lfn(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, &h, &hm, &ep, &tr, dpar, kd, &ierr); - } + int i; // useful counter + int Number_Species = (int) Species.size(); + int Number_Reactions = (int) Reactions.size(); + vector< double > KeyRates; // for mechanism reduction - if(!InitialParameters.UseStiffSolver && !InitialParameters.Jacobian) - { - // hybrid solver with automatic numerical Jacobi matrix computations - dodesol_rkm9mkn(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, &h, &hm, &ep, &tr, dpar, kd, &ierr); - } - if (ierr != 0) - { - printf("\n dodesol_rkm9mkn routine exited with error code %4d\n",ierr); - //return -1; - // Break means it should leave the do loop which would be fine for an error response as it stops the solver - break ; - } + // We have now pre-processed all information, time to set up the ODEs and the solver + // Let us set up the reactions first for the ODE solver + //cout << "On to the solver. \n"; if(InitialParameters.RatesMaxAnalysis) { - MaxRatesAnalysis(RatesAnalysisData,ProductsForRatesAnalysis,ReactantsForReactions,Rates,time_current); - } - - - if(InitialParameters.RatesAnalysisAtTime && - InitialParameters.RatesAnalysisAtTimeData[RatesAnalysisTimepoint] == time_current) - { - using namespace GlobalArrays; - RatesAnalysisAtTimes( - ProductsForRatesAnalysis, - ReactantsForReactions, - Rates, - time_current, - Species, - Reactions - ); - - RatesAnalysisTimepoint = RatesAnalysisTimepoint + 1; + // Initialise array + vector < str_RatesAnalysis > temp((int) Reactions.size()); + for(i=0;i= InitialParameters.TimeEnd[tracker]) - { - cout << "CPU Time: " << ((double) (clock() - cpu_time_current)) / CLOCKS_PER_SEC << " seconds\n"; - cpu_time_current = clock(); - - tracker = tracker + 1; - time_step1 = InitialParameters.TimeStep[tracker]; - time_end = InitialParameters.TimeEnd[tracker]; - cout << "End Time: " << time_end << " Time Step: " << time_step1 << "\n"; + KeyRates.resize(Number_Reactions); } - } while (time_step < time_end); - - - cout << "CPU Time: " << ((double) (clock() - cpu_time_current)) / CLOCKS_PER_SEC << " seconds\n"; - - // close output files - ConcentrationOutput.close(); - ReactionRatesOutput.close(); - - if(InitialParameters.RatesMaxAnalysis) - { - WriteMaxRatesAnalysis(RatesAnalysisData, Species, Reactions,rates_analysis_stream_filename); - } - - - // stop the clock - cpu_time_end = clock(); - cout << "\nTotal CPU time: " << ((double) (cpu_time_end - cpu_time_begin)) / CLOCKS_PER_SEC << " seconds\n"; -} - - - - + cout << "\nHanding Mechanism to Integrator\n" << std::flush; -void ODE_RHS(int*n, double*time_current, double*y, double*f) -{ - // A namespace allows global variables without causing a mess, should be quicker than redefining too - using namespace ODESolverConstant; - using namespace ODESolverVariable; - - int i; + Integrate( + "concentrations.txt", + "reaction_rates.txt", + "PetroOxy-Log.txt", + "", + InitialSpeciesConcentration, + Species, + Reactions, + InitialParameters, + KeyRates + ); - /* 2002 CODATA values */ - //double R = 8.314472e0; - //double Na = 6.0221415e23; - // stability hack - if(GlobalArrays::InitialParameters.EnforceStability) - { - for (i = 0; i <= Number_Species; i++) - { - if(y[i]<0){ - //if(y[i]<1.e-24){ - Concentration[i] = 0; - } - else - { - Concentration[i] = y[i]; - } - } - } - else{ - for (i = 0; i <= Number_Species; i++) + if(InitialParameters.ReduceReactions != 0) { - Concentration[i] = y[i]; - } - } + vector< SingleReactionData > ReducedReactions; + ReducedReactions = ReduceReactionsNew(Species, Reactions, KeyRates); - double time_difference; - - - if(GlobalArrays::InitialParameters.PetroOxy) - { - time_difference = fabs(*time_current - time_previous); - AdjustGasConcentration( - y[OxyGasSpeciesID], - Concentration[Number_Species], - time_difference, - PetroOxyData); - } - - - // Thermodynamic data, Rate Constant, Rates, new Concentrations - Calculate_Thermodynamics(CalculatedThermo, Concentration[Number_Species], Thermodynamics); - Calculate_Rate_Constant(Kf, Kr, Concentration[Number_Species],ReactionParameters, CalculatedThermo, SpeciesLossAll, Delta_N); - CalculateReactionRates(Rates, Concentration, Kf, Kr, ReactantsForReactions, ProductsForReactions); - SpeciesConcentrationChange = SpeciesLossRate(SpeciesLossAll,Number_Species, Rates); + // start a second run only if reduced scheme is not empty and has size different + // to original scheme + if(Number_Reactions > 0 && Number_Reactions != (int) ReducedReactions.size()){ + Number_Reactions = (int) ReducedReactions.size(); + WriteReactions("reduced_scheme.txt", Species, ReducedReactions); - double ctot=0; - double qint=0; - double qtot=0; - - - for (i = 0; i < Number_Species; i++) - { - ctot = ctot + CalculatedThermo[i][2] * Concentration[i]; - } - // ctot = ctot / 1000; // working in moles/l so no Na; - - for (i = 0; i < Number_Reactions; i++) - { - qint = qint + Delta_N[i] * Rates[i]; - } - qtot = -qint / (ctot);//*1000); // scale l to ml and Na not needed for moles/l * Na); //*/ - - // Checked f[i] -> no unexpected rates of change for "inert gases", all 0. - for (i = 0; i < Number_Species; i++) - { - f[i] = SpeciesConcentrationChange[i]; // Species equations // - } - f[Number_Species] = qtot; // Temperature equation // - - //* - if( - Concentration[Number_Species] >= GlobalArrays::InitialParameters.temperature - && - GlobalArrays::InitialParameters.PetroOxy) // fix temperature for Oxy - { - y[Number_Species] = GlobalArrays::InitialParameters.temperature; // ensure temperature is not exceeded - f[Number_Species] = 0; - GlobalArrays::InitialParameters.PetroOxyTemperatureRise = 0; - }//*/ - - if( - GlobalArrays::InitialParameters.PetroOxy - && - GlobalArrays::InitialParameters.PetroOxyTemperatureRise != 0) // fix temperature for Oxy, rise desired - { - // 298K starting temp, GlobalArrays::InitialParameters.temperature is target - // rise time given in minutes - f[Number_Species] = - (GlobalArrays::InitialParameters.temperature - 298) - / - (GlobalArrays::InitialParameters.PetroOxyTemperatureRise); - //std::cout << f[Number_Species] << "\n"; - }//*/ - - // IEEE standard hack to check for NaN - // if temperature blows up, freeze it - if(qtot != qtot) - { - f[Number_Species] = 0; - } + InitialParameters.ReduceReactions = 0; // switch off reduction... + if(InitialParameters.StoichimoetryMatrixForOpt) + { + Write_Stoichiometric_Matrix_For_Opt + ( + "reduced_stoichiometry_matrix.txt" , + // Species, + ReducedReactions + ); + } - // For Oxy limiting - time_previous = *time_current; -} - - - - -void Jacobian_Matrix(int*n,double*t,double*y,double*a) { - - // n -> number of species - // t time - // y concentration - // a Jacobian in column wise order - using namespace ODESolverConstant; - using namespace ODESolverVariable; - using namespace GlobalArrays; + // Option considered experimental, cannot see why it won't work... + if(InitialParameters.RatesMaxAnalysis) + { + // Initialise array + RatesAnalysisData.clear(); // empty for new run + vector < str_RatesAnalysis > temp(Number_Reactions); + for(i=0;i JacobeanColumnWise((Number_Species+1)*(Number_Species+1)); + if(InitialParameters.PetroOxy) + { + PetroOxyData = PetroOxyDataInitial; + OxyGasSpeciesID = InitialParameters.PetroOxyGasSpecies; + PetroOxyOutputHeader("reduced_PetroOxy-Log.txt"); + } - Calculate_Thermodynamics(CalculatedThermo, Concentration[Number_Species], Thermodynamics); - Calculate_Rate_Constant(Kf, Kr, Concentration[Number_Species],ReactionParameters, CalculatedThermo, SpeciesLossAll, Delta_N); + cout << "\nHanding Reduced Mechanism to Integrator\n" << std::flush; - for(i=0;i<(int) JacobianMatrix.size();i++) - { - double temp; + Integrate( + "reduced_concentrations.txt", + "reduced_reaction_rates.txt", + "reduced_PetroOxy-Log.txt", + "reduced_", + InitialSpeciesConcentration, + Species, + ReducedReactions, + InitialParameters, + KeyRates + ); - if(JacobianMatrix[i].IsForward) // Forward - { - if(JacobianMatrix[i].IsProduction) // species are gained - { - temp = Kf[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; - } - if(!JacobianMatrix[i].IsProduction) // species are lost - { - temp = -Kf[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; - } - } + ReportAccuracy( + Number_Species, + Species, + "reduction_accuracy_report.txt", + "concentrations.txt", + "reduced_concentrations.txt" + ); - if(!JacobianMatrix[i].IsForward) // Reverse - { - if(JacobianMatrix[i].IsProduction) // species are gained - { - temp = Kr[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; } - if(!JacobianMatrix[i].IsProduction) // species are lost + else { - temp = -Kr[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; - } - } - + cout << + "Reduced Scheme does not contain any reactions or is identical in size to the original scheme." + " \n Aborting...!!!"; + }; - for(j=0;j<(int) JacobianMatrix[i].Species.size();j++) - { - if(JacobianMatrix[i].Species[j].power != 0) // power 0 = *1 - { - if(JacobianMatrix[i].Species[j].power == 1) // power 1 is simple multiplication - { - temp = temp * Concentration[JacobianMatrix[i].Species[j].SpeciesID]; - } - else - { - temp = temp * - pow(Concentration[JacobianMatrix[i].Species[j].SpeciesID], - JacobianMatrix[i].Species[j].power); - } - } } - - JacobeanColumnWise[JacobianMatrix[i].ColumnWiseArrayPosition] = - JacobeanColumnWise[JacobianMatrix[i].ColumnWiseArrayPosition] + temp; } - for (i = 0; i <= (int) JacobeanColumnWise.size() ; i++) - { - a[i] = JacobeanColumnWise[i]; - } + // And it is done. Return 0 as the code has finished. + return 0; +} - /* - * Debug Output - */ - /* - cout << "\n" << "\n"; - int temp_output = (int) sqrt(JacobeanColumnWise.size()); - for(i=0;i< temp_output;i++) - { - for(j=0;j<(int) JacobeanColumnWise.size();j++) - { - if(j % temp_output == i) - { - cout << a[j] << " "; - //cout << JacobeanColumnWise[j] << " "; - //cout << j % temp_output << " "; - } - } - cout << "\n"; - } - cout << "\n" << "\n"; - //*/ - //cout << "check\n"; -} diff --git a/source/Make_Irreversible_v2.cpp b/source/Make_Irreversible_v2.cpp index cb3e020..56207e3 100644 --- a/source/Make_Irreversible_v2.cpp +++ b/source/Make_Irreversible_v2.cpp @@ -55,11 +55,13 @@ vector< SingleReactionData > Make_Irreversible( // store the k for every reaction and every temperature //cout << "checkpoint 1 \n"; - vector< vector< double > > CalculatedThermo(Number_Species); - for(i=0;i > CalculatedThermo(Number_Species); + vector< CalculatedThermodynamics > CalculatedThermo(Number_Species); + + /*for(i=0;i Kf(Number_Reactions) ; vector< double > Kr(Number_Reactions) ; diff --git a/source/Run_Integrator.cpp b/source/Run_Integrator.cpp new file mode 100644 index 0000000..79a31b4 --- /dev/null +++ b/source/Run_Integrator.cpp @@ -0,0 +1,704 @@ +/* + * Run_Integrator.cpp + * + * Created on: 02.07.2015 + * Author: DetlevCM + */ +#include + + + + +// The actual ODEs, only required locally, therefore defined in here only +extern void ODE_RHS(int*, double*, double*, double*); +extern void Jacobian_Matrix(int*, double*, double*, double*); + +/* It seems a limited number of global variables cannot be avoided. + * The reaction scheme as well as species concentrations will need + * to be provided as a global variable as the ODEs need access to + * them without passing any variables to the relevant void */ + + + + +// Needed for right hand side equation for the solve to function + +namespace ODESolverConstant +{ +// Below Used in ODE RHS +int Number_Species; +int Number_Reactions; +vector< double > Delta_N; +vector< ThermodynamicData > Thermodynamics; +vector< ReactionParameter > ReactionParameters; // tidier than reactions vector +vector< TrackSpecies > ReactantsForReactions; +vector< TrackSpecies > ProductsForReactions; +vector< TrackSpecies > SpeciesLossAll; // vector for recording species loss + +// No more initial parameters globally +ConstantInitRHSODE InitialDataConstants; + +// Calculated Once Only +vector< JacobianData > JacobianMatrix; +} + +namespace ODESolverVariable +{ +vector< CalculatedThermodynamics > CalculatedThermo; +vector< double > Kf; +vector< double > Kr; +vector< double > Rates; +vector< double > Concentration; +vector< double > SpeciesConcentrationChange; + + +//for limited Oxy +double time_previous = 0; +} + +// End solver function + + + +// By making the following global, I avoid vector re-declarations +// This has a positive impact on performance +// by sticking them in a namespace I make them more manageable +// It also allows me to size the array beforehand + + + + + + + + + + + + +// Not a perfect solution, but stick integrator into its own void with global variables via a namespace +void Integrate( + string species_filename, + string rates_filename, + string petrooxy_filename, + string rates_analysis_stream_filename, + vector< double > SpeciesConcentration, + vector< string > Species, + vector< SingleReactionData >& Reactions, + InitParam InitialParameters, + vector< double >& KeyRates +) +{ + + + + using namespace GlobalArrays; + + using namespace ODESolverConstant; + using namespace ODESolverVariable; + + Number_Species = Species.size(); + Number_Reactions = Reactions.size(); + + + ofstream ReactionRatesOutput; + ofstream ConcentrationOutput (species_filename.c_str(),ios::app); + + if(InitialParameters.PrintReacRates) + { + ReactionRatesOutput.open(rates_filename.c_str(),ios::app); + } + + + + int n, ierr, ipar[128], dpar_size; + double h, hm, ep, tr; + int i; + + /* Because we do not know the size of kd and dpar in advance we need to be clever - or at least act as we are. + * By defining kd and dpar as vectors we can assign a size "on the fly" as required. We can further pass a + * pointer to the vector to the solver, thus this is what we do. Rather than arrays, kd and dpar are pointers + * to vectors of doubles in our case. + */ + + + vector kdstep(Number_Species + 1); + + + int* kd = &kdstep[0]; + + n = Number_Species + 1; + if (13 * n > (7 + 2 * n) * n) { + dpar_size = 13 * n; + } else { + dpar_size = (7 + 2 * n) * n; + } + vector dparstep(dpar_size); + double* dpar = &dparstep[0]; + + // For performance assessment, use a clock: + clock_t cpu_time_begin, cpu_time_end, cpu_time_current; + + + + // Some tolerances for the solver: + hm = InitialParameters.hm; // minimal step size for the methods, 1.0e-12 recommended for normalised problems + ep = InitialParameters.rtol ; // relative tolerance. The code cannot guarantee the requested accuracy for ep<1.d-9 + tr = InitialParameters.threshold; // Threshold, absolute tolerance is ep*tr + h = InitialParameters.h; + + + Delta_N = Get_Delta_N(Reactions); // just make sure the Delta_N is current + // Reduce the matrix from a sparse matrix to something more manageable and quicker to use + + + ReactionParameters = Process_Reaction_Parameters(Reactions); + + ReactantsForReactions = Reactants_ForReactionRate(Reactions); + + ProductsForReactions = Products_ForReactionRate(Reactions,false); + + if(InitialParameters.RatesMaxAnalysis || InitialParameters.StreamRatesAnalysis || InitialParameters.RatesAnalysisAtTime) + { + ProductsForRatesAnalysis = Products_ForReactionRate(Reactions,true); + } + + SpeciesLossAll = PrepareSpecies_ForSpeciesLoss(Reactions); // New method of listing species + + // And now it is time to call the solver again... with the right information... + // According to Intel initialise ipar array with zeros before the first call to dodesol + for (i = 0; i < 128; i++) { + ipar[i] = 0; + } + + + + double* y = &SpeciesConcentration[0]; + double time_current, time_step, time_step1, time_end; + + time_current = 0;// -> Solver is designed for t_0 = 0 + time_step1 = InitialParameters.TimeStep[0]; + time_end = InitialParameters.TimeEnd[0]; + int TimeChanges = (int) InitialParameters.TimeStep.size(); + int tracker = 0; + + cout << "\nEnd Time: " << time_end << " Time Step: " << time_step1 << "\n"; + /* -- Initial values at t = 0 -- */ + + Number_Reactions = (int) ReactionParameters.size(); + + CalculatedThermo.resize(Number_Species); + /* + for(i=0;i SpeciesConcentrationChange = SpeciesLossRate(SpeciesLossAll,Number_Species, Rates); + + + /* -- Got values at t = 0 -- */ + + // Resize the vector for the concentrations in the ODE void + Concentration.clear(); // don't think it makes a difference if I clear or not + Concentration.resize(Number_Species + 1); + + // enables reset of Rates Analysis + int RatesAnalysisTimepoint = 0; + + + // start the clock: + cpu_time_begin = cpu_time_current = clock(); + + + do + { + time_step = time_current + time_step1; + + + if(!InitialParameters.UseStiffSolver && InitialParameters.Jacobian) + { + dodesol_rkm9mka(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, (void*) Jacobian_Matrix, &h, &hm, &ep, &tr, dpar, kd, &ierr); + } + + if(InitialParameters.UseStiffSolver && InitialParameters.Jacobian) + { + // stiff solver with automatics numerical Jacobi matric computations + dodesol_mk52lfa(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, (void*) Jacobian_Matrix, &h, &hm, &ep, &tr, dpar, kd, &ierr); + } + + if(InitialParameters.UseStiffSolver && !InitialParameters.Jacobian) + { + // stiff solver with automatics numerical Jacobi matric computations + dodesol_mk52lfn(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, &h, &hm, &ep, &tr, dpar, kd, &ierr); + } + + if(!InitialParameters.UseStiffSolver && !InitialParameters.Jacobian) + { + // hybrid solver with automatic numerical Jacobi matrix computations + dodesol_rkm9mkn(ipar, &n, &time_current, &time_step, y, (void*) ODE_RHS, &h, &hm, &ep, &tr, dpar, kd, &ierr); + } + if (ierr != 0) + { + printf("\n dodesol_rkm9mkn routine exited with error code %4d\n",ierr); + //return -1; + // Break means it should leave the do loop which would be fine for an error response as it stops the solver + break ; + } + + + if(InitialParameters.RatesMaxAnalysis) + { + MaxRatesAnalysis(RatesAnalysisData,ProductsForRatesAnalysis,ReactantsForReactions,Rates,time_current); + } + + + if(InitialParameters.RatesAnalysisAtTime && + InitialParameters.RatesAnalysisAtTimeData[RatesAnalysisTimepoint] == time_current) + { + using namespace GlobalArrays; + RatesAnalysisAtTimes( + ProductsForRatesAnalysis, + ReactantsForReactions, + Rates, + time_current, + Species, + Reactions + ); + + RatesAnalysisTimepoint = RatesAnalysisTimepoint + 1; + } + + + if(InitialParameters.StreamRatesAnalysis) + { + StreamRatesAnalysis( + rates_analysis_stream_filename, + ProductsForRatesAnalysis, + ReactantsForReactions, + Rates, + time_current, + Number_Species + ); + } + + double pressure = 0; + if(InitialParameters.GasPhase) + { + double R = 8.31451; // Gas Constant + /* Pressure Tracking for Gas Phase Kinetics */ + double total_mol = 0; + for(i=0;i= InitialParameters.TimeEnd[tracker]) + { + cout << "CPU Time: " << ((double) (clock() - cpu_time_current)) / CLOCKS_PER_SEC << " seconds\n"; + cpu_time_current = clock(); + + tracker = tracker + 1; + time_step1 = InitialParameters.TimeStep[tracker]; + time_end = InitialParameters.TimeEnd[tracker]; + cout << "End Time: " << time_end << " Time Step: " << time_step1 << "\n"; + } + + + } while (time_step < time_end); + + + cout << "CPU Time: " << ((double) (clock() - cpu_time_current)) / CLOCKS_PER_SEC << " seconds\n"; + + // close output files + ConcentrationOutput.close(); + ReactionRatesOutput.close(); + + if(InitialParameters.RatesMaxAnalysis) + { + WriteMaxRatesAnalysis(RatesAnalysisData, Species, Reactions,rates_analysis_stream_filename); + } + + + // stop the clock + cpu_time_end = clock(); + cout << "\nTotal CPU time: " << ((double) (cpu_time_end - cpu_time_begin)) / CLOCKS_PER_SEC << " seconds\n\n"; +} + + + + + +void ODE_RHS(int*n, double*time_current, double*y, double*f) +{ + // A namespace allows global variables without causing a mess, should be quicker than redefining too + using namespace ODESolverConstant; + using namespace ODESolverVariable; + + int i; + + + /* 2002 CODATA values */ + //double R = 8.314472e0; + //double Na = 6.0221415e23; + + // stability hack + //if(GlobalArrays::InitialParameters.EnforceStability) + if(InitialDataConstants.EnforceStability) + { + for (i = 0; i <= Number_Species; i++) + { + if(y[i]<0){ + //if(y[i]<1.e-24){ + Concentration[i] = 0; + } + else + { + Concentration[i] = y[i]; + } + } + } + else{ + for (i = 0; i <= Number_Species; i++) + { + Concentration[i] = y[i]; + } + } + + double time_difference; + + + //if(GlobalArrays::InitialParameters.PetroOxy) + if(InitialDataConstants.PetroOxy) + { + time_difference = fabs(*time_current - time_previous); + AdjustGasConcentration( + y[OxyGasSpeciesID], + Concentration[Number_Species], + time_difference, + PetroOxyData); + } + + // Thermodynamic data, Rate Constant, Rates, new Concentrations + Calculate_Thermodynamics(CalculatedThermo, Concentration[Number_Species], Thermodynamics); + Calculate_Rate_Constant(Kf, Kr, Concentration[Number_Species],ReactionParameters, CalculatedThermo, SpeciesLossAll, Delta_N); + CalculateReactionRates(Rates, Concentration, Kf, Kr, ReactantsForReactions, ProductsForReactions); + SpeciesConcentrationChange = SpeciesLossRate(SpeciesLossAll,Number_Species, Rates); + + + double ctot=0; + double qint=0; + double qtot=0; + + + for (i = 0; i < Number_Species; i++) + { + ctot = ctot + CalculatedThermo[i].Cv * Concentration[i]; + } + // ctot = ctot / 1000; // working in moles/l so no Na; + + for (i = 0; i < Number_Reactions; i++) + { + qint = qint + Delta_N[i] * Rates[i]; + } + qtot = -qint / (ctot);//*1000); // scale l to ml and Na not needed for moles/l * Na); //*/ + + // Checked f[i] -> no unexpected rates of change for "inert gases", all 0. + for (i = 0; i < Number_Species; i++) + { + f[i] = SpeciesConcentrationChange[i]; // Species equations // + } + f[Number_Species] = qtot; // Temperature equation // + + //* + if( + Concentration[Number_Species] >= InitialDataConstants.temperature + && + InitialDataConstants.PetroOxy) // fix temperature for Oxy + { + y[Number_Species] = InitialDataConstants.temperature; // ensure temperature is not exceeded + InitialDataConstants.PetroOxyTemperatureRise = 0; + + f[Number_Species] = 0; + }//*/ + + if( + InitialDataConstants.PetroOxy + && + InitialDataConstants.PetroOxyTemperatureRise != 0) // fix temperature for Oxy, rise desired + { + // 298K starting temp, GlobalArrays::InitialParameters.temperature is target + // rise time given in minutes + f[Number_Species] = + (InitialDataConstants.temperature - 298) + / + (InitialDataConstants.PetroOxyTemperatureRise); + + //std::cout << f[Number_Species] << "\n"; + }//*/ + + // IEEE standard hack to check for NaN + // if temperature blows up, freeze it + if(qtot != qtot) + { + f[Number_Species] = 0; + } + + + // For Oxy limiting + time_previous = *time_current; +} + + + + +void Jacobian_Matrix(int*n,double*t,double*y,double*a) { + + // n -> number of species + // t time + // y concentration + // a Jacobian in column wise order + + using namespace ODESolverConstant; + using namespace ODESolverVariable; + using namespace GlobalArrays; + + + int i,j; + + // enable force stability? + /* + for (i = 0; i <= Number_Species; i++) + { + if(y[i]<0){ + //if(y[i]<1.e-24){ + Concentration[i] = 0; + } + else + { + Concentration[i] = y[i]; + } + }//*/ + + for (i = 0; i <= Number_Species; i++) + { + Concentration[i] = y[i]; + } + + + + // provides me a fresh array every time :) - ideal + vector< double > JacobeanColumnWise((Number_Species+1)*(Number_Species+1)); + + + + Calculate_Thermodynamics(CalculatedThermo, Concentration[Number_Species], Thermodynamics); + Calculate_Rate_Constant(Kf, Kr, Concentration[Number_Species],ReactionParameters, CalculatedThermo, SpeciesLossAll, Delta_N); + + for(i=0;i<(int) JacobianMatrix.size();i++) + { + double temp; + + + if(JacobianMatrix[i].IsForward) // Forward + { + if(JacobianMatrix[i].IsProduction) // species are gained + { + temp = Kf[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; + } + if(!JacobianMatrix[i].IsProduction) // species are lost + { + temp = -Kf[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; + } + } + + if(!JacobianMatrix[i].IsForward) // Reverse + { + if(JacobianMatrix[i].IsProduction) // species are gained + { + temp = Kr[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; + } + if(!JacobianMatrix[i].IsProduction) // species are lost + { + temp = -Kr[JacobianMatrix[i].ReactionID]*JacobianMatrix[i].coefficient; + } + } + + + for(j=0;j<(int) JacobianMatrix[i].Species.size();j++) + { + if(JacobianMatrix[i].Species[j].power != 0) // power 0 = *1 + { + if(JacobianMatrix[i].Species[j].power == 1) // power 1 is simple multiplication + { + temp = temp * Concentration[JacobianMatrix[i].Species[j].SpeciesID]; + } + else + { + temp = temp * + pow(Concentration[JacobianMatrix[i].Species[j].SpeciesID], + JacobianMatrix[i].Species[j].power); + } + } + } + + + JacobeanColumnWise[JacobianMatrix[i].ColumnWiseArrayPosition] = + JacobeanColumnWise[JacobianMatrix[i].ColumnWiseArrayPosition] + temp; + } + + for (i = 0; i <= (int) JacobeanColumnWise.size() ; i++) + { + a[i] = JacobeanColumnWise[i]; + } + + /* + * Debug Output + */ + + /* + cout << "\n" << "\n"; + int temp_output = (int) sqrt(JacobeanColumnWise.size()); + for(i=0;i< temp_output;i++) + { + for(j=0;j<(int) JacobeanColumnWise.size();j++) + { + if(j % temp_output == i) + { + cout << a[j] << " "; + //cout << JacobeanColumnWise[j] << " "; + //cout << j % temp_output << " "; + } + } + cout << "\n"; + } + cout << "\n" << "\n"; + //*/ + //cout << "check\n"; +} + + diff --git a/source/get_input/Get_Species.cpp b/source/get_input/Get_Species.cpp index 7eade6d..7a2cc1d 100644 --- a/source/get_input/Get_Species.cpp +++ b/source/get_input/Get_Species.cpp @@ -17,13 +17,6 @@ vector< string > Get_Species(string filename) while (Mechanism_Data.good()) { - // http://www.cplusplus.com/forum/general/51349/ - //#ifdef Win32 - // getline (Mechanism_Data,line,r); - //else - // getline (Mechanism_Data,line,n); - //#endif - getline (Mechanism_Data,line1); if(spec_flag && !end_flag) @@ -38,22 +31,6 @@ vector< string > Get_Species(string filename) if(!end_flag){ - /* - // http://stackoverflow.com/questions/236129/splitting-a-string-in-c - istringstream iss(line1); - - do - { - string sub; - iss >> sub; - if(!sub.empty()){ // avoid empty lines - temp_species.push_back(sub); - //printf("%s\r\n",sub.c_str()); // Species read in correctly - } - }while (iss);//*/ - - - char *cstr, *p; vector< string > SplitLine; @@ -83,9 +60,9 @@ vector< string > Get_Species(string filename) temp_species.push_back(p); p=strtok(NULL," "); } - delete[] cstr; - + delete[] cstr; + delete[] p; } } @@ -93,7 +70,6 @@ vector< string > Get_Species(string filename) found = line1.find("SPEC"); if (found!=string::npos && !end_flag) { - //printf("Species found!\n\r"); spec_flag=1; } diff --git a/source/get_input/Handle_Mechanism_Input.cpp b/source/get_input/Handle_Mechanism_Input.cpp index 4ed09d8..401563b 100644 --- a/source/get_input/Handle_Mechanism_Input.cpp +++ b/source/get_input/Handle_Mechanism_Input.cpp @@ -16,7 +16,8 @@ * 5) Extras, e.g. PetroOxyData */ -void Handle_Mechanism_Input( +//void Handle_Mechanism_Input( +bool Handle_Mechanism_Input( string Filename_Mechanism, vector< string >& Species, vector< ThermodynamicData >& Thermodynamics, @@ -42,13 +43,13 @@ void Handle_Mechanism_Input( // check the existence of the 1st input file - the mechanism DataInputFromFile.open(MechanismData.c_str()); if (!DataInputFromFile.is_open()) { - printf("Error opening chem.inp - testing mechanism_data \n"); + cout << "Error opening chem.inp - testing mechanism_data \n"; // for the chemkin users, let the input file be named chem.inp MechanismData = "mechanism_data"; DataInputFromFile.open(MechanismData.c_str()); if (!DataInputFromFile.is_open()) { - printf("Error opening either chem.inp or mechanism_data \n"); - //return -1; + cout << "Error opening either chem.inp or mechanism_data \n"; + return false; } } DataInputFromFile.close(); @@ -57,8 +58,8 @@ void Handle_Mechanism_Input( // check the existence of the 2nd input file - the input data DataInputFromFile.open(InputData.c_str()); if (!DataInputFromFile.is_open()) { - printf("Error opening initial.inp \n"); - //return -1; + cout << "Error opening initial.inp \n"; + return false; } DataInputFromFile.close(); @@ -66,23 +67,24 @@ void Handle_Mechanism_Input( // As we now know that the input files exist, let us continue by reading in // the species list, thermodynamic data and mechanism + + // Get and store Species Information Species = Get_Species(MechanismData); Number_Species = Species.size(); - printf("The Input File contains %u Species. \n", Number_Species); - + cout << "The Input File contains " << Number_Species <<" Species.\n"; Thermodynamics = Get_Thermodynamic_Data_New_Format(MechanismData, Species); - printf("The Input File contains %u Thermodynamic Data Entries. \n", (int) Thermodynamics.size()); - - + cout << "The Input File contains " << Thermodynamics.size() + << " Thermodynamic Data Entries.\n"; // Get and store the Reaction Mechanism data Reactions = Read_Reaction_Matrix(MechanismData, Species); - Number_Reactions = Reactions.size(); - printf("The Input File contains %u Reactions. \n", Number_Reactions); + cout << "The Input File contains " << Number_Reactions << " Reactions.\n"; + + // did the user request species killing? - If yes, execute DataInputFromFile.open("kill.txt"); @@ -144,8 +146,9 @@ void Handle_Mechanism_Input( InitalSpecies ); // new function for improved input reading - printf("Initial concentrations are supplied for %u species as follow: \n", - (int) InitalSpecies.size());//InitialParameters.nrspec); + + cout << "Initial concentrations are supplied for " << InitalSpecies.size() + << " species as follow:\n"; /* @@ -237,13 +240,13 @@ void Handle_Mechanism_Input( // and the Initial Temperature InitialSpeciesConcentration[Number_Species] = InitialParameters.temperature; //Input[1][0]; - printf("\n"); + cout << "\n"; // Did the user request an irreversible scheme? if (InitialParameters.irrev) // contains true of false { - printf("Transformation to irreversible scheme requested! \n"); + cout << "Transformation to irreversible scheme requested!\n"; Reactions = Make_Irreversible(Reactions, Thermodynamics); WriteReactions("irreversible_scheme.txt", Species, Reactions); } @@ -278,7 +281,7 @@ void Handle_Mechanism_Input( /* Handle Input In here if it exists */ SpeciesMapping = Get_Combine_Species_Mapping("species_mapping.txt"); - printf("Species Mapping Information Read in. \n"); + cout << "Species Mapping Information Read in. \n"; /* * By supplying the mapping, the user has to agreed to try and combine the species. @@ -322,29 +325,16 @@ void Handle_Mechanism_Input( InitialSpeciesConcentration.clear(); // we overwrite the reactions, so why not the initial concentrations, can be used later for reset InitialSpeciesConcentration = temp; - printf("Reduced Reaction matrix consists of %i reactions and %i species. \n", - (int) Reactions.size(), (int) Reactions[0].Reactants.size()); + cout << "Reduced Reaction matrix consists of " << Reactions.size() + << " reactions and " << Reactions[0].Reactants.size() << " species.\n"; } DataInputFromFile.close(); - // did the user request a PetroOxy modification? - /* Data is: - * 0) Sample Size <- m^3 - * p - 1) Initial Pressure - * T - 2) Target Temperature - * pmax - 3) Maximum Pressure - * 4) O2 space in PetroOxy <- m^3 - * 5) Gas Species - * 6) mol of gas species at 25 degree celsius - * 7) O2 derived pressure - * 8) Vapour Pressure solvent component - * 9) solubility of gas at 298K, 1atm mol/L - * 10) "k" as Henry's Law Constant, M / L - */ - /* New PetroOxy handling */ + + // did the user request a PetroOxy modification? if(InitialParameters.PetroOxy) { @@ -353,12 +343,6 @@ void Handle_Mechanism_Input( double R = 8.314472e0; //double Na = 6.0221415e23; - //PetroOxyData[0] = InitialParameters.PetroOxySampleSize; - //PetroOxyData[1] = InitialParameters.PetroOxyInitPressure; - //PetroOxyData[2] = InitialParameters.temperature; - //PetroOxyData[3] = InitialParameters.PetroOxyMaxPressure; - //PetroOxyData[5] = InitialParameters.PetroOxyGasSpecies; - //PetroOxyData[9] = InitialParameters.PetroOxyGasSolubility; PetroOxyData.SampleSize = InitialParameters.PetroOxySampleSize; @@ -386,14 +370,12 @@ void Handle_Mechanism_Input( // Mod for limitied diffusion PetroOxyData.HenryLawDiffusionLimitSet = InitialParameters.HenryLawDiffusionLimitSet; PetroOxyData.HenryLawDiffusionLimit = InitialParameters.HenryLawDiffusionLimit; - //cout << PetroOxyData.HenryLawDiffusionLimitSet << "\n"; - //cout << PetroOxyData.HenryLawDiffusionLimit << "\n"; } - + return true; } diff --git a/source/get_input/Read_Input_Data_v2.cpp b/source/get_input/Read_Input_Data_v2.cpp index 845066d..6b8aad6 100644 --- a/source/get_input/Read_Input_Data_v2.cpp +++ b/source/get_input/Read_Input_Data_v2.cpp @@ -101,7 +101,6 @@ void Read_Input_Data_v2( string str; size_t found; - //printf("Line In: %s \n\n",line1.c_str()); // while the file is open while (Input_Data.good()) @@ -235,7 +234,6 @@ void Read_Input_Data_v2( line1.clear(); delete[] cstr; SetupParam.ReduceReactions = LineIn[0]; - //printf("Test - Temperature is: %.3e \n",LineIn[0]); } @@ -260,7 +258,6 @@ void Read_Input_Data_v2( line1.clear(); delete[] cstr; SetupParam.temperature = LineIn[0]; - //printf("Test - Temperature is: %.3e \n",LineIn[0]); } } @@ -283,7 +280,6 @@ void Read_Input_Data_v2( line1.clear(); delete[] cstr; //SetupParam.RatesAnalysisAtTimeData = LineIn; - //printf("Test - Temperature is: %.3e \n",LineIn[0]); } @@ -303,7 +299,6 @@ void Read_Input_Data_v2( line1.clear(); delete[] cstr; SetupParam.hm = LineIn[0]; - //printf("Test - Temperature is: %.3e \n",LineIn[0]); } found = line1.find("initialh"); @@ -322,7 +317,6 @@ void Read_Input_Data_v2( line1.clear(); delete[] cstr; SetupParam.h = LineIn[0]; - //printf("Test - Temperature is: %.3e \n",LineIn[0]); } @@ -343,7 +337,6 @@ void Read_Input_Data_v2( delete[] cstr; SetupParam.rtol = LineIn[0]; SetupParam.threshold = LineIn[1]; - //printf("Test - Tolerance is: %.3e and Threshhold is: %.3e \n",LineIn[0],LineIn[1]); } found = line1.find("Threshold"); @@ -361,7 +354,6 @@ void Read_Input_Data_v2( line1.clear(); delete[] cstr; SetupParam.threshold = LineIn[0]; - //printf("Test - Tolerance is: %.3e and Threshhold is: %.3e \n",LineIn[0],LineIn[1]); } found = line1.find("RTOL"); @@ -379,7 +371,6 @@ void Read_Input_Data_v2( line1.clear(); delete[] cstr; SetupParam.rtol = LineIn[0]; - //printf("Test - Tolerance is: %.3e and Threshhold is: %.3e \n",LineIn[0],LineIn[1]); } @@ -401,31 +392,9 @@ void Read_Input_Data_v2( // allows users to provide multiple time points SetupParam.TimeEnd.push_back(LineIn[0]); SetupParam.TimeStep.push_back(LineIn[1]); - //printf("Test - End Time is: %.3e and Time Step is; %.3e \n",LineIn[0],LineIn[1]); } - // Legacy function - /* - found = line1.find("TimeStepChange"); - if (found!=string::npos) - { - LineIn.clear(); // make sure storage array is empty - - p=strtok (cstr," \t"); // break at space or tab - p=strtok(NULL," \t"); // break again as first is the keyword - - while(p!=NULL){ // only read remainder is something is left - LineIn.push_back(strtod(p,NULL)); - p=strtok(NULL," \t"); - } - line1.clear(); - delete[] cstr; - // Not implemented with new system - //printf("Test - Time Step Change Time is: %.3e and New Time Step is; %.3e \n",LineIn[0],LineIn[1]); - }//*/ - - /* * Gas Phase Code Extension */ @@ -495,7 +464,6 @@ void Read_Input_Data_v2( while (p!=NULL) { Readout.push_back(p); - //cout << p << "\n"; p=strtok(NULL,"<>"); } delete[] cstr; @@ -522,7 +490,6 @@ void Read_Input_Data_v2( while (p!=NULL) { Readout.push_back(p); - //cout << p << "\n"; p=strtok(NULL,"<>"); } delete[] cstr; @@ -545,7 +512,6 @@ void Read_Input_Data_v2( while (p!=NULL) { Readout.push_back(p); - //cout << p << "\n"; p=strtok(NULL,"<>"); } delete[] cstr; @@ -568,7 +534,6 @@ void Read_Input_Data_v2( while (p!=NULL) { Readout.push_back(p); - //cout << p << "\n"; p=strtok(NULL,"<>"); } delete[] cstr; @@ -592,7 +557,6 @@ void Read_Input_Data_v2( while (p!=NULL) { Readout.push_back(p); - //cout << p << "\n"; p=strtok(NULL,"<>"); } delete[] cstr; @@ -600,12 +564,9 @@ void Read_Input_Data_v2( for(int i = 0;i<(int)Species.size();i++) { - //if(strcmp(species_name.c_str(),Comparator.c_str()) == 0){ if(strcmp(Readout[1].c_str(),Species[i].c_str()) == 0){ SetupParam.PetroOxyGasSpecies = i; OxyGasSpeciesDefined = true; - //PetroOxyData[5] = i; //strtod(Readout[1].c_str(),NULL); - //cout << i << "\n"; } } @@ -624,7 +585,6 @@ void Read_Input_Data_v2( while (p!=NULL) { Readout.push_back(p); - //cout << p << "\n"; p=strtok(NULL,"<>"); } delete[] cstr; @@ -673,7 +633,6 @@ void Read_Input_Data_v2( SetupParam.HenryLawDiffusionLimitSet = true; Readout.clear(); - //cout << "Henry Law Limiter is Set! \n"; // checked, works } diff --git a/source/get_input/Read_Reaction_Matrix.cpp b/source/get_input/Read_Reaction_Matrix.cpp index b07ddc3..b384ecd 100644 --- a/source/get_input/Read_Reaction_Matrix.cpp +++ b/source/get_input/Read_Reaction_Matrix.cpp @@ -19,6 +19,7 @@ // Pass vector by reference as input needs not be changed + struct SpeciesWithCoefficient { int SpeciesID; @@ -32,7 +33,7 @@ SpeciesWithCoefficient Return_Species_With_Coefficient(string , const vector< s vector< SingleReactionData > Read_Reaction_Matrix(string filename, const vector< string >& Species){ - //vector< vector< vector< double > > > Get_Reactions(string filename, const vector< string >& Species) + /* Struct Reaction Matrix * 1) Vector Ractants @@ -41,7 +42,6 @@ vector< SingleReactionData > Read_Reaction_Matrix(string filename, const vector< */ - vector< SingleReactionData > Reaction_Matrix; ifstream Mechanism_Data; diff --git a/source/get_input/Read_Thermodynamic_Data_New_Format.cpp b/source/get_input/Read_Thermodynamic_Data_New_Format.cpp index 7af8923..c87d505 100644 --- a/source/get_input/Read_Thermodynamic_Data_New_Format.cpp +++ b/source/get_input/Read_Thermodynamic_Data_New_Format.cpp @@ -23,12 +23,11 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( int Number_Species = Species.size(); int i; + // Vector to store all data as a struct (more efficient than vector< vector> > vector< ThermodynamicData > read_in_thermodynamics(Number_Species); - //vector< vector< double > > temp_thermodynamics_sorted; - //* if (Mechanism_Data.is_open()) { @@ -107,21 +106,18 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( // now the 3 temperature readings getline (Mechanism_Data,line1); - //char *cstr, *p; - //vector< string > temp_split_line; str = line1; cstr = new char [str.size()+1]; strcpy (cstr, str.c_str()); p=strtok (cstr," "); // split by space and tab - //cout << "Checpoint 1\n"; - while (p!=NULL) { temp_split_line.push_back(p); p=strtok(NULL," "); } + delete[] cstr; temp_read_in_single_species.TLow = strtod(temp_split_line[0].c_str(),NULL); @@ -143,6 +139,7 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( temp_split_line.push_back(p); p=strtok(NULL," "); } + delete[] cstr; temp_read_in_single_species.NasaLow1 = strtod(temp_split_line[0].c_str(),NULL); @@ -168,7 +165,9 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( temp_split_line.push_back(p); p=strtok(NULL," "); } + delete[] cstr; + delete[] p; temp_read_in_single_species.NasaHigh1 = strtod(temp_split_line[0].c_str(),NULL); temp_read_in_single_species.NasaHigh2 = strtod(temp_split_line[1].c_str(),NULL); @@ -184,20 +183,19 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( read_in_thermodynamics[temp_species_id] = temp_read_in_single_species; //cout << read_in_thermodynamics.size() << "\n"; //cout << temp_species_id << "\n"; + } } - //cout << "Checkpoint 1\n"; + // Check if the Thermodynamics Data Begins, if yes set a start flag found = line1.find("ThermData"); if (found!=string::npos && !end_flag) { begin_flag=true; - //cout << "check!\n"; } - //cout << "Checkpoint 2 \n"; } @@ -205,16 +203,14 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( } - //*/ - //cout << read_in_thermodynamics.size() << "\n"; + Mechanism_Data.open (filename.c_str()); // standard data structure if (Mechanism_Data.is_open()) { - //cout << "Getting Thermo \n"; string line1; size_t found; @@ -226,18 +222,14 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( getline (Mechanism_Data,line1); - //cout << "checkpoint 1 \n"; - if(begin_flag && !end_flag) { - //cout << "checkpoint 2"; found = line1.find("END"); // need to check for end in loop for 4 line blocks if (found!=string::npos && begin_flag) { end_flag = true; } - //cout << "Checkpoint \2 n"; if( !end_flag && @@ -360,47 +352,5 @@ vector< ThermodynamicData > Get_Thermodynamic_Data_New_Format( } - - - - - //cout << read_in_thermodynamics.size() << "\n"; - - - // now sort the data to fit the legacy input - - /* - for(i=0;i temp(17); - - temp[0] = read_in_thermodynamics[i].TLow; - temp[1] = read_in_thermodynamics[i].THigh; - temp[2] = read_in_thermodynamics[i].TChange; - - temp[3] = read_in_thermodynamics[i].NasaLow1; - temp[4] = read_in_thermodynamics[i].NasaLow2; - temp[5] = read_in_thermodynamics[i].NasaLow3; - temp[6] = read_in_thermodynamics[i].NasaLow4; - temp[7] = read_in_thermodynamics[i].NasaLow5; - temp[8] = read_in_thermodynamics[i].NasaLow6; - temp[9] = read_in_thermodynamics[i].NasaLow7; - - - temp[10] = read_in_thermodynamics[i].NasaHigh1; - temp[11] = read_in_thermodynamics[i].NasaHigh2; - temp[12] = read_in_thermodynamics[i].NasaHigh3; - temp[13] = read_in_thermodynamics[i].NasaHigh4; - temp[14] = read_in_thermodynamics[i].NasaHigh5; - temp[15] = read_in_thermodynamics[i].NasaHigh6; - temp[16] = read_in_thermodynamics[i].NasaHigh7; - - temp_thermodynamics_sorted.push_back(temp); - }//*/ - - - //return temp_thermodynamics_sorted; - - return read_in_thermodynamics; } diff --git a/source/mechanism_reduction/Combine_Species-k-based.cpp b/source/mechanism_reduction/Combine_Species-k-based.cpp index 96f7cab..e2490df 100644 --- a/source/mechanism_reduction/Combine_Species-k-based.cpp +++ b/source/mechanism_reduction/Combine_Species-k-based.cpp @@ -24,8 +24,6 @@ vector< SingleReactionData > Process_Species_Combination_Reactions_v2( * as well as a count of unmapped species and species classes * , we can begin to compact the reaction matrix */ - //vector< vector< vector< double > > > temp_reactions1; // output - //vector< vector< double > > SingleReactionData; vector< SingleReactionData > temp_reactions1; @@ -35,7 +33,7 @@ vector< SingleReactionData > Process_Species_Combination_Reactions_v2( vector< double > ProductData; // Product Information ProductData.resize(Number_Species_Classes); vector< double > ReactionData; // Reaction parameters such as Arrhenius parameters and whether irreversible or not - //ReactionData.resize(4); + vector< int > SpeciesClassesSize; SpeciesClassesSize.resize(Number_Species_Classes); @@ -47,7 +45,8 @@ vector< SingleReactionData > Process_Species_Combination_Reactions_v2( // Per reaction: - printf("%i reactions and %i species classes to process. \n",Number_Reactions,Number_Species_Classes); + cout << Number_Reactions << " reactions and " << Number_Species_Classes << " species classes to process.\n"; + for(i=0;i Process_Species_Combination_Reactions_v2( // Reaction Parameters Unaffected - //ReactionData = Reactions[i][2]; // reassemble output SingleReactionData.Reactants = ReactantData; @@ -189,8 +187,6 @@ vector< SingleReactionData > Process_Species_Combination_Reactions_v2( vector< double > ReactantData(temp_reactions3[i].Reactants); // Reactant Information vector< double > ProductData(temp_reactions3[i].Products); // Product Information - //vector< double > ReactionData; // Reaction parameters such as Arrhenius parameters and whether irreversible or not - //vector< vector< double > > SingleReactionData; // 1 reaction SingleReactionData SingleReaction; @@ -256,12 +252,6 @@ vector< SingleReactionData > Process_Species_Combination_Reactions_v2( - //ReactionData.push_back(Average_A); - //ReactionData.push_back(Average_n); - //ReactionData.push_back(Average_Ea); - //ReactionData.push_back(1); // Set reaction to irrvervesible which it shoud be anyway - //ReactionData.resize(4); - SingleReaction.Reactants = ReactantData; SingleReaction.Products = ProductData; SingleReaction.paramA = Average_A; @@ -274,8 +264,6 @@ vector< SingleReactionData > Process_Species_Combination_Reactions_v2( ReactantData.clear(); ProductData.clear(); - //ReactionData.clear(); - //SingleReactionData.clear(); // finally, remove the processed stuff: int delete_positions = Reaction_Grouping[i]; diff --git a/source/mechanism_reduction/Compare_Concentrations.cpp b/source/mechanism_reduction/Compare_Concentrations.cpp index b6cc4ec..3f460e9 100644 --- a/source/mechanism_reduction/Compare_Concentrations.cpp +++ b/source/mechanism_reduction/Compare_Concentrations.cpp @@ -92,116 +92,5 @@ vector< double > Compare_Concentrations( EvaluationResult[1] = large; } - - - /* - - vector< double > OldCheckPoints; - vector< double > NewCheckPoints; - - for(i=0;i 1.e-12 && gradient_new < gradient_old) // stopped getting steeper concentration less than 10^-12 is meaningless - { - OldCheckPoints.push_back(OldTimePoints[j]); - } - // negative gradient - if(gradient_new < -1.e-12 && gradient_new > gradient_old) // stopped getting steeper concentration less than 10^-12 is meaningless - { - OldCheckPoints.push_back(OldTimePoints[j]); - } - } - - for(j=1;j 1.e-12 && gradient_new < gradient_old) // stopped getting steeper concentration less than 10^-12 is meaningless - { - NewCheckPoints.push_back(NewTimePoints[j]); - } - // negative gradient - if(gradient_new < -1.e-12 && gradient_new > gradient_old) // stopped getting steeper concentration less than 10^-12 is meaningless - { - NewCheckPoints.push_back(NewTimePoints[j]); - } - } - - // OldTimePoints & NewTimePoints - - - printf("Species %u: ", i); - for(j=0;j<(int)OldCheckPoints.size();j++) - { - printf("%u ",OldCheckPoints[j]); - } - printf("\n"); - OldCheckPoints.clear(); - - printf("Species %u: ", i); - for(j=0;j<(int)NewCheckPoints.size();j++) - { - printf("%u ",NewCheckPoints[j]); - } - printf("\n"); - NewCheckPoints.clear(); - } -//*/ - - - - - //} - - /* - // now we need to compare the difference, why not an average? - //double PercentageSum = 0; - vector< double > EvaluationResult(2); - vector< double > Evaluation(Number_Species); - int large = 0; - for(i=0;i1.e-15) // only register discrepancies for relevant species - //{ - large = large + 1; - //} - } - } - //*/ - - - //EvaluationResult[0] = EvaluationResult[0] / (Number_Species-large); - //EvaluationResult = PercentageSum / (Number_Species+1); // dont't forget temperature.... -> so +1 - - - - //EvaluationResult[0] = 1; - //EvaluationResult[1] = 1; - return EvaluationResult; } diff --git a/source/mechanism_reduction/Get_Combine_Species_Mapping.cpp b/source/mechanism_reduction/Get_Combine_Species_Mapping.cpp index 1a4814a..a21f9f3 100644 --- a/source/mechanism_reduction/Get_Combine_Species_Mapping.cpp +++ b/source/mechanism_reduction/Get_Combine_Species_Mapping.cpp @@ -34,7 +34,7 @@ vector< vector< string > > Get_Combine_Species_Mapping(string filename) found = line.find("END"); if (found!=string::npos && continue_flag) { - printf("End found\n\r"); + //cout << "End found\n"; end_flag = 1; } @@ -92,7 +92,7 @@ vector< vector< string > > Get_Combine_Species_Mapping(string filename) found = line.find("MAPPING"); if (found!=string::npos) { - printf("Mapping found!\n\r"); + //cout << "Mapping found!\n"; continue_flag=1; } diff --git a/source/mechanism_reduction/Remove_Species.cpp b/source/mechanism_reduction/Remove_Species.cpp index 9b08feb..cf92201 100644 --- a/source/mechanism_reduction/Remove_Species.cpp +++ b/source/mechanism_reduction/Remove_Species.cpp @@ -35,7 +35,7 @@ vector< bool > Read_Kill_List( // first I need to read in my species to remove if (InputFile.is_open()) { - cout << "Let us start the cull! Killing selected species \n"; + cout << "Killing selected species\n"; while (InputFile.good()) { @@ -178,7 +178,7 @@ void Reduce_Species_Thermo_Mechanism( Thermodynamics.clear(); - cout << "The Scheme now contains the following counts \n" << + cout << "\nThe Scheme now contains the following counts\n" << "Species: " << NewSpecies.size() << "\n" << "Thermodynamic Entries: " << NewThermodynamics.size() << "\n" << "Reactions: " << NewReactions.size() << "\n\n"; diff --git a/source/mechanism_reduction/Species_Class_Mapping.cpp b/source/mechanism_reduction/Species_Class_Mapping.cpp index d2a4e54..5c98e55 100644 --- a/source/mechanism_reduction/Species_Class_Mapping.cpp +++ b/source/mechanism_reduction/Species_Class_Mapping.cpp @@ -34,7 +34,7 @@ vector< int > Map_Species_Classes( stringstream MappingParameter(SpeciesMapping[i][1]); if( (MappingParameter >> OrderedSpeciesMapping[j]).fail() ) { - printf("Failure in Species Ordering... \n"); + cout << "Failure in Species Ordering... \n"; } //OrderedSpeciesMapping[j] = strtod(SpeciesMapping[i][j].c_str(),NULL); // compare class ID to int @@ -63,9 +63,10 @@ vector< int > Map_Species_Classes( } } NumberOfClasses = CountUnmappedSpecies + BiggestClass; // sum of unmapped species and the number of classes - printf("The Mapping consists of %u unmapped Species, " - "and %u dedicated class mappings, resulting in %u Species Classes. \n" - ,CountUnmappedSpecies,BiggestClass,NumberOfClasses); + + cout << "The Mapping consists of " << CountUnmappedSpecies << " unmapped Species, and " + << BiggestClass << "dedicated class mappings, resulting in " << NumberOfClasses + << "Species Classes. \n"; // Classed species have a number, unclassed species get identified via a negative ID: // This will only work if I coded the rest well as position 0 drops out @@ -78,13 +79,11 @@ vector< int > Map_Species_Classes( OrderedSpeciesMapping[i] = -UnclassedSpeciesGrouping; UnclassedSpeciesGrouping = UnclassedSpeciesGrouping + 1; } - // printf("%i \n",OrderedSpeciesMapping[i]); } for(i=0;i& Kr, const double Temperature, const vector< ReactionParameter >& ReactionParameters, - const vector< vector< double > >& CalculatedThermo, + //const vector< vector< double > >& CalculatedThermo, + const vector< CalculatedThermodynamics >& CalculatedThermo, const vector< TrackSpecies >& SpeciesAll, const vector< double >& Delta_N ) @@ -62,12 +63,12 @@ void Calculate_Rate_Constant( delta_H[SpeciesAll[i].ReactionID] = delta_H[SpeciesAll[i].ReactionID] + - CalculatedThermo[SpeciesAll[i].SpeciesID][0] * + CalculatedThermo[SpeciesAll[i].SpeciesID].Hf * SpeciesAll[i].coefficient; delta_S[SpeciesAll[i].ReactionID] = delta_S[SpeciesAll[i].ReactionID] + - CalculatedThermo[SpeciesAll[i].SpeciesID][3] * + CalculatedThermo[SpeciesAll[i].SpeciesID].S * SpeciesAll[i].coefficient; } diff --git a/source/solver_calculations/Calculate_Thermodynamics.cpp b/source/solver_calculations/Calculate_Thermodynamics.cpp index 64b71e2..cc0d277 100644 --- a/source/solver_calculations/Calculate_Thermodynamics.cpp +++ b/source/solver_calculations/Calculate_Thermodynamics.cpp @@ -21,7 +21,8 @@ */ void Calculate_Thermodynamics( - vector< vector< double > >& CalculatedThermo, + //vector< vector< double > >& CalculatedThermo, + vector< CalculatedThermodynamics >& CalculatedThermo, const double& CurrentTemp, const vector< ThermodynamicData >& Thermodynamics ) @@ -51,7 +52,7 @@ void Calculate_Thermodynamics( //Hf - CalculatedThermo[i][0] = R*CurrentTemp*( + CalculatedThermo[i].Hf = R*CurrentTemp*( Thermodynamics[i].NasaLow1 + Thermodynamics[i].NasaLow2*CurrentTemp*0.5 + Thermodynamics[i].NasaLow3*CurrentTemp2/3 + @@ -60,7 +61,7 @@ void Calculate_Thermodynamics( Thermodynamics[i].NasaLow6/CurrentTemp); //Cp - CalculatedThermo[i][1] = R*( + CalculatedThermo[i].Cp = R*( Thermodynamics[i].NasaLow1+ Thermodynamics[i].NasaLow2*CurrentTemp + Thermodynamics[i].NasaLow3*CurrentTemp2 + @@ -68,10 +69,10 @@ void Calculate_Thermodynamics( Thermodynamics[i].NasaLow5*CurrentTemp4); //Cv = Cp - R; - CalculatedThermo[i][2] = CalculatedThermo[i][1] - R;//Cv; + CalculatedThermo[i].Cp = CalculatedThermo[i].Cp - R;//Cv; //S - CalculatedThermo[i][3] = R*( + CalculatedThermo[i].S = R*( Thermodynamics[i].NasaLow1*logCurrentTemp + Thermodynamics[i].NasaLow2*CurrentTemp + Thermodynamics[i].NasaLow3*CurrentTemp2*0.5 + @@ -83,7 +84,7 @@ void Calculate_Thermodynamics( else { //Hf = - CalculatedThermo[i][0] = R*CurrentTemp*( + CalculatedThermo[i].Hf = R*CurrentTemp*( Thermodynamics[i].NasaHigh1+ Thermodynamics[i].NasaHigh2*CurrentTemp*0.5 + Thermodynamics[i].NasaHigh3*CurrentTemp2/3 + @@ -92,7 +93,7 @@ void Calculate_Thermodynamics( Thermodynamics[i].NasaHigh6/CurrentTemp); //Cp - CalculatedThermo[i][1] = R*( + CalculatedThermo[i].Cp = R*( Thermodynamics[i].NasaHigh1+ Thermodynamics[i].NasaHigh2*CurrentTemp + Thermodynamics[i].NasaHigh3*CurrentTemp2 + @@ -100,10 +101,10 @@ void Calculate_Thermodynamics( Thermodynamics[i].NasaHigh5*CurrentTemp4); //Cv = Cp - R; - CalculatedThermo[i][2] = CalculatedThermo[i][1] - R; + CalculatedThermo[i].Cv = CalculatedThermo[i].Cp - R; //S - CalculatedThermo[i][3] = R*( + CalculatedThermo[i].S = R*( Thermodynamics[i].NasaHigh1*logCurrentTemp + Thermodynamics[i].NasaHigh2*CurrentTemp + Thermodynamics[i].NasaHigh3*CurrentTemp2*0.5 + diff --git a/source/write_output/WriteSpecies.cpp b/source/write_output/WriteSpecies.cpp index 5281db7..db2af4a 100644 --- a/source/write_output/WriteSpecies.cpp +++ b/source/write_output/WriteSpecies.cpp @@ -35,7 +35,7 @@ void WriteSpecies(string filename ,const vector< string >& Species) } else cout << "Unable to open file"; - printf("File %s written. \n",filename.c_str()); + cout << "File " << filename << " written.\n"; } diff --git a/source/write_output/Write_Mechanism.cpp b/source/write_output/Write_Mechanism.cpp index f512cc3..f13bcab 100644 --- a/source/write_output/Write_Mechanism.cpp +++ b/source/write_output/Write_Mechanism.cpp @@ -43,7 +43,7 @@ void WriteReactions( } else cout << "Unable to open file"; - printf("File %s written. \n",filename.c_str()); + cout << "File " << filename << " written.\n"; } diff --git a/source/write_output/Write_Stoichiometry_Matrix_For_Opt.cpp b/source/write_output/Write_Stoichiometry_Matrix_For_Opt.cpp index 7fc2c36..99d0ff8 100644 --- a/source/write_output/Write_Stoichiometry_Matrix_For_Opt.cpp +++ b/source/write_output/Write_Stoichiometry_Matrix_For_Opt.cpp @@ -10,7 +10,7 @@ void Write_Stoichiometric_Matrix_For_Opt ( string filename , - const vector< string >& Species, + //const vector< string >& Species, const vector< SingleReactionData >& Reactions ) {