From 2b117727777211e3a97c1c16e231ef59b5fec58b Mon Sep 17 00:00:00 2001 From: David L Parkhurst Date: Fri, 10 Jun 2011 22:00:46 +0000 Subject: [PATCH] Rewrote run_cells to do all of the calculations that would be done if a series of USE and SAVE. Modified the logic when time step is defined in run_cells. Equivalent to time_step in nmax steps, where nmax is max steps in kinetics, reaction, reaction_temperature. Wrote dump for chart and chart_handler. chart_handler calls dump for all user_graph. Fixed bug in chart_object with log definition. Had wrong logic when looking for t/T/l/L. initial_total_time initialized twice in mainsubs and twice in phreeqc.cpp. Fixed heading logic for mixing graph_x, graph_y, graph_sy with plot_xy. Moved NA to header file. git-svn-id: svn://136.177.114.72/svn_GW/phreeqcpp/trunk@5431 1feff8c3-07ed-0310-ac33-dd36852eb9cd --- ChartHandler.cpp | 11 +++ ChartHandler.h | 2 +- ChartObject.cpp | 180 +++++++++++++++++++++++++++++++++++++++- ChartObject.h | 5 ++ Phreeqc.cpp | 1 - ReadClass.cxx | 209 +++++++++++++++++++++++++++++++++++++++++++++++ runner.cpp | 28 ++++--- runner.h | 3 + 8 files changed, 425 insertions(+), 14 deletions(-) diff --git a/ChartHandler.cpp b/ChartHandler.cpp index 88d9362e..36264e90 100644 --- a/ChartHandler.cpp +++ b/ChartHandler.cpp @@ -162,5 +162,16 @@ ChartHandler::End_timer(PHREEQC_PTR_ARG) return true; } +bool +ChartHandler::dump(std::ostream & oss, unsigned int indent) +{ + std::map::iterator it = this->chart_map.begin(); + for ( ; it != chart_map.end(); it++) + { + size_t i = 0; + it->second->dump(oss, indent); + } + return true; +} #endif diff --git a/ChartHandler.h b/ChartHandler.h index 40eb295e..14bca879 100644 --- a/ChartHandler.h +++ b/ChartHandler.h @@ -29,7 +29,7 @@ public: bool Read(PHREEQC_PTR_ARG_COMMA CParser &parser); void Punch_user_graph(PHREEQC_PTR_ARG); bool End_timer(PHREEQC_PTR_ARG); - + bool dump(std::ostream & oss, unsigned int indent); protected: std::map chart_map; int current_chart_n_user; diff --git a/ChartObject.cpp b/ChartObject.cpp index d776b15d..5a7a8823 100644 --- a/ChartObject.cpp +++ b/ChartObject.cpp @@ -205,7 +205,7 @@ ChartObject::Set_axis_scale(CParser & parser) if (string_vector.size() == 6) { std::string s = string_vector[5]; - if (s[0] != 't' || s[0] != 'T' || s[0] != 'l' || s[0] != 'L') + if (s[0] == 't' || s[0] == 'T' || s[0] == 'l' || s[0] == 'L') { scale_ptr[4] = 10.0; if (((fabs(scale_ptr[0] - NA) > 1e-3) && scale_ptr[0] <=0) || @@ -330,7 +330,7 @@ ChartObject::Read(CParser & parser) { this->new_headings.push_back(token); } - + this->headings_original = this->new_headings; break; case 4: /* chart title */ { @@ -407,6 +407,7 @@ ChartObject::Read(CParser & parser) parser.get_rest_of_line(file_name); file_name = trim(file_name); this->OpenCSVFile(file_name); + this->csv_file_names.push_back(file_name); } break; case 13: /* clear */ @@ -423,11 +424,13 @@ ChartObject::Read(CParser & parser) if (!new_command_lines) { this->rate_command_list.clear(); + this->rate_command_list_original.clear(); new_command_lines = true; } this->rate_new_def = true; /* read command */ std::string cmd(parser.line()); + this->rate_command_list_original.push_back(cmd); std::string cmd_lower = cmd; Utilities::str_tolower(cmd_lower); if ((cmd_lower.find("graph_y") != std::string::npos) || @@ -1097,4 +1100,177 @@ ChartObject::Set_rate_new_def(bool tf) } } } +void +ChartObject::dump(std::ostream & oss, unsigned int indent) +{ + size_t i; + oss.precision(DBL_DIG - 1); + std::string indent0(""), indent1(""); + for (i = 0; i < indent; ++i) + indent0.append(Utilities::INDENT); + for (i = 0; i < indent + 1; ++i) + indent1.append(Utilities::INDENT); + oss << indent0 << "USER_GRAPH " << this->n_user << " " << this->description << std::endl; + + // chart title + oss << indent1 << "-chart_title \"" << this->chart_title << "\"" << std::endl; + + // axis titles + if (this->axis_titles.size() > 0) + { + oss << indent1 << "-axis_titles "; + + for (i = 0; i < this->axis_titles.size(); i++) + { + oss << "\"" << axis_titles[i] << "\" "; + } + oss << std::endl; + } + + // axis_scale_x + double *scale_ptr = this->axis_scale_x; + { + oss << indent1 << "-axis_scale x_axis "; + for (i = 0; i < 4; i++) + { + if (scale_ptr[i] == NA) + { + oss << " auto"; + } + else + { + oss << " " << scale_ptr[i]; + } + } + if (scale_ptr[4] == 10.0) + { + oss << " log"; + } + oss << std::endl; + } + // axis_scale_y + scale_ptr = this->axis_scale_y; + { + oss << indent1 << "-axis_scale y_axis "; + for (i = 0; i < 4; i++) + { + if (scale_ptr[i] == NA) + { + oss << " auto"; + } + else + { + oss << " " << scale_ptr[i]; + } + } + if (scale_ptr[4] == 10.0) + { + oss << " log"; + } + oss << std::endl; + } + // axis_scale_sy + scale_ptr = this->axis_scale_y2; + { + oss << indent1 << "-axis_scale sy_axis "; + for (i = 0; i < 4; i++) + { + if (scale_ptr[i] == NA) + { + oss << " auto"; + } + else + { + oss << " " << scale_ptr[i]; + } + } + if (scale_ptr[4] == 10.0) + { + oss << " log"; + } + oss << std::endl; + } + // chart type + if (this->chart_type == 0) + { + oss << indent1 << "-plot_concentration_vs x" << std::endl; + } + else + { + oss << indent1 << "-plot_concentration_vs t" << std::endl; + } + // graph_initial_solutions + if (this->graph_initial_solutions) + { + oss << indent1 << "-initial_solutions true" << std::endl; + } + else + { + oss << indent1 << "-initial_solutions false" << std::endl; + } + // connect_simulations + if (this->connect_simulations) + { + oss << indent1 << "-connect_simulations true" << std::endl; + } + else + { + oss << indent1 << "-connect_simulations false" << std::endl; + } + // csv files + for (i = 0; i < this->csv_file_names.size(); i++) + { + oss << indent1 << "-plot_csv_file " << this->csv_file_names[i] << std::endl; + } + + // headings + if (this->headings_original.size() > 0) + { + oss << indent1 << "-headings "; + for (i = 0; i < this->headings_original.size(); i++) + { + oss << this->headings_original[i] << " "; + } + oss << std::endl; + } + + // commands + oss << indent1 << "-start" << std::endl; + std::list::iterator it = rate_command_list_original.begin(); + for (; it != rate_command_list_original.end(); it++) + { + oss << *it << std::endl; + } + oss << indent1 << "-end" << std::endl; + + /* + struct rate *user_graph; + // C++ for rate struct + std::string rate_name; + std::list rate_command_list; + bool rate_new_def; + + int default_symbol; + int default_symbol_csv; + int default_color; + int default_color_csv; + + // temporary storage before stored graph_x/y/sy data are stored in curves + // Initialize_graph_pts and Finalize_graph_pts use this storage. + double graph_x; + std::map graph_y; + std::map secondary_y; + + // temporary plotxy curve definitions before stored in curves + // a plotxy curve is copied to Curves when cmdplotxy is encountered + // this keeps order correct between plotxy and graph_x/y/sy + std::vector new_plotxy_curves; + + // temporary headings until stored during basic_run + std::vector new_headings; + bool active; + bool detach; + bool form_started; + */ +} #endif // MULTICHART \ No newline at end of file diff --git a/ChartObject.h b/ChartObject.h index c890b759..99d64554 100644 --- a/ChartObject.h +++ b/ChartObject.h @@ -309,6 +309,7 @@ class ChartObject:public cxxNumKeyword double symbol_size = 6.0, int y_axis = 1, std::string color = ""); + void dump(std::ostream & s_oss, unsigned int indent); protected: bool new_ug; @@ -341,6 +342,7 @@ class ChartObject:public cxxNumKeyword bool end_timer; bool done; + std::vector csv_file_names; std::vector CurvesCSV; std::vector Curves; bool curve_added; @@ -350,6 +352,7 @@ class ChartObject:public cxxNumKeyword // C++ for rate struct std::string rate_name; std::list rate_command_list; + std::list rate_command_list_original; bool rate_new_def; int default_symbol; @@ -370,6 +373,7 @@ class ChartObject:public cxxNumKeyword // temporary headings until stored during basic_run std::vector new_headings; + std::vector headings_original; bool active; bool detach; bool form_started; @@ -377,6 +381,7 @@ class ChartObject:public cxxNumKeyword class Phreeqc * p_instance1; #endif + public: int usingResource; diff --git a/Phreeqc.cpp b/Phreeqc.cpp index 6f592cd9..7ae6e2c5 100644 --- a/Phreeqc.cpp +++ b/Phreeqc.cpp @@ -875,7 +875,6 @@ void Phreeqc::init(void) rate_sim_time_end = 0; rate_sim_time = 0; rate_moles = 0; - initial_total_time = 0; /* * user_print, user_punch diff --git a/ReadClass.cxx b/ReadClass.cxx index 7a89eab0..8a70cde2 100644 --- a/ReadClass.cxx +++ b/ReadClass.cxx @@ -2383,6 +2383,7 @@ run_as_cells(void) return (OK); } #endif +#ifdef SKIP2 /* ---------------------------------------------------------------------- */ int CLASS_QUALIFIER run_as_cells(void) @@ -2456,7 +2457,215 @@ run_as_cells(void) /* last_model.force_prep = TRUE; */ return (OK); } +#endif +/* ---------------------------------------------------------------------- */ +int CLASS_QUALIFIER +run_as_cells(void) +/* ---------------------------------------------------------------------- */ +{ + struct save save_data; + LDBLE kin_time; + int count_steps, use_mix, m; + char token[2 * MAX_LENGTH]; + struct kinetics *kinetics_ptr; + state = REACTION; + if (run_info.Get_cells().Get_numbers().size() == 0 || + !(run_info.Get_cells().Get_defined())) return(OK); + + // running cells + run_info.Set_run_cells(true); + + dup_print("Beginning of run as cells.", TRUE); + double initial_total_time_save; + if (run_info.Get_start_time() != NA) + { + initial_total_time_save = run_info.Get_start_time(); + } + else + { + initial_total_time_save = initial_total_time; + } + + std::set < int >::iterator it = run_info.Get_cells().Get_numbers().begin(); + + for ( ; it != run_info.Get_cells().Get_numbers().end(); it++) + { + int i = *it; + if (i < 0) continue; + initial_total_time = initial_total_time_save; + set_advection(i, TRUE, TRUE, i); +/* + * Run reaction step + */ + /* + * Find maximum number of steps + */ + dup_print("Beginning of batch-reaction calculations.", TRUE); + count_steps = 1; + if (use.irrev_in == TRUE && use.irrev_ptr != NULL) + { + if (abs(use.irrev_ptr->count_steps) > count_steps) + count_steps = abs(use.irrev_ptr->count_steps); + } + if (use.kinetics_in == TRUE && use.kinetics_ptr != NULL) + { + if (abs(use.kinetics_ptr->count_steps) > count_steps) + count_steps = abs(use.kinetics_ptr->count_steps); + } + if (use.temperature_in == TRUE && use.temperature_ptr != NULL) + { + if (abs(use.temperature_ptr->count_t) > count_steps) + count_steps = abs(use.temperature_ptr->count_t); + } + count_total_steps = count_steps; + /* + * save data for saving solutions + */ + memcpy(&save_data, &save, sizeof(struct save)); + /* + *Copy everything to -2 + */ + copy_use(-2); + rate_sim_time_start = 0; + rate_sim_time = 0; + for (reaction_step = 1; reaction_step <= count_steps; reaction_step++) + { + sprintf(token, "Reaction step %d.", reaction_step); + if (reaction_step > 1 && incremental_reactions == FALSE) + { + copy_use(-2); + } + set_initial_moles(-2); + dup_print(token, FALSE); + /* + * Determine time step for kinetics + */ + kin_time = 0.0; + if (use.kinetics_in == TRUE) + { + // runner kin_time + // equivalent to kin_time in count_steps + if (run_info.Get_time_step() != NA) + { + if (incremental_reactions == FALSE) + { + /* not incremental reactions */ + kin_time = reaction_step * run_info.Get_time_step() / ((LDBLE) count_steps); + } + else + { + /* incremental reactions */ + kin_time = run_info.Get_time_step() / ((LDBLE) count_steps); + } + } + // runner kin_time not defined + else + { + kinetics_ptr = kinetics_bsearch(-2, &m); + if (incremental_reactions == FALSE) + { + if (kinetics_ptr->count_steps > 0) + { + if (reaction_step > kinetics_ptr->count_steps) + { + kin_time = kinetics_ptr->steps[kinetics_ptr->count_steps - 1]; + } + else + { + kin_time = kinetics_ptr->steps[reaction_step - 1]; + } + } + else if (kinetics_ptr->count_steps < 0) + { + if (reaction_step > -kinetics_ptr->count_steps) + { + kin_time = kinetics_ptr->steps[0]; + } + else + { + kin_time = reaction_step * kinetics_ptr->steps[0] / ((LDBLE) (-kinetics_ptr->count_steps)); + } + } + } + else + { + /* incremental reactions */ + if (kinetics_ptr->count_steps > 0) + { + if (reaction_step > kinetics_ptr->count_steps) + { + kin_time = kinetics_ptr->steps[kinetics_ptr->count_steps - 1]; + } + else + { + kin_time = kinetics_ptr->steps[reaction_step - 1]; + } + } + else if (kinetics_ptr->count_steps < 0) + { + if (reaction_step > -kinetics_ptr->count_steps) + { + kin_time = 0; + } + else + { + kin_time = + kinetics_ptr->steps[0] / ((LDBLE) (-kinetics_ptr->count_steps)); + } + } + } + } + } + if (incremental_reactions == FALSE || + (incremental_reactions == TRUE && reaction_step == 1)) + { + use_mix = TRUE; + } + else + { + use_mix = FALSE; + } + /* + * Run reaction step + */ + run_reactions(-2, kin_time, use_mix, 1.0); + if (incremental_reactions == TRUE) + { + rate_sim_time_start += kin_time; + rate_sim_time = rate_sim_time_start; + } + else + { + rate_sim_time = kin_time; + } + if (state != ADVECTION) + { + punch_all(); + print_all(); + } + /* saves back into -2 */ + if (reaction_step < count_steps) + { + saver(); + } + } + /* + * save end of reaction + */ + memcpy(&save, &save_data, sizeof(struct save)); + if (use.kinetics_in == TRUE) + { + kinetics_duplicate(-2, use.n_kinetics_user); + } + saver(); + } + initial_total_time += rate_sim_time; + run_info.Get_cells().Set_defined(false); + // not running cells + run_info.Set_run_cells(false); + return (OK); +} /* ---------------------------------------------------------------------- */ void CLASS_QUALIFIER dump_ostream(std::ostream& os) diff --git a/runner.cpp b/runner.cpp index 1a8952bf..622e7063 100644 --- a/runner.cpp +++ b/runner.cpp @@ -1,14 +1,18 @@ #include "runner.h" #include "Parser.h" +#include "NA.h" runner::runner(void) { - this->time_step = 0; - this->start_time = 0; + this->time_step = NA; + this->start_time = NA; + this->run_cells = false; + } runner::runner(CParser & parser) { - this->time_step = 0; - this->start_time = 0; + this->time_step = NA; + this->start_time = NA; + this->run_cells = false; this->Read(parser); } @@ -27,6 +31,9 @@ bool runner::Read(CParser & parser) vopts.push_back("cells"); vopts.push_back("start_time"); vopts.push_back("time_step"); + vopts.push_back("time_steps"); + vopts.push_back("step"); + vopts.push_back("steps"); } std::istream::pos_type ptr; @@ -52,13 +59,11 @@ bool runner::Read(CParser & parser) opt_save = opt; } + //// Read dump entity list of numbers or number ranges for line, store in item + //if (opt >= 0 && opt <= 1) + //{ - - // Read dump entity list of numbers or number ranges for line, store in item - if (opt >= 0 && opt <= 1) - { - - } + //} // Process other identifiers std::set < int >::iterator it; @@ -99,6 +104,9 @@ bool runner::Read(CParser & parser) } break; case 3: //time_step + case 4: //time_steps + case 5: //step + case 6: //steps if (!(parser.get_iss() >> this->time_step)) { parser.error_msg("Expected time_step for RUN_CELLS.", diff --git a/runner.h b/runner.h index 11073164..033f659c 100644 --- a/runner.h +++ b/runner.h @@ -16,10 +16,13 @@ public: StorageBinListItem & Get_cells(void) { return(this->cells); }; double Get_time_step() { return(this->time_step); }; double Get_start_time() { return(this->start_time); }; + bool Get_run_cells() { return(this->run_cells); }; + void Set_run_cells(bool tf) { this->run_cells = tf; }; protected: double time_step; double start_time; StorageBinListItem cells; + bool run_cells; }; #endif // !defined(RUNNER_H_INCLUDED)