diff --git a/phreeqcpp/Makefile.old b/phreeqcpp/Makefile.old index 6d7279c2..b780ea89 100644 --- a/phreeqcpp/Makefile.old +++ b/phreeqcpp/Makefile.old @@ -158,7 +158,7 @@ ifeq ($(CFG), CLASS_DEBUG_64) CL1MP_OBJS=cl1mp.o CL1MP_LIB=libgmp.a endif - DEFINES = -DUSE_PHRQ_ALLOC $(DEFINE_INVERSE_CL1MP) # -DPHREEQC2 + DEFINES = -DUSE_PHRQ_ALLOC $(DEFINE_INVERSE_CL1MP) -DTEST_COPY_OPERATOR VPATH = ..:../PhreeqcKeywords:../common INCLUDES = -I.. -I../PhreeqcKeywords -I../common CXX = g++ diff --git a/phreeqcpp/Phreeqc.cpp b/phreeqcpp/Phreeqc.cpp index 6860488b..f7e499c3 100644 --- a/phreeqcpp/Phreeqc.cpp +++ b/phreeqcpp/Phreeqc.cpp @@ -277,7 +277,7 @@ size_t Phreeqc::list_SolidSolutions(std::list &list_comps, std::lis { std::string ssname = ssit->second.Get_name(); std::set accumulator_phases; - for (int i = 0; i < ssit->second.Get_ss_comps().size(); i++) + for (size_t i = 0; i < ssit->second.Get_ss_comps().size(); i++) { std::string pname = ssit->second.Get_ss_comps()[i].Get_name(); int j; @@ -348,7 +348,7 @@ size_t Phreeqc::list_Surfaces(std::list &list_surftype, std::listsecond; std::vector &scomps = entity.Get_surface_comps(); - // std::vector &scharges = entity.Get_surface_charges(); + //std::vector &scharges = entity.Get_surface_charges(); for (size_t i = 0; i < scomps.size(); i++) { std::pair p(scomps[i].Get_master_element(), scomps[i].Get_charge_name()); @@ -687,7 +687,7 @@ void Phreeqc::init(void) * Transport data *---------------------------------------------------------------------- */ count_cells = 1; - cell_data_max_cells = count_cells; + cell_data_max_cells = 1; // count_cells; count_shifts = 1; ishift = 1; bcon_first = bcon_last = 3; @@ -1063,7 +1063,8 @@ void Phreeqc::init(void) V_solutes = 0.0; viscos = 0.0; viscos_0 = 0.0; - rho_0 = 0; + viscos_0_25 = 0.0; + rho_0 = 0.0; kappa_0 = 0.0; p_sat = 0.0; eps_r = EPSILON; @@ -1395,9 +1396,10 @@ void Phreeqc::init(void) /* utilities.cpp ------------------------------- */ spinner = 0; // keycount; + keycount.resize(Keywords::KEY_COUNT_KEYWORDS); for (int i = 0; i < Keywords::KEY_COUNT_KEYWORDS; i++) { - keycount.push_back(0); + keycount[i] = 0; } return; @@ -1465,6 +1467,21 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) change_surf_count = 0; change_surf = NULL; */ + change_surf_count = pSrc->change_surf_count; + change_surf = change_surf_alloc(change_surf_count + 1); + if (change_surf_count > 0) + { + for (int ii = 0; ii < change_surf_count; ii++) + { + change_surf[ii].comp_name = string_hsave(pSrc->change_surf[ii].comp_name); + change_surf[ii].fraction = pSrc->change_surf[ii].fraction; + change_surf[ii].new_comp_name = string_hsave(pSrc->change_surf[ii].new_comp_name); + change_surf[ii].new_Dw = pSrc->change_surf[ii].new_Dw; + change_surf[ii].cell_no = pSrc->change_surf[ii].cell_no; + change_surf[ii].next = pSrc->change_surf[ii].next; + } + } + /* ---------------------------------------------------------------------- * Exchange * ---------------------------------------------------------------------- */ @@ -1502,19 +1519,31 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) *---------------------------------------------------------------------- */ // Should be empty after each END // auto Rxn_mix_map; + Rxn_mix_map = pSrc->Rxn_mix_map; // auto Dispersion_mix_map; + Dispersion_mix_map = pSrc->Dispersion_mix_map; // auto Rxn_solution_mix_map; + Rxn_solution_mix_map = pSrc->Rxn_solution_mix_map; // auto Rxn_exchange_mix_map; + Rxn_exchange_mix_map = pSrc->Rxn_exchange_mix_map; // auto Rxn_gas_phase_mix_map; + Rxn_gas_phase_mix_map = pSrc->Rxn_gas_phase_mix_map; // auto Rxn_kinetics_mix_map; + Rxn_kinetics_mix_map = pSrc->Rxn_kinetics_mix_map; // auto Rxn_pp_assemblage_mix_map; + Rxn_pp_assemblage_mix_map = pSrc->Rxn_pp_assemblage_mix_map; // auto Rxn_ss_assemblage_mix_map; + Rxn_ss_assemblage_mix_map = pSrc->Rxn_ss_assemblage_mix_map; // auto Rxn_surface_mix_map; + Rxn_surface_mix_map = pSrc->Rxn_surface_mix_map; + /* + * List new definitions + */ + // Assume no new definitions /*---------------------------------------------------------------------- * Irreversible reaction *---------------------------------------------------------------------- */ Rxn_reaction_map = pSrc->Rxn_reaction_map; - run_cells_one_step = pSrc->run_cells_one_step; /*---------------------------------------------------------------------- * Gas phase *---------------------------------------------------------------------- */ @@ -1570,6 +1599,9 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) *---------------------------------------------------------------------- */ /* title_x = NULL; + */ + last_title_x = pSrc->last_title_x; + /* new_x = FALSE; description_x = NULL; tc_x = 0; @@ -1609,7 +1641,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) * Transport data *---------------------------------------------------------------------- */ count_cells = pSrc->count_cells; - cell_data_max_cells = pSrc->cell_data_max_cells; + cell_data_max_cells = 1; //pSrc->cell_data_max_cells; count_shifts = pSrc->count_shifts; ishift = pSrc->ishift; bcon_first = pSrc->bcon_first; @@ -1623,6 +1655,8 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) cell = pSrc->cell; mcd_substeps = pSrc->mcd_substeps; /* stag_data */ + stag_data = (struct stag_data *) free_check_null(stag_data); + stag_data = (struct stag_data *) PHRQ_malloc(sizeof(struct stag_data)); memcpy(stag_data, pSrc->stag_data, sizeof(struct stag_data)); print_modulus = pSrc->print_modulus; punch_modulus = pSrc->punch_modulus; @@ -1630,19 +1664,31 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) dump_modulus = pSrc->dump_modulus; transport_warnings = pSrc->transport_warnings; /* cell_data */ - if (count_cells > 0) - { - cell_data = (struct cell_data *) free_check_null(cell_data); - cell_data = (struct cell_data *) PHRQ_malloc((size_t) ((count_cells + 2) * sizeof(struct cell_data))); - if (cell_data == NULL) malloc_error(); - memcpy(cell_data, pSrc->cell_data, ((size_t) ((count_cells + 2) * sizeof(struct cell_data)))); - } old_cells = pSrc->old_cells; max_cells = pSrc->max_cells; + + if (stag_data->count_stag > 0) + { + max_cells = (max_cells - 2) / (1 + stag_data->count_stag); + } + all_cells = pSrc->all_cells; + cell_data_max_cells = 1; + if (count_cells > 0) + { + //cell_data = (struct cell_data *) free_check_null(cell_data); + //cell_data = (struct cell_data *) PHRQ_malloc((size_t) ((count_cells + 2) * sizeof(struct cell_data))); + //if (cell_data == NULL) malloc_error(); + //memcpy(cell_data, pSrc->cell_data, ((size_t) ((count_cells + 2) * sizeof(struct cell_data + int all_cells_now = max_cells * (1 + stag_data->count_stag) + 2; + space((void **)((void *)&cell_data), all_cells_now, &cell_data_max_cells, sizeof(struct cell_data)); + memcpy(cell_data, pSrc->cell_data, ((size_t)(all_cells_now * sizeof(struct cell_data)))); + } + max_cells = pSrc->max_cells; multi_Dflag = pSrc->multi_Dflag; interlayer_Dflag = pSrc->interlayer_Dflag; default_Dw = pSrc->default_Dw; + correct_Dw = pSrc->correct_Dw; multi_Dpor = pSrc->multi_Dpor; interlayer_Dpor = pSrc->interlayer_Dpor; multi_Dpor_lim = pSrc->multi_Dpor_lim; @@ -1650,6 +1696,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) multi_Dn = pSrc->multi_Dn; interlayer_tortf = pSrc->interlayer_tortf; cell_no = pSrc->cell_no; + mixrun = pSrc->mixrun; fix_current = pSrc->fix_current; /*---------------------------------------------------------------------- * Advection data @@ -1700,6 +1747,12 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) /*---------------------------------------------------------------------- * Elements *---------------------------------------------------------------------- */ + //max_elements = pSrc->max_elements; + //elements = (struct element **) free_check_null(elements); + //elements = (struct element **) PHRQ_malloc((size_t)max_elements * sizeof(struct element)); + space((void **)((void *)&elements), pSrc->max_elements, &max_elements, + sizeof(struct element *)); + count_elements = 0; for (int i = 0; i < pSrc->count_elements; i++) { string_hsave(pSrc->elements[i]->name); @@ -1722,6 +1775,11 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) max_elts = MAX_ELTS; */ /*---------------------------------------------------------------------- + * Reaction + *---------------------------------------------------------------------- */ + //bool run_cells_one_step; + run_cells_one_step = pSrc->run_cells_one_step; + /*---------------------------------------------------------------------- * Species *---------------------------------------------------------------------- */ /* @@ -1746,6 +1804,15 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) s_o2 = NULL; */ // logk + + space((void **)((void *)&logk), pSrc->max_logk, &max_logk, sizeof(struct logk *)); + //for (int i = 0; i < count_logk; i++) + //{ + // logk[i] = (struct logk *) free_check_null(logk[i]); + //} + //logk = (struct logk **) free_check_null(logk); + //max_logk = pSrc->max_logk; + //logk = (struct logk **) PHRQ_malloc((size_t) max_logk * sizeof(struct logk *)); for (int i = 0; i < pSrc->count_logk; i++) { char * name = string_duplicate(pSrc->logk[i]->name); @@ -1756,6 +1823,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) logk_ptr->add_logk = NULL; if (logk_ptr->count_add_logk > 0) { + logk_ptr->add_logk = (struct name_coef *) free_check_null(logk_ptr->add_logk); logk_ptr->add_logk = (struct name_coef *) PHRQ_malloc((size_t) pSrc->logk[i]->count_add_logk * sizeof(struct name_coef)); if (logk[i]->add_logk == NULL) malloc_error(); for (int j = 0; j < logk_ptr->count_add_logk; j++) @@ -1765,7 +1833,15 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) } } } + count_logk = pSrc->count_logk; // s, species + count_s = 0; + //max_s = pSrc->max_s; + + //s = (struct species **) free_check_null(s); + //s = (struct species **) PHRQ_malloc(sizeof(struct species *)*size_t(max_s)); + + space((void **)((void *)&s), pSrc->max_s, &max_s, sizeof(struct species *)); for (int i = 0; i < pSrc->count_s; i++) { struct species *s_ptr = s_store(pSrc->s[i]->name, pSrc->s[i]->z, FALSE); @@ -1800,10 +1876,15 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) } //next_secondary s_ptr->next_secondary = NULL; - if (pSrc->s[i]->next_secondary) + if (pSrc->s[i]->next_secondary && pSrc->s[i]->mole_balance) { - cxxNameDouble next_secondary(pSrc->s[i]->next_secondary); - s_ptr->next_secondary = NameDouble2elt_list(next_secondary); + count_elts = 0; + paren_count = 0; + char * string = string_duplicate(s_ptr->mole_balance); + char * ptr = string; + get_secondary_in_species(&ptr, 1.0); + s_ptr->next_secondary = elt_list_save(); + free_check_null(string); } //next_sys_total s_ptr->next_sys_total = NULL; @@ -1850,6 +1931,11 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) count_phases = 0; max_phases = MAX_PHASES; */ + //max_phases = pSrc->max_phases; + //phases = (struct phase **) PHRQ_malloc((size_t)max_phases * sizeof(struct phase)); + //space((void **)((void *)&phases), INIT, &max_phases, + // sizeof(struct phase *)); + count_phases = 0; for (int i = 0; i < pSrc->count_phases; i++) { struct phase *phase_ptr = phase_store(pSrc->phases[i]->name); @@ -1891,14 +1977,14 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) phase_ptr->rxn = cxxChemRxn2rxn(rxn); } //rxn_s - phase_ptr->rxn_s = NULL; + //phase_ptr->rxn_s = NULL; if (pSrc->phases[i]->rxn_s != NULL) { cxxChemRxn rxn_s(pSrc->phases[i]->rxn_s); phase_ptr->rxn_s = cxxChemRxn2rxn(rxn_s); } //rxn_x - phase_ptr->rxn_x = NULL; + //phase_ptr->rxn_x = NULL; if (pSrc->phases[i]->rxn_x != NULL) { cxxChemRxn rxn_x(pSrc->phases[i]->rxn_x); @@ -1915,9 +2001,11 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) max_master = MAX_MASTER; */ count_master = pSrc->count_master; - max_master = pSrc->max_master; - master = (struct master **) free_check_null(master); - master = (struct master **) PHRQ_malloc((size_t) max_master * sizeof(struct master *)); + //max_master = pSrc->max_master; + //master = (struct master **) free_check_null(master); + //master = (struct master **) PHRQ_malloc((size_t) max_master * sizeof(struct master *)); + space((void **)((void *)&master), pSrc->max_master, &max_master, + sizeof(struct master *)); if (master == NULL) malloc_error(); dbg_master = master; for (int i = 0; i < count_master; i++) @@ -2054,6 +2142,9 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) rate_sim_time = 0; rate_moles = 0; initial_total_time = 0; + */ + initial_total_time = pSrc->initial_total_time; + /* // auto rate_p count_rate_p = 0; */ @@ -2076,6 +2167,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) /* ---------------------------------------------------------------------- * USER PRINT COMMANDS * ---------------------------------------------------------------------- */ + /* user_print = NULL; */ @@ -2144,6 +2236,9 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) /* error_string = NULL; simulation = 0; + */ + simulation = pSrc->simulation; + /* int state = INITIALIZE; reaction_step = 0; transport_step = 0; @@ -2151,6 +2246,9 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) advection_step = 0; stop_program = FALSE; incremental_reactions = FALSE; + */ + incremental_reactions = pSrc->incremental_reactions; + /* count_strings = 0; array = NULL; delta = NULL; @@ -2172,6 +2270,11 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) debug_diffuse_layer = FALSE; debug_inverse = FALSE; */ + debug_model = pSrc->debug_model; + debug_prep = pSrc->debug_prep; + debug_set = pSrc->debug_set; + debug_diffuse_layer = pSrc->debug_diffuse_layer; + debug_inverse = pSrc->debug_inverse; inv_tol_default = pSrc->inv_tol_default; itmax = pSrc->itmax; max_tries = pSrc->max_tries; @@ -2191,8 +2294,8 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) censor = pSrc->censor; aqueous_only = pSrc->aqueous_only; negative_concentrations = pSrc->negative_concentrations; - calculating_deriv = FALSE; - numerical_deriv = FALSE; + calculating_deriv = pSrc->calculating_deriv; + numerical_deriv = pSrc->numerical_deriv; count_total_steps = 0; phast = FALSE; /* @@ -2249,8 +2352,29 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) } // Not implemented for now - //SelectedOutput_map = pSrc->SelectedOutput_map; - SelectedOutput_map.clear(); + SelectedOutput_map = pSrc->SelectedOutput_map; + { + std::map::iterator it = SelectedOutput_map.begin(); + for (; it != SelectedOutput_map.end(); it++) + { + //phrq_io->punch_open(it->second.Get_file_name().c_str()); + //it->second.Set_punch_ostream(phrq_io->Get_punch_ostream()); + //phrq_io->Set_punch_ostream(NULL); + it->second.Set_punch_ostream(NULL); + } + } + //SelectedOutput_map.clear(); + + UserPunch_map = pSrc->UserPunch_map; + { + std::map::iterator it = UserPunch_map.begin(); + for (; it != UserPunch_map.end(); it++) + { + struct rate *rate_new = rate_copy(it->second.Get_rate()); + it->second.Set_rate(rate_new); + } + } + //selected_output_file_name = NULL; //dump_file_name = NULL; @@ -2280,6 +2404,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) master_isotope = NULL; max_master_isotope = MAX_ELTS; */ + for (int i = 0; i < pSrc->count_master_isotope; i++) { struct master_isotope *master_isotope_ptr = master_isotope_store(pSrc->master_isotope[i]->name, FALSE); @@ -2318,18 +2443,18 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) for (int i = 0; i < pSrc->count_calculate_value; i++) { struct calculate_value *calculate_value_ptr = calculate_value_store(pSrc->calculate_value[i]->name, FALSE); - memcpy(calculate_value_ptr, pSrc->calculate_value[i], sizeof(struct calculate_value)); + //memcpy(calculate_value_ptr, pSrc->calculate_value[i], sizeof(struct calculate_value)); calculate_value_ptr->value = pSrc->calculate_value[i]->value; - calculate_value_ptr->commands = NULL; + //calculate_value_ptr->commands = NULL; if (pSrc->calculate_value[i]->commands) { calculate_value_ptr->commands = string_duplicate(pSrc->calculate_value[i]->commands); } - calculate_value_ptr->new_def = TRUE; - calculate_value_ptr->calculated = FALSE; - calculate_value_ptr->linebase = NULL; - calculate_value_ptr->varbase = NULL; - calculate_value_ptr->loopbase = NULL; + //calculate_value_ptr->new_def = TRUE; + //calculate_value_ptr->calculated = FALSE; + //calculate_value_ptr->linebase = NULL; + //calculate_value_ptr->varbase = NULL; + //calculate_value_ptr->loopbase = NULL; } /* count_isotope_ratio = 0; @@ -2363,196 +2488,140 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) user_database = string_duplicate(pSrc->user_database); //have_punch_name = pSrc->have_punch_name; print_density = pSrc->print_density; + print_viscosity = pSrc->print_viscosity; #ifdef SKIP - zeros = NULL; - zeros_max = 1; - cell_pore_volume = 0; - cell_volume = 0; - cell_porosity = 0; - cell_saturation = 0; - sys = NULL; - count_sys = 0; - max_sys = 0; - sys_tot = 0; + LDBLE *zeros; + int zeros_max; +#endif + viscos = pSrc->viscos; + viscos_0 = pSrc->viscos_0; + viscos_0_25 = pSrc->viscos_0_25; // viscosity of the solution, of pure water, of pure water at 25 C +#ifdef SKIP + LDBLE cell_pore_volume; + LDBLE cell_porosity; + LDBLE cell_volume; + LDBLE cell_saturation; + struct system_species *sys; + int count_sys, max_sys; + LDBLE sys_tot; - V_solutes = 0.0; - rho_0 = 0; - kappa_0 = 0.0; - p_sat = 0.0; - eps_r = EPSILON; - DH_A = 0.0; - DH_B = 0.0; - DH_Av = 0.0; - QBrn = 0.0; - ZBrn = 0.0; - dgdP = 0.0; + LDBLE V_solutes, rho_0, rho_0_sat, kappa_0, p_sat/*, ah2o_x0*/; + LDBLE SC; // specific conductance mS/cm + LDBLE eps_r; // relative dielectric permittivity + LDBLE DH_A, DH_B, DH_Av; // Debye-Hueckel A, B and Av + LDBLE QBrn; // Born function d(ln(eps_r))/dP / eps_r * 41.84004, for supcrt calc'n of molal volume + LDBLE ZBrn; // Born function (-1/eps_r + 1) * 41.84004, for supcrt calc'n of molal volume + LDBLE dgdP; // dg / dP, pressure derivative of g-function, for supcrt calc'n of molal volume + + int need_temp_msg; + LDBLE solution_mass, solution_volume; - need_temp_msg = 0; - solution_mass = 0; - solution_volume = 0; /* phqalloc.cpp ------------------------------- */ - s_pTail = NULL; + PHRQMemHeader *s_pTail; + /* Basic */ - basic_interpreter = NULL; + PBasic * basic_interpreter; + double(*basic_callback_ptr) (double x1, double x2, const char *str, void *cookie); + void *basic_callback_cookie; +#ifdef IPHREEQC_NO_FORTRAN_MODULE + double(*basic_fortran_callback_ptr) (double *x1, double *x2, char *str, size_t l); +#else + double(*basic_fortran_callback_ptr) (double *x1, double *x2, const char *str, int l); +#endif +#if defined(SWIG) || defined(SWIG_IPHREEQC) + class BasicCallback *basicCallback; + void SetCallback(BasicCallback *cb) { basicCallback = cb; } +#endif + /* cl1.cpp ------------------------------- */ - x_arg = NULL; - res_arg = NULL; - scratch = NULL; - x_arg_max = 0; - res_arg_max = 0; - scratch_max = 0; + LDBLE *x_arg, *res_arg, *scratch; + int x_arg_max, res_arg_max, scratch_max; +#ifdef SKIP /* dw.cpp ------------------------------- */ - /* COMMON /QQQQ/ */ - Q0 = 0; - Q5 = 0; - GASCON = 0.461522e0; - TZ = 647.073e0; - AA = 1.e0; - Z = 0; - DZ = 0; - Y = 0; - G1 = 11.e0; - G2 = 44.333333333333e0; - GF = 3.5e0; - B1 = 0; - B2 = 0; - B1T = 0; - B2T = 0; - B1TT = 0; - B2TT = 0; + /* COMMON /QQQQ/ */ + LDBLE Q0, Q5; + LDBLE GASCON, TZ, AA; + LDBLE Z, DZ, Y; + LDBLE G1, G2, GF; + LDBLE B1, B2, B1T, B2T, B1TT, B2TT; +#endif /* gases.cpp ------------------------------- */ - a_aa_sum = 0; - b2 = 0; - b_sum = 0; - R_TK = 0; + LDBLE a_aa_sum, b2, b_sum, R_TK; + /* input.cpp ------------------------------- */ - check_line_return = 0; - reading_db = FALSE; + int check_line_return; + int reading_db; + /* integrate.cpp ------------------------------- */ - midpoint_sv = 0; - z_global = 0; - xd_global = 0; - alpha_global = 0; + LDBLE midpoint_sv; + LDBLE z_global, xd_global, alpha_global; + /* inverse.cpp ------------------------------- */ - max_row_count = 50; - max_column_count = 50; - carbon = FALSE; - col_name = NULL; - row_name = NULL; - count_rows = 0; - count_optimize = 0; - col_phases = 0; - col_redox = 0; - col_epsilon = 0; - col_ph = 0; - col_water = 0; - col_isotopes = 0; - col_phase_isotopes = 0; - row_mb = 0; - row_fract = 0; - row_charge = 0; - row_carbon = 0; - row_isotopes = 0; - row_epsilon = 0; - row_isotope_epsilon = 0; - row_water = 0; - inv_zero = NULL; - array1 = 0; - inv_res = NULL; - inv_delta1 = NULL; - delta2 = NULL; - delta3 = NULL; - inv_cu = NULL; - delta_save = NULL; - min_delta = NULL; - max_delta = NULL; - inv_iu = NULL; - inv_is = NULL; - klmd = 0; - nklmd = 0; - n2d = 0; - kode = 0; - iter = 0; - toler = 0; - error = 0; - max_pct = 0; - scaled_error = 0; - master_alk = NULL; - row_back = NULL; - col_back = NULL; - good = NULL; - bad = NULL; - minimal = NULL; - max_good = 0; - max_bad = 0; - max_minimal = 0; - count_good = 0; - count_bad = 0; - count_minimal = 0; - count_calls = 0; - soln_bits = 0; - phase_bits = 0; - current_bits = 0; - temp_bits = 0; - netpath_file = NULL; - count_inverse_models = 0; - count_pat_solutions = 0; - for (int i = 0; i < 32; i++) - { - min_position[i] = 0; - max_position[i] = 0; - now[i] = 0; - } + int max_row_count, max_column_count; + int carbon; + const char **col_name, **row_name; + int count_rows, count_optimize; + int col_phases, col_redox, col_epsilon, col_ph, col_water, + col_isotopes, col_phase_isotopes; + int row_mb, row_fract, row_charge, row_carbon, row_isotopes, + row_epsilon, row_isotope_epsilon, row_water; + LDBLE *inv_zero, *array1, *inv_res, *inv_delta1, *delta2, *delta3, *inv_cu, + *delta_save; + LDBLE *min_delta, *max_delta; + int *inv_iu, *inv_is; + int klmd, nklmd, n2d, kode, iter; + LDBLE toler, error, max_pct, scaled_error; + struct master *master_alk; + int *row_back, *col_back; + unsigned long *good, *bad, *minimal; + int max_good, max_bad, max_minimal; + int count_good, count_bad, count_minimal, count_calls; + unsigned long soln_bits, phase_bits, current_bits, temp_bits; + FILE *netpath_file; + int count_inverse_models, count_pat_solutions; + int min_position[32], max_position[32], now[32]; + std::vector inverse_heading_names; + /* kinetics.cpp ------------------------------- */ - count_pp = count_pg = count_ss = 0; - cvode_kinetics_ptr = NULL; - cvode_test = FALSE; - cvode_error = FALSE; - cvode_n_user = -99; - cvode_n_reactions = -99; - cvode_step_fraction = 0.0; - cvode_rate_sim_time = 0.0; - cvode_rate_sim_time_start = 0.0; - cvode_last_good_time = 0.0; - cvode_prev_good_time = 0.0; - cvode_last_good_y = NULL; - cvode_prev_good_y = NULL; - kinetics_machEnv = NULL; - kinetics_y = NULL; - kinetics_abstol = NULL; - kinetics_cvode_mem = NULL; - cvode_pp_assemblage_save= NULL; - cvode_ss_assemblage_save= NULL; - m_original = NULL; - m_temp = NULL; - rk_moles = NULL; - set_and_run_attempt = 0; - x0_moles = NULL; +public: + int count_pp, count_pg, count_ss; + void *cvode_kinetics_ptr; + int cvode_test; + int cvode_error; + int cvode_n_user; + int cvode_n_reactions; + realtype cvode_step_fraction; + realtype cvode_rate_sim_time; + realtype cvode_rate_sim_time_start; + realtype cvode_last_good_time; + realtype cvode_prev_good_time; + N_Vector cvode_last_good_y; + N_Vector cvode_prev_good_y; + M_Env kinetics_machEnv; + N_Vector kinetics_y, kinetics_abstol; + void *kinetics_cvode_mem; + cxxSSassemblage *cvode_ss_assemblage_save; + cxxPPassemblage *cvode_pp_assemblage_save; +protected: + LDBLE *m_original; + LDBLE *m_temp; + LDBLE *rk_moles; + int set_and_run_attempt; + LDBLE *x0_moles; + /* model.cpp ------------------------------- */ - gas_in = FALSE; - min_value = 1e-10; - normal = NULL; - ineq_array = NULL; - res = NULL; - cu = NULL; - zero = NULL; - delta1 = NULL; - iu = NULL; - is = NULL; - back_eq = NULL; - normal_max = 0; - ineq_array_max = 0; - res_max = 0; - cu_max = 0; - zero_max = 0; - delta1_max = 0; - iu_max = 0; - is_max = 0; - back_eq_max = 0; + int gas_in; + LDBLE min_value; + LDBLE *normal, *ineq_array, *res, *cu, *zero, *delta1; + int *iu, *is, *back_eq; + int normal_max, ineq_array_max, res_max, cu_max, zero_max, + delta1_max, iu_max, is_max, back_eq_max; + /* phrq_io_output.cpp ------------------------------- */ - forward_output_to_log = 0; + int forward_output_to_log; + /* phreeqc_files.cpp ------------------------------- */ - default_data_base = string_duplicate("phreeqc.dat"); + char *default_data_base; #ifdef PHREEQ98 int outputlinenr; char *LogFileNameC; @@ -2572,27 +2641,54 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) //AW = 0; //VP = 0; //DW0 = 0; - + full_pitzer = pSrc->full_pitzer; + always_full_pitzer = pSrc->always_full_pitzer; ICON = pSrc->ICON; + IC = pSrc->IC; /* pitz_params = NULL; count_pitz_param = 0; max_pitz_param = 100; */ + for (int i = 0; i < pSrc->count_pitz_param; i++) { pitz_param_store(pSrc->pitz_params[i], true); } - + pitz_param_map = pSrc->pitz_param_map; // auto pitz_param_map /* theta_params = 0; count_theta_param = 0; max_theta_param = 100; use_etheta = TRUE; + */ + count_theta_param = pSrc->count_theta_param; + max_theta_param = count_theta_param; + space((void **)((void *)&theta_params), count_theta_param, &max_theta_param, + sizeof(struct theta_param *)); + if (pSrc->theta_params != NULL) + { + //theta_params = (struct theta_param **) malloc((size_t)count_theta_param * sizeof(struct theta_param *)); + for (int i = 0; i < count_theta_param; i++) + { + theta_params[i] = theta_param_alloc(); + memcpy(theta_params[i], pSrc->theta_params[i], sizeof(struct theta_param)); + } + } + use_etheta = pSrc->use_etheta; + /* OTEMP = -100.0; OPRESS = -100.0; - A0 = 0; + A0 = 0; + aphi + */ + if (pSrc->aphi != NULL) + { + aphi = (struct pitz_param *) malloc(sizeof(struct pitz_param)); + memcpy(aphi, pSrc->aphi, sizeof(struct pitz_param)); + } + /* spec = NULL; cations = NULL; anions = NULL; @@ -2671,11 +2767,13 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) sit_M = NULL; sit_LGAMMA = NULL; */ - + count_sit_param = 0; //pSrc->count_sit_param; + max_sit_param = 1; // count_sit_param; for (int i = 0; i < pSrc->count_sit_param; i++) { sit_param_store(pSrc->sit_params[i], true); } + sit_param_map = pSrc->sit_param_map; /* tidy.cpp ------------------------------- */ //a0 = 0; //a1 = 0; @@ -2714,6 +2812,9 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) heat_mix_f_m = pSrc->heat_mix_f_m; warn_MCD_X = pSrc->warn_MCD_X; warn_fixed_Surf = pSrc->warn_fixed_Surf; + current_x = pSrc->current_x; + current_A = pSrc->current_A; + fix_current = pSrc->fix_current; #ifdef PHREEQ98 int AutoLoadOutputFile, CreateToC; @@ -2729,6 +2830,11 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc) //{ // keycount.push_back(0); //} + spinner = pSrc->spinner; + gfw_map = pSrc->gfw_map; + rates_map = pSrc->rates_map; + sum_species_map = pSrc->sum_species_map; + sum_species_map_db = pSrc->sum_species_map_db; // make sure new_model gets set this->keycount[Keywords::KEY_SOLUTION_SPECIES] = 1; @@ -2752,7 +2858,10 @@ Phreeqc &Phreeqc::operator=(const Phreeqc &rhs) } // copy Phreeqc object to this - this->phrq_io = rhs.phrq_io; + //this->phrq_io = rhs.phrq_io; + //this->phrq_io = new PHRQ_io; + this->phrq_io->Set_output_ostream(&std::cout); + this->phrq_io->Set_error_ostream(&std::cerr); this->init(); this->initialize(); this->InternalCopy(&rhs); diff --git a/phreeqcpp/Phreeqc.h b/phreeqcpp/Phreeqc.h index f9177850..eca08283 100644 --- a/phreeqcpp/Phreeqc.h +++ b/phreeqcpp/Phreeqc.h @@ -912,6 +912,7 @@ protected: public: struct rate *rate_bsearch(char *ptr, int *j); int rate_free(struct rate *rate_ptr); + struct rate * rate_copy(struct rate *rate_ptr); struct rate *rate_search(const char *name, int *n); int rate_sort(void); struct reaction *rxn_alloc(int ntokens); diff --git a/phreeqcpp/Surface.cxx b/phreeqcpp/Surface.cxx index b2e7ef2f..99a986de 100644 --- a/phreeqcpp/Surface.cxx +++ b/phreeqcpp/Surface.cxx @@ -24,6 +24,7 @@ cxxSurface::cxxSurface(PHRQ_io *io) : cxxNumKeyword(io) { new_def = false; + tidied = false; type = DDL; dl_type = NO_DL; sites_units = SITES_ABSOLUTE; @@ -42,6 +43,7 @@ cxxNumKeyword(io) { this->n_user = this->n_user_end = l_n_user; this->new_def = false; + this->tidied = true; type = DDL; dl_type = NO_DL; sites_units = SITES_ABSOLUTE; @@ -215,6 +217,8 @@ cxxSurface::dump_raw(std::ostream & s_oss, unsigned int indent, int *n_out) cons s_oss << indent1; s_oss << "-new_def " << this->new_def << "\n"; s_oss << indent1; + s_oss << "-tidied " << this->tidied << "\n"; + s_oss << indent1; s_oss << "-sites_units " << this->sites_units << "\n"; s_oss << indent1; s_oss << "-solution_equilibria " << this->solution_equilibria << "\n"; @@ -243,6 +247,7 @@ cxxSurface::read_raw(CParser & parser, bool check) // Read surface number and description this->read_number_description(parser); this->Set_new_def(false); + this->Set_tidied(true); bool only_counter_ions_defined(false); bool thickness_defined(false); @@ -515,6 +520,15 @@ cxxSurface::read_raw(CParser & parser, bool check) PHRQ_io::OT_CONTINUE); } break; + case 18: // tidied + if (!(parser.get_iss() >> this->tidied)) + { + this->tidied = false; + parser.incr_input_error(); + parser.error_msg("Expected boolean value for tidied.", + PHRQ_io::OT_CONTINUE); + } + break; } if (opt == CParser::OPT_EOF || opt == CParser::OPT_KEYWORD) break; @@ -768,6 +782,7 @@ cxxSurface::Serialize(Dictionary & dictionary, std::vector < int >&ints, } } ints.push_back(this->new_def ? 1 : 0); + ints.push_back(this->tidied ? 1 : 0); ints.push_back((int) this->type); ints.push_back((int) this->dl_type); ints.push_back((int) this->sites_units); @@ -813,6 +828,7 @@ cxxSurface::Deserialize(Dictionary & dictionary, std::vector < int >&ints, } } this->new_def = (ints[ii++] != 0); + this->tidied = (ints[ii++] != 0); this->type = (SURFACE_TYPE) ints[ii++]; this->dl_type = (DIFFUSE_LAYER_TYPE) ints[ii++]; this->sites_units = (SITES_UNITS) ints[ii++]; @@ -847,6 +863,7 @@ const std::vector< std::string >::value_type temp_vopts[] = { std::vector< std::string >::value_type("new_def"), // 14 std::vector< std::string >::value_type("solution_equilibria"), // 15 std::vector< std::string >::value_type("n_solution"), // 16 - std::vector< std::string >::value_type("totals") // 17 + std::vector< std::string >::value_type("totals"), // 17 + std::vector< std::string >::value_type("tidied") // 18 }; const std::vector< std::string > cxxSurface::vopts(temp_vopts, temp_vopts + sizeof temp_vopts / sizeof temp_vopts[0]); diff --git a/phreeqcpp/Surface.h b/phreeqcpp/Surface.h index 4a197453..7f97bce3 100644 --- a/phreeqcpp/Surface.h +++ b/phreeqcpp/Surface.h @@ -49,6 +49,8 @@ public: void Set_surface_charges(std::vector < cxxSurfaceCharge > &sc) {this->surface_charges = sc;} bool Get_new_def(void) {return new_def;} void Set_new_def(bool tf) {new_def = tf;} + bool Get_tidied(void) { return tidied; } + void Set_tidied(bool tf) { tidied = tf; } SURFACE_TYPE Get_type(void) const {return this->type;} void Set_type(SURFACE_TYPE t) {this->type = t;} DIFFUSE_LAYER_TYPE Get_dl_type(void) const {return dl_type;} @@ -80,6 +82,7 @@ protected: std::vector < cxxSurfaceComp > surface_comps; std::vector < cxxSurfaceCharge > surface_charges; bool new_def; + bool tidied; SURFACE_TYPE type; DIFFUSE_LAYER_TYPE dl_type; SITES_UNITS sites_units; diff --git a/phreeqcpp/basicsubs.cpp b/phreeqcpp/basicsubs.cpp index 525dc18a..003726aa 100644 --- a/phreeqcpp/basicsubs.cpp +++ b/phreeqcpp/basicsubs.cpp @@ -1392,8 +1392,9 @@ get_calculate_value(const char *name) { error_string = sformatf( "CALC_VALUE Basic function, %s not found.", name); - error_msg(error_string, CONTINUE); - input_error++; + //error_msg(error_string, CONTINUE); + //input_error++; + warning_msg(error_string); return (MISSING); } if (name == NULL) diff --git a/phreeqcpp/inverse.cpp b/phreeqcpp/inverse.cpp index 88f50cb1..16c9c1a2 100644 --- a/phreeqcpp/inverse.cpp +++ b/phreeqcpp/inverse.cpp @@ -1522,7 +1522,7 @@ solve_with_mask(struct inverse *inv_ptr, unsigned long cur_bits) } kode = 1; - iter = 1000; + iter = 100000; count_calls++; #ifdef INVERSE_CL1MP diff --git a/phreeqcpp/isotopes.cpp b/phreeqcpp/isotopes.cpp index 4ced82a5..9ed21205 100644 --- a/phreeqcpp/isotopes.cpp +++ b/phreeqcpp/isotopes.cpp @@ -1730,7 +1730,10 @@ calculate_value_init(struct calculate_value *calculate_value_ptr) if (calculate_value_ptr) { calculate_value_ptr->name = NULL; + calculate_value_ptr->value = 0.0; calculate_value_ptr->commands = NULL; + calculate_value_ptr->new_def = TRUE; + calculate_value_ptr->calculated = FALSE; calculate_value_ptr->linebase = NULL; calculate_value_ptr->varbase = NULL; calculate_value_ptr->loopbase = NULL; diff --git a/phreeqcpp/mainsubs.cpp b/phreeqcpp/mainsubs.cpp index 62aed0a3..f0512ccf 100644 --- a/phreeqcpp/mainsubs.cpp +++ b/phreeqcpp/mainsubs.cpp @@ -2445,7 +2445,38 @@ run_simulations(void) */ for (simulation = 1;; simulation++) { - +#ifdef TEST_COPY_OPERATOR + { + //int simulation_save = simulation; + Phreeqc phreeqc_new; + phreeqc_new = *this; + PHRQ_io *temp_io = this->phrq_io; + std::vector so_ostreams; + { + std::map::iterator so_it = this->SelectedOutput_map.begin(); + for (; so_it != this->SelectedOutput_map.end(); so_it++) + { + so_ostreams.push_back(so_it->second.Get_punch_ostream()); + so_it->second.Set_punch_ostream(NULL); + } + } + this->clean_up(); + this->init(); + this->initialize(); + this->phrq_io = temp_io; + this->InternalCopy(&phreeqc_new); + { + size_t i = 0; + std::map::iterator so_it = this->SelectedOutput_map.begin(); + for (; so_it != this->SelectedOutput_map.end(); so_it++) + { + so_it->second.Set_punch_ostream(so_ostreams[i++]); + } + } + //this->simulation = simulation_save; + //delete phreeqc_new.Get_phrq_io(); + } +#endif #if defined PHREEQ98 AddSeries = !connect_simulations; #endif diff --git a/phreeqcpp/prep.cpp b/phreeqcpp/prep.cpp index 8deef0a1..9ae513a1 100644 --- a/phreeqcpp/prep.cpp +++ b/phreeqcpp/prep.cpp @@ -5601,7 +5601,8 @@ calc_lk_phase(phase *p_ptr, LDBLE TK, LDBLE pa) } } } - else if (s_x[i]->millero[0]) + //else if (s_x[i]->millero[0]) + else if (s_ptr->millero[0]) { /* Millero volume at I = 0... */ d_v += s_ptr->millero[0] + tc * (s_ptr->millero[1] + tc * s_ptr->millero[2]); diff --git a/phreeqcpp/print.cpp b/phreeqcpp/print.cpp index ad04d4aa..2808c85d 100644 --- a/phreeqcpp/print.cpp +++ b/phreeqcpp/print.cpp @@ -12,6 +12,7 @@ #include "SSassemblage.h" #include "cxxKinetics.h" #include "Solution.h" +#include "Surface.h" /* ---------------------------------------------------------------------- */ int Phreeqc:: array_print(LDBLE * array_l, int row_count, int column_count, diff --git a/phreeqcpp/structures.cpp b/phreeqcpp/structures.cpp index 67c39f2a..0e8b7000 100644 --- a/phreeqcpp/structures.cpp +++ b/phreeqcpp/structures.cpp @@ -1671,6 +1671,26 @@ rate_free(struct rate *rate_ptr) return (OK); } +/* ---------------------------------------------------------------------- */ +struct rate * Phreeqc:: +rate_copy(struct rate *rate_ptr) +/* ---------------------------------------------------------------------- */ +{ + /* + * Copies a rate to new allocated space + */ + if (rate_ptr == NULL) + return (NULL); + struct rate * rate_new = (struct rate *) PHRQ_malloc(sizeof(struct rate)); + if (rate_new == NULL) malloc_error(); + rate_new->commands = string_duplicate(rate_ptr->commands); + rate_new->new_def = TRUE; + rate_new->linebase = NULL; + rate_new->varbase = NULL; + rate_new->loopbase = NULL; + return (rate_new); +} + /* ---------------------------------------------------------------------- */ struct rate * Phreeqc:: rate_search(const char *name_in, int *n) diff --git a/phreeqcpp/tidy.cpp b/phreeqcpp/tidy.cpp index 865f01c9..72e55366 100644 --- a/phreeqcpp/tidy.cpp +++ b/phreeqcpp/tidy.cpp @@ -3066,6 +3066,8 @@ tidy_surface(void) } //if (!kit->second.Get_new_def()) continue; surface_ptr = &(kit->second); + if (surface_ptr->Get_tidied()) continue; + surface_ptr->Set_tidied(true); // ccm incompatible with Donnan or diffuse_layer if (surface_ptr->Get_type() == cxxSurface::CCM) {