From c7111f77b203eb081111dc8c47aca5b31b8a9292 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Tue, 28 May 2019 17:05:02 -0600 Subject: [PATCH 1/4] Sort of works, still bugs and serious errors compared to explicit --- Phreeqc.cpp | 2 + Phreeqc.h | 7 +- global_structures.h | 4 +- kinetics.cpp | 57 ++- phqalloc.cpp | 1 + readtr.cpp | 30 +- transport.cpp | 1132 ++++++++++++++++++++++++++++++++++++------- 7 files changed, 1024 insertions(+), 209 deletions(-) diff --git a/Phreeqc.cpp b/Phreeqc.cpp index f7e499c3..00e13a2b 100644 --- a/Phreeqc.cpp +++ b/Phreeqc.cpp @@ -711,6 +711,8 @@ void Phreeqc::init(void) all_cells = 0; multi_Dflag = FALSE; interlayer_Dflag = FALSE; + implicit = FALSE; + max_mixf = 1.0; default_Dw = 0; correct_Dw = 0; multi_Dpor = 0; diff --git a/Phreeqc.h b/Phreeqc.h index 0ca33f0b..f1f03b64 100644 --- a/Phreeqc.h +++ b/Phreeqc.h @@ -1065,10 +1065,11 @@ public: LDBLE viscosity(void); LDBLE calc_vm_Cl(void); int multi_D(LDBLE DDt, int mobile_cell, int stagnant); - int find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant); + LDBLE find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant); + void diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant); int fill_spec(int cell_no); void define_ct_structures(void); - int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec); + int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec, int i); static int sort_species_name(const void *ptr1, const void *ptr2); int disp_surf(LDBLE stagkin_time); int diff_stag_surf(int mobile_cell); @@ -1426,6 +1427,8 @@ protected: int old_cells, max_cells, all_cells; int multi_Dflag; /* signals calc'n of multicomponent diffusion */ int interlayer_Dflag; /* multicomponent diffusion and diffusion through interlayer porosity */ + int implicit; /* implicit calculation of diffusion */ + LDBLE max_mixf; /* the maximum value of the implicit mixfactor = De * Dt / (Dx^2) */ LDBLE default_Dw; /* default species diffusion coefficient in water at 25oC, m2/s */ int correct_Dw; /* if true, Dw is adapted in calc_SC */ LDBLE multi_Dpor; /* uniform porosity of free porewater in solid medium */ diff --git a/global_structures.h b/global_structures.h index 90cb786d..7b759984 100644 --- a/global_structures.h +++ b/global_structures.h @@ -1052,12 +1052,12 @@ struct sol_D struct J_ij { const char *name; - LDBLE tot1, tot2; /* species change in cells i and j */ + LDBLE tot1, tot2, charge; /* species change in cells i and j */ }; struct M_S { const char *name; - LDBLE tot1, tot2; + LDBLE tot1, tot2, charge; /* master species transport in cells i and j */ }; // Pitzer definitions typedef enum diff --git a/kinetics.cpp b/kinetics.cpp index 9b938aee..c13af7d6 100644 --- a/kinetics.cpp +++ b/kinetics.cpp @@ -2478,7 +2478,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction) * Rates and moles of each reaction are calculated in calc_kinetic_reaction * Total number of moles in reaction is stored in kinetics[i].totals */ - + //int increase_tol = 0; // appt int converge, m_iter; int pr_all_save; int nsaver; @@ -2529,9 +2529,10 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction) converge = set_and_run_wrapper(i, use_mix, FALSE, nsaver, step_fraction); if (converge == MASS_BALANCE) - error_msg - ("Negative concentration in system. Stopping calculation.", - STOP); + { + error_string = sformatf("Negative concentration in solution %d. Stopping calculation.", cell_no); + error_msg(error_string, STOP); + } run_reactions_iterations += iterations; } else @@ -2603,9 +2604,10 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction) else converge = set_and_run_wrapper(i, use_mix, FALSE, i, step_fraction); if (converge == MASS_BALANCE) - error_msg - ("Negative concentration in system. Stopping calculation.", - STOP); + { + error_string = sformatf("Negative concentration in solution %d. Stopping calculation.", cell_no); + error_msg(error_string, STOP); + } saver(); pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, i); ss_assemblage_ptr = Utilities::Rxn_find(Rxn_ss_assemblage_map, i); @@ -2675,6 +2677,7 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction) cvode_mem = CVodeMalloc(n_reactions, f, 0.0, y, ADAMS, FUNCTIONAL, SV, &reltol, abstol, NULL, NULL, FALSE, iopt, ropt, machEnv); iopt[MXSTEP] is maximum number of steps that CVODE tries. */ + //iopt[SLDET] = TRUE; // appt iopt[MXSTEP] = kinetics_ptr->Get_cvode_steps(); iopt[MAXORD] = kinetics_ptr->Get_cvode_order(); kinetics_cvode_mem = @@ -2704,30 +2707,49 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction) */ m_iter = 0; sum_t = 0; - RESTART: + RESTART: while (flag != SUCCESS) { sum_t += cvode_last_good_time; - if (state != TRANSPORT) { - error_string = sformatf( - "CVode incomplete at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d, continuing...\n", - (int)iopt[NST], cell_no, (double)sum_t, m_iter + 1); - warning_msg(error_string); + error_string = sformatf("CV_ODE: Time: %8.2e s. Delta t: %8.2e s. Calls: %d.", (double)(sum_t), (double) cvode_last_good_time, m_iter); + status(0, error_string, true); } + //if (state != TRANSPORT) + //{ + // error_string = sformatf( + // "CVode incomplete at cvode_steps %d. Cell: %d. Time: %8.2e s. Cvode calls: %d, continuing...\n", + // (int)iopt[NST], cell_no, (double)sum_t, m_iter + 1); + // warning_msg(error_string); + //} #ifdef DEBUG_KINETICS if (m_iter > 5) dump_kinetics_stderr(cell_no); #endif + //if (m_iter > 0.5 * kinetics_ptr->Get_bad_step_max() && + // (cvode_last_good_time < 1e-6 || cvode_last_good_time < 1e-6 * tout)) // appt + //{ + // if (increase_tol < 3) + // { + // increase_tol += 1; + // for (size_t j = 0; j < kinetics_ptr->Get_kinetics_comps().size(); j++) + // { + // cxxKineticsComp * kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]); + // LDBLE tr = kinetics_comp_ptr->Get_tol() * 10.0; + // kinetics_comp_ptr->Set_tol(tr); + // tr += 0; + // } + // } + //} cvode_last_good_time = 0; if (++m_iter >= kinetics_ptr->Get_bad_step_max()) { m_temp = (LDBLE *) free_check_null(m_temp); m_original = (LDBLE *) free_check_null(m_original); error_string = sformatf( - "Restart of integration at cvode_steps %d. Cell: %d\tTime: %e\tCvode calls: %d.", - (int)iopt[NST], cell_no, (double)sum_t, m_iter - 1); + "CVode is at maximum calls: %d. Cell: %d. Time: %8.2e s\nERROR: Please increase the maximum calls with -bad_step_max.", + m_iter, cell_no, (double)sum_t); error_msg(error_string, STOP); } tout1 = tout - sum_t; @@ -2810,6 +2832,11 @@ run_reactions(int i, LDBLE kin_time, int use_mix, LDBLE step_fraction) free_cvode(); use.Set_mix_in(use_save.Get_mix_in()); use.Set_mix_ptr(use_save.Get_mix_ptr()); + + error_string = sformatf("CV_ODE: Final Delta t: %8.2e s. Calls: %d. ", (double)cvode_last_good_time, m_iter); + status(0, error_string, true); + + //status(0, NULL); } rate_sim_time = rate_sim_time_start + kin_time; diff --git a/phqalloc.cpp b/phqalloc.cpp index 3823f397..8b7ce14b 100644 --- a/phqalloc.cpp +++ b/phqalloc.cpp @@ -142,6 +142,7 @@ PHRQ_calloc(size_t num, size_t size assert((s_pTail == NULL) || (s_pTail->pNext == NULL)); p = (PHRQMemHeader *) malloc(sizeof(PHRQMemHeader) + size * num); + //p = (PHRQMemHeader *) calloc(1, sizeof(PHRQMemHeader) + size * num); // appt if (p == NULL) return NULL; diff --git a/readtr.cpp b/readtr.cpp index 566fb781..fee12eea 100644 --- a/readtr.cpp +++ b/readtr.cpp @@ -94,9 +94,10 @@ read_transport(void) "porosities", /* 42 */ "porosity", /* 43 */ "fix_current", /* 44 */ - "current" /* 45 */ + "current", /* 45 */ + "implicit" /* 46 */ }; - int count_opt_list = 46; + int count_opt_list = 47; strcpy(file_name, "phreeqc.dmp"); /* @@ -682,6 +683,31 @@ read_transport(void) } opt_save = OPTION_DEFAULT; break; + case 46: /* implicit diffusion */ + copy_token(token, &next_char, &l); + str_tolower(token); + if (strstr(token, "f") == token) + implicit = FALSE; + else if (strstr(token, "t") == token) + implicit = TRUE; + else + { + input_error++; + error_msg + ("Expected flag for implicit diffusion calc`s: 'true' or 'false'.", + CONTINUE); + } + if (copy_token(token, &next_char, &l) == DIGIT) + { + sscanf(token, SCANFORMAT, &max_mixf); + } + else + { + //warning_msg("Expected the maximal value for the mixfactor (= D * Dt / Dx^2) in implicit calc`s of diffusion."); + max_mixf = 1.0; + } + opt_save = OPTION_DEFAULT; + break; } if (return_value == EOF || return_value == KEYWORD) break; diff --git a/transport.cpp b/transport.cpp index 9f6d9c6b..2230b10a 100644 --- a/transport.cpp +++ b/transport.cpp @@ -10,9 +10,18 @@ #include LDBLE F_Re3 = F_C_MOL / (R_KJ_DEG_MOL * 1e3); -LDBLE tk_x2; // average tx_x of icell and jcell +LDBLE tk_x2; // average tk_x of icell and jcell LDBLE dV_dcell; // difference in Volt among icell and jcell int find_current; +char token[MAX_LENGTH]; + +// implicit... +std::set dif_spec_names; +std::set dif_els_names; +std::map> neg_moles; +std::map els; +double *y, *Ct1, *Ct2, *l_tk_x2, **A, **mixf, mixf_comp_size = 0; + struct CURRENT_CELLS { LDBLE dif, ele, R; // diffusive and electric components, relative cell resistance @@ -25,11 +34,14 @@ struct V_M // For calculating Vinograd and McBain's zero-charge, diffusive tra }; struct CT /* summed parts of V_M and mcd transfer in a timestep for all cells, for free + DL water */ { - LDBLE dl_s, Dz2c, Dz2c_dl, visc1, visc2, J_ij_sum; + LDBLE kgw, dl_s, Dz2c, visc1, visc2, J_ij_sum; LDBLE A_ij_il, Dz2c_il, mixf_il; int J_ij_count_spec, J_ij_il_count_spec; struct V_M *v_m, *v_m_il; struct J_ij *J_ij, *J_ij_il; + int count_m_s; + struct M_S *m_s; + int v_m_size, J_ij_size, m_s_size; } *ct = NULL; struct MOLES_ADDED /* total moles added to balance negative conc's */ { @@ -73,14 +85,11 @@ transport(void) { input_error++; error_string = sformatf( - "Solution %d is needed for transport, but is not defined.", - i); + "Solution %d is needed for transport, but is not defined.", i); error_msg(error_string, CONTINUE); } else - { cell_data[i].temp = use.Get_solution_ptr()->Get_tc(); - } } if (multi_Dflag) @@ -88,36 +97,6 @@ transport(void) sol_D = (struct sol_D *) PHRQ_malloc((size_t) (all_cells)* sizeof(struct sol_D)); if (sol_D == NULL) malloc_error(); - //sol_D_dbg = sol_D; - - ct = (struct CT *) PHRQ_malloc((size_t) (all_cells)* sizeof(struct CT)); - if (ct == NULL) - malloc_error(); - { - for (int i = 0; i < all_cells; i++) - { - ct[i].dl_s = 0.0; - ct[i].Dz2c = 0.0; - ct[i].Dz2c_dl = 0.0; - ct[i].visc1 = 0.0; - ct[i].visc2 = 0.0; - ct[i].J_ij_sum = 0.0; - ct[i].A_ij_il = 0.0; - ct[i].Dz2c_il = 0.0; - ct[i].mixf_il = 0.0; - ct[i].J_ij_count_spec = 0; - ct[i].J_ij_il_count_spec = 0; - ct[i].v_m = NULL; - ct[i].v_m_il = NULL; - ct[i].J_ij = NULL; - ct[i].J_ij_il = NULL; - } - } - - moles_added = (struct MOLES_ADDED *) PHRQ_malloc((size_t) (count_elements)* sizeof(struct MOLES_ADDED)); - if (moles_added == NULL) - malloc_error(); - for (i = 0; i < all_cells; i++) { sol_D[i].count_spec = 0; @@ -125,14 +104,45 @@ transport(void) sol_D[i].exch_total = 0; sol_D[i].x_max = 0; sol_D[i].spec = NULL; + sol_D[i].spec_size = 0; + sol_D[i].tk_x = 298.15; } + //sol_D_dbg = sol_D; + + ct = (struct CT *) PHRQ_malloc((size_t) (all_cells)* sizeof(struct CT)); + if (ct == NULL) + malloc_error(); + for (i = 0; i < all_cells; i++) + { + ct[i].kgw = 1.0; + ct[i].dl_s = 0.0; + ct[i].Dz2c = 0.0; + ct[i].visc1 = 0.0; + ct[i].visc2 = 0.0; + ct[i].J_ij_sum = 0.0; + ct[i].A_ij_il = 0.0; + ct[i].Dz2c_il = 0.0; + ct[i].mixf_il = 0.0; + ct[i].J_ij_count_spec = 0; + ct[i].J_ij_il_count_spec = 0; + ct[i].v_m = NULL; + ct[i].v_m_il = NULL; + ct[i].J_ij = NULL; + ct[i].J_ij_il = NULL; + ct[i].m_s = NULL; + ct[i].v_m_size = ct[i].J_ij_size = ct[i].m_s_size = 0; + } + + moles_added = (struct MOLES_ADDED *) PHRQ_malloc((size_t) (count_elements)* sizeof(struct MOLES_ADDED)); + if (moles_added == NULL) + malloc_error(); for (i = 0; i < count_elements; i++) { moles_added[i].name = NULL; moles_added[i].moles = 0; } } - try + try { /* check solution 0 */ use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, 0)); @@ -188,7 +198,7 @@ transport(void) } } - if (fix_current && !dV_dcell) + if (fix_current && !dV_dcell && !implicit) { warning_msg("fix_current (A) was defined, but potential in a boundary cell was not.\n\tUsing 1e-8 V in cell 0 for calculating the potentials."); cell_data[0].potV = 1e-8; @@ -202,7 +212,6 @@ transport(void) error_string = sformatf( "Electro-diffusion cannot be combined with advective transport."); error_msg(error_string, CONTINUE); - free_check_null(sol_D); } if (!multi_Dflag) { @@ -210,7 +219,6 @@ transport(void) error_string = sformatf( "Electrical Field (potential) was defined, but needs -multi_D."); error_msg(error_string, CONTINUE); - free_check_null(sol_D); } else { @@ -228,7 +236,7 @@ transport(void) "Electrical Field (potential) was defined, assuming constant boundary condition for last cell."); warning_msg(error_string); } - current_cells = (struct CURRENT_CELLS *) PHRQ_malloc((size_t) + current_cells = (struct CURRENT_CELLS *) PHRQ_malloc((size_t) (count_cells + 1) * sizeof(struct CURRENT_CELLS)); if (current_cells == NULL) malloc_error(); @@ -240,6 +248,34 @@ transport(void) } } } + if (implicit && interlayer_Dflag) + { + input_error++; + error_string = sformatf( + "Interlayer diffusion cannot be calculated implicitly. Please remove -implicit true"); + error_msg(error_string, CONTINUE); + } + if (implicit && !multi_Dflag) + { + input_error++; + error_string = sformatf( + "Implicit diffusion needs diffusion coefficients for individual species. Please add -multi_d true"); + error_msg(error_string, CONTINUE); + } + if (implicit && current_cells == NULL) + { + current_cells = (struct CURRENT_CELLS *) PHRQ_malloc((size_t) + (count_cells + 1) * sizeof(struct CURRENT_CELLS)); + if (current_cells == NULL) + malloc_error(); + for (int i = 0; i < count_cells + 1; i++) + { + current_cells[i].dif = 0.0; + current_cells[i].ele = 0.0; + current_cells[i].R = 0.0; + } + } + /* * First equilibrate solutions */ @@ -247,9 +283,9 @@ transport(void) transport_step = 0; for (i = 0; i <= count_cells + 1; i++) { - if ((bcon_first == 2 && i == 0) || - (bcon_last == 2 && i == count_cells + 1)) - continue; + //if ((bcon_first == 2 && i == 0) || + // (bcon_last == 2 && i == count_cells + 1)) + // continue; set_initial_moles(i); cell_no = i; set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0); @@ -378,12 +414,12 @@ transport(void) if (Utilities::Rxn_find(Rxn_solution_map, j) == NULL) error_msg ("Could not find mobile cell solution in TRANSPORT.", - STOP); + STOP); if (Utilities::Rxn_find(Rxn_solution_map, j_imm) == NULL) //error_msg //("Could not find immobile cell solution in TRANSPORT.", //STOP); - continue; + continue; water_m = Utilities::Rxn_find(Rxn_solution_map, j)->Get_mass_water(); water_imm = Utilities::Rxn_find(Rxn_solution_map, j_imm)->Get_mass_water(); /* @@ -432,8 +468,12 @@ transport(void) /* * Now transport */ - sprintf(token, "\nCalculating transport: %d (mobile) cells, %d shifts, %d mixruns...\n\n", - count_cells, count_shifts - transport_start + 1, nmix); + if (implicit) + sprintf(token, "\nCalculating implicit transport: %d (mobile) cells, %d shifts, %d mixruns, max. mixf = %g.\n\n", + count_cells, count_shifts - transport_start + 1, nmix, max_mixf); + else + sprintf(token, "\nCalculating transport: %d (mobile) cells, %d shifts, %d mixruns...\n\n", + count_cells, count_shifts - transport_start + 1, nmix); screen_msg(token); max_iter = 0; for (transport_step = transport_start; transport_step <= count_shifts; @@ -488,7 +528,7 @@ transport(void) dup_print(token, FALSE); } - if (heat_nmix > 0) + if (heat_nmix > 0 && !implicit) { heat_mix(heat_nmix); /* equilibrate again ... */ @@ -506,14 +546,16 @@ transport(void) { if (disp_surf(stagkin_time) == ERROR) error_msg("Error in surface transport, stopping.", - STOP); + STOP); } - if (multi_Dflag) + if (implicit) + diffuse_implicit(1, stagkin_time, FALSE); + else if (multi_Dflag) multi_D(stagkin_time, 1, FALSE); for (i = 0; i <= count_cells + 1; i++) { - if (!dV_dcell && (i == 0 || i == count_cells + 1)) + if (!dV_dcell && (i == 0 || i == count_cells + 1 && !implicit)) continue; if (overall_iterations > max_iter) max_iter = overall_iterations; @@ -521,16 +563,19 @@ transport(void) mixrun = j; if (multi_Dflag) sprintf(token, - "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", - transport_step, j, i, max_iter); + "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", + transport_step, j, i, max_iter); else sprintf(token, - "Transport step %3d. Mixrun %3d. Cell %3d. (Max. iter %3d)", - transport_step, j, i, max_iter); + "Transport step %3d. Mixrun %3d. Cell %3d. (Max. iter %3d)", + transport_step, j, i, max_iter); status(0, token); - if (i == 0 || i == count_cells + 1) + if (implicit || i == 0 || i == count_cells + 1) + { run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i + //print_punch(i, true); // appt debug + } else run_reactions(i, kin_time, DISP, step_fraction); // nsaver = -2 if (multi_Dflag) @@ -539,7 +584,7 @@ transport(void) /* punch and output file */ if (ishift == 0 && j == nmix && stag_data->count_stag == 0) print_punch(i, true); - if (i > 1) + if (!implicit && i > 1) Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1); saver(); @@ -565,7 +610,7 @@ transport(void) } } - if (!dV_dcell) + if (!dV_dcell && !implicit) Utilities::Rxn_copy(Rxn_solution_map, -2, count_cells); /* Stagnant zone mixing after completion of each diffusive/dispersive step ... */ @@ -646,7 +691,7 @@ transport(void) { if ((Utilities::Rxn_find(Rxn_surface_map, i) != NULL) && ((i == 0 && bcon_first == 3) - || (i == count_cells + 1 && bcon_last == 3))) + || (i == count_cells + 1 && bcon_last == 3))) { Rxn_surface_map.erase(i); } @@ -677,7 +722,7 @@ transport(void) /* * thermal diffusion when nmix = 0... */ - if ((nmix == 0) && (heat_nmix > 0)) + if ((nmix == 0) && (heat_nmix > 0)) // must add disp in implicit appt ?? { heat_mix(heat_nmix); /* equilibrate again ... */ @@ -699,12 +744,12 @@ transport(void) mixrun = 0; if (multi_Dflag) sprintf(token, - "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", - transport_step, 0, i, max_iter); + "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", + transport_step, 0, i, max_iter); else sprintf(token, - "Transport step %3d. Mixrun %3d. Cell %3d. (Max. iter %3d)", - transport_step, 0, i, max_iter); + "Transport step %3d. Mixrun %3d. Cell %3d. (Max. iter %3d)", + transport_step, 0, i, max_iter); status(0, token); run_reactions(i, kin_time, NOMIX, step_fraction); if (multi_Dflag == TRUE) @@ -785,7 +830,7 @@ transport(void) rate_sim_time_start += kin_time; rate_sim_time = rate_sim_time_start + kin_time; - if (heat_nmix > 0) + if (heat_nmix > 0 && !implicit) { heat_mix(heat_nmix); /* equilibrate again ... */ @@ -804,29 +849,37 @@ transport(void) error_msg("Error in surface transport, stopping.", STOP); } - if (multi_Dflag == TRUE) + if (implicit) + diffuse_implicit(1, stagkin_time, FALSE); + else if (multi_Dflag) multi_D(stagkin_time, 1, FALSE); /* for each cell in column */ for (i = 0; i <= count_cells + 1; i++) { if (!dV_dcell && (i == 0 || i == count_cells + 1)) + { + if (j == nmix && stag_data->count_stag == 0 && + (cell_data[0].print || cell_data[0].punch || + cell_data[count_cells + 1].print || cell_data[count_cells + 1].punch)) + print_punch(i, false); continue; + } if (overall_iterations > max_iter) max_iter = overall_iterations; cell_no = i; mixrun = j; if (multi_Dflag) sprintf(token, - "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", - transport_step, j, i, max_iter); + "Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)", + transport_step, j, i, max_iter); else sprintf(token, - "Transport step %3d. Mixrun %3d. Cell %3d. (Max. iter %3d)", - transport_step, j, i, max_iter); + "Transport step %3d. Mixrun %3d. Cell %3d. (Max. iter %3d)", + transport_step, j, i, max_iter); status(0, token); - if (i == 0 || i == count_cells + 1) + if (implicit || i == 0 || i == count_cells + 1) run_reactions(i, kin_time, NOMIX, step_fraction); else run_reactions(i, kin_time, DISP, step_fraction); @@ -834,7 +887,7 @@ transport(void) fill_spec(i); if (j == nmix && stag_data->count_stag == 0) print_punch(i, true); - if (i > 1) + if (!implicit && i > 1) Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1); saver(); @@ -896,6 +949,35 @@ transport(void) } screen_msg("\n"); + if (implicit) + { + std::map> ::iterator it1; + std::map ::iterator it2; + for (i = 0; i <= all_cells; i++) + { + if ((it1 = neg_moles.find(i)) != neg_moles.end() && (els = it1->second).size()) + { + for (it2 = els.begin(); it2 != els.end(); it2++) + { + for (int i1 = 0; i1 < count_elements; i1++) + { + if (moles_added[i1].name && !strcmp(moles_added[i1].name, it2->first.c_str())) + { + moles_added[i1].moles -= it2->second; + break; + } + else if (!moles_added[i1].moles) + { + moles_added[i1].name = string_duplicate(it2->first.c_str()); + moles_added[i1].moles -= it2->second; + break; + } + } + } + } + } + } + if (multi_Dflag && moles_added[0].moles > 0) { sprintf(token, @@ -959,6 +1041,7 @@ transport_cleanup(void) ct[i].v_m_il = (struct V_M *) free_check_null(ct[i].v_m_il); ct[i].J_ij = (struct J_ij *) free_check_null(ct[i].J_ij); ct[i].J_ij_il = (struct J_ij *) free_check_null(ct[i].J_ij_il); + ct[i].m_s = (struct M_S *) free_check_null(ct[i].m_s); } ct = (struct CT *) free_check_null(ct); for (int i = 0; i < count_elements; i++) @@ -967,6 +1050,28 @@ transport_cleanup(void) } moles_added = (struct MOLES_ADDED *) free_check_null(moles_added); } + if (implicit) + { + y = (LDBLE *)free_check_null(y); + Ct1 = (LDBLE *)free_check_null(Ct1); + Ct2 = (LDBLE *)free_check_null(Ct2); + l_tk_x2 = (LDBLE *)free_check_null(l_tk_x2); + for (i = 0; i < count_cells + 2; i++) + { + A[i] = (LDBLE *)free_check_null(A[i]); + mixf[i] = (LDBLE *)free_check_null(mixf[i]); + if (!dV_dcell) + { + cell_data[i].potV = 0.0; + use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, i)); + use.Get_solution_ptr()->Set_potV(0); + } + } + A = (LDBLE **)free_check_null(A); + mixf = (LDBLE **)free_check_null(mixf); + dif_spec_names.clear(); + mixf_comp_size = 0; + } current_cells = (struct CURRENT_CELLS *) free_check_null(current_cells); } /* ---------------------------------------------------------------------- */ @@ -1106,6 +1211,19 @@ init_mix(void) if (mcd_substeps > 1 && stag_data->count_stag > 0) l_nmix = (int) ceil(mcd_substeps); } + else if (implicit) + { + l_nmix = 1; + if (maxmix > max_mixf) + l_nmix = 1 + (int)floor(maxmix / max_mixf); + if (ishift != 0 && (bcon_first == 1 || bcon_last == 1)) + { // must be adjusted appt + if (l_nmix < 2) + l_nmix = 2; + } + if (mcd_substeps > 1) + l_nmix = (int)ceil(l_nmix * mcd_substeps); + } else { if (2.25 * maxmix + 1.0 > (double)INT_MAX) @@ -1113,7 +1231,7 @@ init_mix(void) m = (LDBLE *)free_check_null(m); m1 = (LDBLE *)free_check_null(m1); char token[MAX_LENGTH]; - sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please decrease time_step, or increase cell-lengths.", 2.25 * maxmix); + sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please set implicit true, or decrease time_step, or increase cell-lengths.", 2.25 * maxmix); error_msg(token, STOP); } if (bcon_first == 1 || bcon_last == 1) @@ -1234,7 +1352,7 @@ init_mix(void) m = (LDBLE *)free_check_null(m); m1 = (LDBLE *)free_check_null(m1); char token[MAX_LENGTH]; - sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please decrease time_step, or increase cell-lengths.", 1.5 * maxmix); + sprintf(token, "Calculated number of mixes %g, is beyond program limit,\nERROR: please set implicit true, or decrease time_step, or increase cell-lengths.", 1.5 * maxmix); error_msg(token, STOP); } l_nmix = 1 + (int) floor(1.5 * maxmix); @@ -1342,6 +1460,10 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) error_msg("Error in surface transport, stopping.", STOP); } + //if (implicit) // work this out appt + // diffuse_implicit(1, stagkin_time, FALSE); + //else if (multi_Dflag) + // multi_D(stagkin_time, 1, FALSE); if (multi_Dflag == TRUE) multi_D(1.0, i, TRUE); set_and_run_wrapper(i, STAG, FALSE, -2, 0.0); @@ -1419,19 +1541,23 @@ int Phreeqc:: init_heat_mix(int l_nmix) /* ---------------------------------------------------------------------- */ { - LDBLE lav, mixf, maxmix, corr_disp; + LDBLE lav, mixf, maxmix, corr_disp, l_diffc; int i, k, n; int l_heat_nmix; LDBLE t0; /* * Check for need to model thermal diffusion... */ - if (heat_diffc <= diffc) + if (heat_diffc <= diffc && !implicit) return (0); if (count_cells < 2) return (0); l_heat_nmix = 0; + if (implicit) + l_diffc = heat_diffc; + else + l_diffc = heat_diffc - diffc_tr; t0 = Utilities::Rxn_find(Rxn_solution_map, 0)->Get_tc(); for (i = 1; i <= count_cells; i++) { @@ -1488,15 +1614,13 @@ init_heat_mix(int l_nmix) if (bcon_last == 3) corr_disp += 1. / count_cells; } - if (l_nmix > 0) + if (l_nmix > 0 && !implicit) corr_disp /= l_nmix; maxmix = 0.0; for (i = 1; i < count_cells; i++) { lav = (cell_data[i + 1].length + cell_data[i].length) / 2; - mixf = - (heat_diffc - - diffc_tr) * timest * corr_disp / tempr / (lav * lav); + mixf = (l_diffc) * timest * corr_disp / tempr / (lav * lav); if (mixf > maxmix) maxmix = mixf; heat_mix_array[i + 1] = mixf; /* m[i] has mixf with lower cell */ @@ -1507,9 +1631,7 @@ init_heat_mix(int l_nmix) if (bcon_first == 1) { lav = cell_data[1].length; - mixf = - (heat_diffc - - diffc_tr) * timest * corr_disp / tempr / (lav * lav); + mixf = (l_diffc) * timest * corr_disp / tempr / (lav * lav); if (2 * mixf > maxmix) maxmix = 2 * mixf; heat_mix_array[1] = 2 * mixf; @@ -1520,9 +1642,7 @@ init_heat_mix(int l_nmix) if (bcon_last == 1) { lav = cell_data[count_cells].length; - mixf = - (heat_diffc - - diffc_tr) * timest * corr_disp / tempr / (lav * lav); + mixf = (l_diffc) * timest * corr_disp / tempr / (lav * lav); if (2 * mixf > maxmix) maxmix = 2 * mixf; heat_mix_array[count_cells + 1] = 2 * mixf; @@ -1536,9 +1656,26 @@ init_heat_mix(int l_nmix) l_heat_nmix = 0; else { - l_heat_nmix = 1 + (int) floor(3.0 * maxmix); - for (i = 1; i <= count_cells + 1; i++) - heat_mix_array[i] /= l_heat_nmix; + if (implicit) + { + // must correct Dw for temp... + fill_spec(i); + LDBLE viscos_f; + l_heat_nmix = l_nmix; + for (i = 1; i <= count_cells + 1; i++) + { + heat_mix_array[i - 1] = heat_mix_array[i] / l_heat_nmix; /* for implicit, m[i] has mixf with higher cell */ + viscos_f = sol_D[i - 1].viscos_f * exp(heat_diffc / sol_D[i - 1].tk_x - heat_diffc / 298.15); + viscos_f += sol_D[i].viscos_f * exp(heat_diffc / sol_D[i].tk_x - heat_diffc / 298.15); + heat_mix_array[i - 1] *= (viscos_f / 2); + } + } + else + { + l_heat_nmix = 1 + (int)floor(3.0 * maxmix); + for (i = 1; i <= count_cells + 1; i++) + heat_mix_array[i] /= l_heat_nmix; + } } return (l_heat_nmix); @@ -1694,7 +1831,7 @@ fill_spec(int l_cell_no) { /* copy species activities into sol_D.spec... */ - int i, i2, count_spec, count_exch_spec; + int i, i1, i2, i3, count_spec, count_exch_spec, size_xt; char token[MAX_LENGTH]; const char * name; struct species *s_ptr, *s_ptr2; @@ -1703,34 +1840,46 @@ fill_spec(int l_cell_no) LDBLE lm; LDBLE por, por_il, viscos_f, viscos_il_f, viscos; bool x_max_done = false; + std::set loc_spec_names; s_ptr2 = NULL; - sol_D[l_cell_no].spec = - (struct spec *) free_check_null(sol_D[l_cell_no].spec); - sol_D[l_cell_no].spec = - (struct spec *) PHRQ_malloc((size_t) count_species_list * - sizeof(struct spec)); + // for implicit + std::set ::iterator it; + std::pair ::iterator, bool> name_ret; + size_xt = 5; + + //sol_D[l_cell_no].spec = (struct spec *) free_check_null(sol_D[l_cell_no].spec); + if (sol_D[l_cell_no].spec == NULL) + { + sol_D[l_cell_no].spec = (struct spec *) PHRQ_malloc((size_t)(count_species_list + size_xt) * sizeof(struct spec)); + sol_D[l_cell_no].spec_size = count_species_list + size_xt; + } + else if (count_species_list + size_xt > sol_D[l_cell_no].spec_size) + { + sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(count_species_list + size_xt) * sizeof(struct spec)); + sol_D[l_cell_no].spec_size = count_species_list + size_xt; + } if (sol_D[l_cell_no].spec == NULL) malloc_error(); - sol_D[l_cell_no].spec_size = count_species_list; + + for (i = 0; i < sol_D[l_cell_no].spec_size; i++) { - for (int i = 0; i < count_species_list; i++) - { - sol_D[l_cell_no].spec[i].name = NULL; - sol_D[l_cell_no].spec[i].aq_name = NULL; - sol_D[l_cell_no].spec[i].type = -1; - sol_D[l_cell_no].spec[i].a = 0.0; - sol_D[l_cell_no].spec[i].lm = 0.0; - sol_D[l_cell_no].spec[i].lg = 0.0; - sol_D[l_cell_no].spec[i].c = 0.0; - sol_D[l_cell_no].spec[i].z = 0.0; - sol_D[l_cell_no].spec[i].Dwt = 0.0; - sol_D[l_cell_no].spec[i].dw_t = 0.0; - sol_D[l_cell_no].spec[i].erm_ddl = 0.0; - } + sol_D[l_cell_no].spec[i].name = NULL; + sol_D[l_cell_no].spec[i].aq_name = NULL; + sol_D[l_cell_no].spec[i].type = -1; + sol_D[l_cell_no].spec[i].a = 0.0; + sol_D[l_cell_no].spec[i].lm = 0.0; + sol_D[l_cell_no].spec[i].lg = 0.0; + sol_D[l_cell_no].spec[i].c = 0.0; + sol_D[l_cell_no].spec[i].z = 0.0; + sol_D[l_cell_no].spec[i].Dwt = 0.0; + sol_D[l_cell_no].spec[i].dw_t = 0.0; + sol_D[l_cell_no].spec[i].erm_ddl = 0.0; + sol_D[l_cell_no].count_exch_spec = sol_D[l_cell_no].count_spec = 0; } + sol_D[l_cell_no].tk_x = tk_x; viscos_f = viscos_il_f = 1.0; @@ -1847,8 +1996,7 @@ fill_spec(int l_cell_no) //string_hsave(s_ptr2->name); sol_D[l_cell_no].spec[count_spec].z = s_ptr2->z; if (s_ptr2->dw == 0) - sol_D[l_cell_no].spec[count_spec].Dwt = - default_Dw * viscos_il_f; + sol_D[l_cell_no].spec[count_spec].Dwt = default_Dw * viscos_il_f; else { if (s_ptr2->dw_t) @@ -1860,19 +2008,92 @@ fill_spec(int l_cell_no) else sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr2->dw * viscos_il_f; } + + //if (implicit) // && name_ret.second && (l_cell_no > 1 || (l_cell_no == 1 && bcon_first != 2))) + //{ + // // name_ret = dif_spec_names.insert(s_ptr->name); // but not in implicit now... + // // must fill in the spec in previous cells, order the names, see below for aqueous species... + //} count_exch_spec++; count_spec++; } - continue; } lm = s_ptr->lm; + if (lm > MIN_LM) { + if (implicit) + { + name_ret = dif_spec_names.insert(s_ptr->name); + loc_spec_names.insert(s_ptr->name); + i2 = 0; + for (it = dif_spec_names.begin(); it != dif_spec_names.end(); it++) + { + if (*it == s_ptr->name) + break; + i2++; + } + if (i2 > count_spec) + { + if (l_cell_no == 0) + { + for (i1 = i2; i1 > count_spec; i1--) + { + it = dif_spec_names.find(s_ptr->name); + dif_spec_names.erase(--it); + } + } + else + { + // there are species before s_ptr->name that must be included first + i1 = 0; + for (it = loc_spec_names.begin(); it != loc_spec_names.end(); it++) + { + if (*it == s_ptr->name) + break; + i1++; + } + i3 = i2 - i1; + for (i1; i1 < i2; i1++) + { + if (i3 > sol_D[l_cell_no].spec_size) + { + sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(i3 + size_xt) * sizeof(struct spec)); + if (sol_D[l_cell_no].spec == NULL) + malloc_error(); + sol_D[l_cell_no].spec_size = i3 + size_xt; + } + memmove(&sol_D[l_cell_no].spec[i1], &sol_D[l_cell_no - 1].spec[i1], sizeof(struct spec)); + sol_D[l_cell_no].spec[i1].c = 0.0; + sol_D[l_cell_no].spec[i1].a = 0.0; + sol_D[l_cell_no].spec[i1].lm = 0.0; + sol_D[l_cell_no].spec[i1].lg = 0.0; + // may have to redefine Dwt appt + + loc_spec_names.insert(sol_D[l_cell_no].spec[i1].name); + count_spec++; + if (count_spec >= sol_D[l_cell_no].spec_size) + { + sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(count_spec + size_xt) * sizeof(struct spec)); + if (sol_D[l_cell_no].spec == NULL) + malloc_error(); + sol_D[l_cell_no].spec_size = count_spec + size_xt; + } + } + } + } + } + if (count_spec >= sol_D[l_cell_no].spec_size) + { + sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(count_spec + size_xt) * sizeof(struct spec)); + if (sol_D[l_cell_no].spec == NULL) + malloc_error(); + sol_D[l_cell_no].spec_size = count_spec + size_xt; + } sol_D[l_cell_no].spec[count_spec].name = s_ptr->name; sol_D[l_cell_no].spec[count_spec].type = AQ; - sol_D[l_cell_no].spec[count_spec].c = - s_ptr->moles / mass_water_aq_x; + sol_D[l_cell_no].spec[count_spec].c = s_ptr->moles / mass_water_aq_x; sol_D[l_cell_no].spec[count_spec].a = under(lm + s_ptr->lg); sol_D[l_cell_no].spec[count_spec].lm = lm; sol_D[l_cell_no].spec[count_spec].lg = s_ptr->lg; @@ -1899,39 +2120,89 @@ fill_spec(int l_cell_no) diffc_max = sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn); sol_D[l_cell_no].spec[count_spec].erm_ddl = s_ptr->erm_ddl; + if (implicit) + { + if (name_ret.second && l_cell_no) + { + // must fill in the spec in previous cells, order the names... + i3 = 0; + for (it = dif_spec_names.begin(); it != dif_spec_names.end(); ++it) + { + if (it == name_ret.first) + break; + i3++; + } + for (i1 = 0; i1 < l_cell_no; i1++) + { + i2 = sol_D[i1].count_spec + 1; + if (i2 > sol_D[i1].spec_size) + { + sol_D[i1].spec = (struct spec *) PHRQ_realloc(sol_D[i1].spec, (size_t)(i2 + size_xt) * sizeof(struct spec)); + if (sol_D[i1].spec == NULL) + malloc_error(); + sol_D[i1].spec_size = i2 + size_xt; + } + i2--; + for (i2; i2 > i3; i2--) + sol_D[i1].spec[i2] = sol_D[i1].spec[i2 - 1]; + + memmove(&sol_D[i1].spec[i2], &sol_D[l_cell_no].spec[i2], sizeof(struct spec)); + sol_D[i1].spec[i2].a = 0.0; + sol_D[i1].spec[i2].lm = 0.0; + sol_D[i1].spec[i2].lg = 0.0; + sol_D[i1].spec[i2].c = 0.0; + sol_D[i1].count_spec += 1; + } + } + } count_spec++; } } - sol_D[l_cell_no].spec = - (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, - (size_t) count_spec * - sizeof(struct spec)); - - if (sol_D[l_cell_no].spec == NULL) - malloc_error(); - { - if (count_spec > sol_D[l_cell_no].spec_size) - { - for (int i = sol_D[l_cell_no].spec_size; i < count_spec; i++) - { - sol_D[l_cell_no].spec[i].name = NULL; - sol_D[l_cell_no].spec[i].aq_name = NULL; - sol_D[l_cell_no].spec[i].type = -1; - sol_D[l_cell_no].spec[i].a = 0.0; - sol_D[l_cell_no].spec[i].lm = 0.0; - sol_D[l_cell_no].spec[i].lg = 0.0; - sol_D[l_cell_no].spec[i].c = 0.0; - sol_D[l_cell_no].spec[i].z = 0.0; - sol_D[l_cell_no].spec[i].Dwt = 0.0; - sol_D[l_cell_no].spec[i].dw_t = 0.0; - sol_D[l_cell_no].spec[i].erm_ddl = 0.0; - } - } - sol_D[l_cell_no].spec_size = count_spec; - } sol_D[l_cell_no].count_spec = count_spec; sol_D[l_cell_no].count_exch_spec = count_exch_spec; + if (implicit && loc_spec_names.size() < dif_spec_names.size()) + { + i3 = (int) dif_spec_names.size(); + if (i3 > sol_D[l_cell_no].spec_size) + { + sol_D[l_cell_no].spec = (struct spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (size_t)(i3 + size_xt) * sizeof(struct spec)); + if (sol_D[l_cell_no].spec == NULL) + malloc_error(); + sol_D[l_cell_no].spec_size = i3 + size_xt; + } + for (i1 = count_spec; i1 < i3; i1++) + { + memmove(&sol_D[l_cell_no].spec[i1], &sol_D[l_cell_no - 1].spec[i1], sizeof(struct spec)); + sol_D[l_cell_no].spec[i1].c = 0.0; + sol_D[l_cell_no].spec[i1].a = 0.0; + sol_D[l_cell_no].spec[i1].lm = 0.0; + sol_D[l_cell_no].spec[i1].lg = 0.0; + sol_D[l_cell_no].count_spec += 1; + loc_spec_names.insert(sol_D[l_cell_no].spec[i1].name); + count_spec++; + } + //if (loc_spec_names != dif_spec_names) + //{ + // error_string = sformatf( + // "Species in implicit diffusion in cell %d are %; different from previous cells with %d species.", + // l_cell_no, loc_spec_names.size(), dif_spec_names.size()); + // error_msg(error_string, CONTINUE); + //} + } + if (implicit && l_cell_no) + { + for (i = 0; i < count_spec; i++) + { + if (strcmp(sol_D[l_cell_no].spec[i].name, sol_D[l_cell_no - 1].spec[i].name)) + { + error_string = sformatf( + "Implicit diffusion: species %s in cell %d differs from species %s in previous cells.", + sol_D[l_cell_no].spec[i].name, l_cell_no, sol_D[l_cell_no - 1].spec[i].name); + error_msg(error_string, CONTINUE); + } + } + } return (OK); } @@ -1947,6 +2218,423 @@ sort_species_name(const void *ptr1, const void *ptr2) return (strcmp(nptr1->s->name, nptr2->s->name)); } +/* ---------------------------------------------------------------------- */ +void Phreeqc:: +diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) +/* ---------------------------------------------------------------------- */ +{ + // The structure comes from example 11.3 in A&P and Appelo 2017, CCR 101, 102: + // Ct2 is implicitly solved for the individual species, corrected for electro-neutrality, this always gives Ct2 > 0. + // Transport of aqueous species is summarized into master species. + // With electro-migration, transport of anions and cations is calculated in opposite directions since (sign) J = - z * dV. + // Only available moles are transported, thus are > 0, but if concentrations oscillate, + // change max_mixf in input file: -implicit true 1 # max_mixf = 1 (default). + int i, ifirst, ilast, cp, comp; + double mfr, mfr1, max_b = 0, b, grad, dVc, j_0e; + cxxSolution *sptr1, *sptr2; + + comp = sol_D[1].count_spec - sol_D[1].count_exch_spec; + if (heat_nmix) + comp += 1; + + if (y == NULL) + y = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); + if (y == NULL) malloc_error(); + if (Ct1 == NULL) + Ct1 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); + if (Ct1 == NULL) malloc_error(); + if (Ct2 == NULL) + Ct2 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); + if (Ct2 == NULL) malloc_error(); + if (l_tk_x2 == NULL) + l_tk_x2 = (LDBLE *) PHRQ_malloc((size_t) (count_cells + 2) * sizeof(LDBLE)); + if (l_tk_x2 == NULL) malloc_error(); + //for (i = 0; i < count_cells + 2; i++) + //{ + // y[i] = 0.0; Ct1[i] = 0.0; Ct2[i] = 0.0; l_tk_x2[i] = 0.0; + //} + + if (A == NULL) + { + A = (LDBLE **)PHRQ_malloc((size_t)(count_cells + 2) * sizeof(LDBLE *)); + if (A == NULL) malloc_error(); + for (i = 0; i < count_cells + 2; i++) + { + A[i] = (LDBLE *)PHRQ_malloc((size_t)3 * sizeof(LDBLE)); + if (A[i] == NULL) malloc_error(); + } + } + if (mixf == NULL) + mixf = (LDBLE **)PHRQ_malloc((size_t)(count_cells + 2) * sizeof(LDBLE *)); + if (mixf == NULL) malloc_error(); + if (comp + 2 > mixf_comp_size) + { + for (i = 0; i < count_cells + 2; i++) + { + mixf[i] = (LDBLE *)PHRQ_malloc((size_t)(comp + 2) * sizeof(LDBLE)); + if (mixf[i] == NULL) malloc_error(); + mixf_comp_size = comp + 2; + } + } + //for (i = 0; i < count_cells + 2; i++) + //{ + // for (int i1 = 0; i1 < 3; i1++) + // A[i][i1] = 0.0; + // for (int i1 = 0; i1 < comp; i1++) + // mixf[i][i1] = 0.0; + //} + + if (dV_dcell) + find_current = 1; + // obtain b_ij... + for (i = 0; i <= count_cells; i++) + { + if (!heat_nmix) + l_tk_x2[i] = (sol_D[i].tk_x + sol_D[i + 1].tk_x) / 2; + b = find_J(i, i + 1, 1, 1, stagnant); + if (b > max_b) + max_b = b; + } + if (dV_dcell) + find_current = 0; + + ifirst = (bcon_first == 2 ? 1 : 0); + ilast = (bcon_last == 2 ? count_cells - 1 : count_cells); + + for (cp = 0; cp < comp; cp++) + { + for (i = 0; i <= count_cells + 1; i++) + { + if (heat_nmix && cp == comp - 1) + Ct2[i] = Ct1[i] = sol_D[i].tk_x; + else + Ct2[i] = Ct1[i] = sol_D[i].spec[cp].c; + } + // calculate diffuson... + //for (int ss = 0; ss < substeps; ss++) + { + // fill coefficient matrix A ... + // boundary cells ... + if (heat_nmix && cp == comp - 1) + { + mfr = mixf[0][cp] = heat_mix_array[0]; + mfr1 = mixf[1][cp] = heat_mix_array[1]; + } + else + { + mfr = mixf[0][cp] = DDt * ct[0].v_m[cp].b_ij; // maybe add retardation appt + mfr1 = mixf[1][cp] = DDt * ct[1].v_m[cp].b_ij; + } + if (bcon_first == 2) + { + A[0][1] = 1; A[0][2] = 0; + A[1][0] = 0; A[1][1] = 1 + mfr1; A[1][2] = -mfr1; + } + else + { + if (dV_dcell) + { + A[0][1] = 1 + mfr; A[0][2] = -mfr; + } + else + { + A[0][1] = 1; A[0][2] = 0; + } + A[1][0] = -mfr; A[1][1] = 1 + mfr + mfr1; A[1][2] = -mfr1; + } + if (heat_nmix && cp == comp - 1) + { + mfr = mixf[count_cells - 1][cp] = heat_mix_array[count_cells - 1]; + mfr1 = mixf[count_cells][cp] = heat_mix_array[count_cells]; + } + else + { + mfr = mixf[count_cells - 1][cp] = DDt * ct[count_cells - 1].v_m[cp].b_ij; + mfr1 = mixf[count_cells][cp] = DDt * ct[count_cells].v_m[cp].b_ij; + } + if (bcon_last == 2) + { + A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr; A[count_cells][2] = 0; + A[count_cells + 1][0] = 0; A[count_cells + 1][1] = 1; + } + else + { + A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr + mfr1; A[count_cells][2] = -mfr1; + if (dV_dcell) + { + A[count_cells + 1][0] = -mfr1; A[count_cells + 1][1] = 1 + mfr1; + } + else + { + A[count_cells + 1][0] = 0; A[count_cells + 1][1] = 1; + } + } + + // inner cells ... + for (i = 2; i < count_cells; i++) + { + if (heat_nmix && cp == comp - 1) + { + mfr = mixf[i - 1][cp] = heat_mix_array[i - 1]; + mfr1 = mixf[i][cp] = heat_mix_array[i]; + } + else + { + mfr = mixf[i - 1][cp] = DDt * ct[i - 1].v_m[cp].b_ij; + mfr1 = mixf[i][cp] = DDt * ct[i].v_m[cp].b_ij; + } + A[i][0] = -mfr; A[i][1] = 1 + mfr + mfr1; A[i][2] = -mfr1; + } + // decompose A in LU : store L in A[..][0..1] and U in A[..][2] ... + for (i = 1; i <= count_cells + 1; i++) + { + A[i - 1][2] = A[i - 1][2] / A[i - 1][1]; + A[i][1] = A[i][1] - A[i][0] * A[i - 1][2]; + } + // solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1 ... + // First, find y in L.y = Ct1 + y[0] = Ct2[0]; + for (i = 1; i <= count_cells + 1; i++) + y[i] = (Ct2[i] - A[i][0] * y[i - 1]) / A[i][1]; + // Now obtain Ct2 in U.Ct2 = y ... + Ct2[count_cells + 1] = y[count_cells + 1]; + for (i = count_cells; i >= 0; i--) + Ct2[i] = y[i] - A[i][2] * Ct2[i + 1]; + } + // Moles transported by concentration gradient from cell [i] to [i + 1] go in tot1. + // Correct for electro-neutrality... + for (i = ifirst; i <= ilast + 1; i++) + { + if (heat_nmix && cp == comp - 1) + { + l_tk_x2[i] = (Ct2[i + 1] + Ct2[i]) / 2; + cell_data[i].temp = Ct2[i] - 273.15; + sptr1 = Utilities::Rxn_find(Rxn_solution_map, i); + sptr1->Set_tc(Ct2[i] - 273.15); + continue; + } + else if (i == ilast + 1) + continue; + + if (cp == 0) + current_cells[i].ele = current_cells[i].dif = 0; + grad = Ct2[i + 1] - Ct2[i]; + ct[i].J_ij[cp].name = sol_D[i].spec[cp].name; + ct[i].J_ij[cp].tot1 = -mixf[i][cp] * grad; + ct[i].J_ij[cp].charge = ct[i].v_m[cp].z; + if (ct[i].v_m[cp].z) + { + if (i == 0) + ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i + 1]; + else if (i == ilast) + ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i]; + else + ct[i].v_m[cp].zc = ct[i].v_m[cp].z * (Ct2[i] + Ct2[i + 1]) / 2; + + current_cells[i].dif -= ct[i].v_m[cp].z * mixf[i][cp] * grad; + current_cells[i].ele -= ct[i].v_m[cp].z * mixf[i][cp] * ct[i].v_m[cp].zc; + } + } + // define element list + if (cp == 0) + dif_els_names.clear(); + if (heat_nmix && cp == comp - 1) + { + comp -= 1; + break; + } + char * temp_name = string_duplicate(ct[1].J_ij[cp].name); + char * ptr = temp_name; + count_elts = 0; + get_elts_in_species(&ptr, 1); + free_check_null(temp_name); + for (int k = 0; k < count_elts; k++) + { + if (!strcmp(elt_list[k].elt->name, "X")) continue; + dif_els_names.insert(elt_list[k].elt->name); + } + } + count_m_s = (int) dif_els_names.size(); + sum_R = sum_Rd = 0; + for (i = ifirst; i <= ilast; i++) + { + ct[i].J_ij_count_spec = comp; + current_cells[i].R = l_tk_x2[i] / (current_cells[i].ele * F_Re3); + if (dV_dcell && !fix_current) + { + sum_R += current_cells[i].R; + if (i > 1) + sum_Rd += (current_cells[0].dif - current_cells[i].dif) * current_cells[i].R; + } + } + if (dV_dcell) + { + if (fix_current) + { + int sign = (current_x >= 0 ? 1 : -1); + current_x = sign * fix_current * DDt / F_C_MOL; + j_0e = current_x - current_cells[ifirst].dif; + } + else + { + j_0e = (dV_dcell * count_cells - sum_Rd) / sum_R; + current_x = j_0e + current_cells[ifirst].dif; + } + } + else + { + current_x = 1e-10 / F_C_MOL; + cell_data[ifirst].potV = 1e-8; + j_0e = current_x - current_cells[ifirst].dif; + } + dVc = j_0e * current_cells[ifirst].R; + cell_data[ifirst + 1].potV = cell_data[ifirst].potV + dVc; + for (i = ifirst + 1; i < ilast; i++) + { + dVc = current_cells[i].R * (current_x - current_cells[i].dif); + cell_data[i + 1].potV = cell_data[i].potV + dVc; + } + if (!dV_dcell || fix_current) + { + dVc = current_cells[i].R * (current_x - current_cells[i].dif); + cell_data[i + 1].potV = cell_data[i].potV + dVc; + } + + for (cp = 0; cp < comp; cp++) + { + for (i = ifirst; i <= ilast; i++) + { + ct[i].J_ij_sum = 0; + if (!ct[i].v_m[cp].z) + continue; + dVc = (cell_data[i + 1].potV - cell_data[i].potV) * F_Re3 / l_tk_x2[i]; + ct[i].J_ij[cp].tot1 -= mixf[i][cp] * ct[i].v_m[cp].zc * dVc; + // With dVc known, charge transport = 0 if dV_dcell = 0... + if (!dV_dcell) + ct[i].J_ij_sum += ct[i].v_m[cp].z * ct[i].J_ij[cp].tot1; + } + } + current_A = current_x / DDt * F_C_MOL; + + + for (i = ifirst; i <= ilast + 1; i++) + { + // preserve the potentials... + sptr1 = Utilities::Rxn_find(Rxn_solution_map, i); + sptr1->Set_potV(cell_data[i].potV); + if (i == ilast + 1) + continue; + + // Translate transport of the solute species into master species... + ct[i].count_m_s = count_m_s; + if (ct[i].m_s_size == 0 && ct[i].m_s != NULL) + ct[i].m_s = (struct M_S *) free_check_null(ct[i].m_s); + if (ct[i].m_s == NULL) + { + ct[i].m_s = (struct M_S *) PHRQ_malloc((size_t)(count_m_s) * sizeof(struct M_S)); + ct[i].m_s_size = count_m_s; + } + else if (count_m_s > ct[i].m_s_size) + { + ct[i].m_s = (struct M_S *) PHRQ_realloc(ct[i].m_s, (size_t)(count_m_s) * sizeof(struct M_S)); + ct[i].m_s_size = count_m_s; + } + if (ct[i].m_s == NULL) + malloc_error(); + std::set ::iterator it = dif_els_names.begin(); + for (int i1 = 0; i1 < count_m_s; i1++) + { + ct[i].m_s[i1].tot1 = 0.0; + ct[i].m_s[i1].charge = 0.0; + ct[i].m_s[i1].name = (*it).c_str(); + it++; + } + fill_m_s(ct[i].J_ij, ct[i].J_ij_count_spec, i); + } + /* + * 3. find the solutions, add or subtract the moles... + */ + std::map> ::iterator it1; + std::map ::iterator it2; + int icell, if1, il1, incr; + LDBLE dum1, dum2; + for (cp = 0; cp < count_m_s; cp++) + { + if (!dV_dcell || dV_dcell * ct[1].m_s[cp].charge <= 0) + { + if1 = ifirst; il1 = ilast + 1; incr = 1; + } + else + { + if1 = ilast; il1 = ifirst - 1; incr = -1; + } + for (icell = if1; icell != il1; icell += incr) + { + dum1 = dum2 = 0; + sptr1 = Utilities::Rxn_find(Rxn_solution_map, icell); + sptr2 = Utilities::Rxn_find(Rxn_solution_map, icell + 1); + if (!strcmp(ct[icell].m_s[cp].name, "H")) + { + if (dV_dcell || (icell > 0 && icell <= ilast)) + sptr1->Set_total_h(sptr1->Get_total_h() - ct[icell].m_s[cp].tot1); + if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last ==2)) + sptr2->Set_total_h(sptr2->Get_total_h() + ct[icell].m_s[cp].tot1); + continue; + } + if (!strcmp(ct[icell].m_s[cp].name, "O")) + { + if (dV_dcell || (icell > 0 && icell <= ilast)) + sptr1->Set_total_o(sptr1->Get_total_o() - ct[icell].m_s[cp].tot1); + if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) + sptr2->Set_total_o(sptr2->Get_total_o() + ct[icell].m_s[cp].tot1); + continue; + } + // see if icell - incr has negative moles, subtract it, so that it is not transported... + if (icell > 0 && icell <= ilast && + (it1 = neg_moles.find(icell - incr)) != neg_moles.end() + && (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end()) + { + ct[icell].m_s[cp].tot1 += it2->second; + // the negative moles are canceled... + neg_moles.erase(it1); + els.erase(it2); + neg_moles.insert(std::make_pair(icell - incr, els)); + } + dum1 = sptr1->Get_totals()[ct[icell].m_s[cp].name]; + dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name]; + { + if (ct[icell].m_s[cp].tot1 > dum1) + ct[icell].m_s[cp].tot1 = dum1; + else if (-ct[icell].m_s[cp].tot1 > dum2) + ct[icell].m_s[cp].tot1 = -dum2; + } + if (dV_dcell || (icell > 0 && icell <= ilast)) + { + dum1 -= ct[icell].m_s[cp].tot1; + sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : 1e-16); + } + if (dum1 < -1e-12 && incr > 0) + { + els.insert(std::make_pair(ct[icell].m_s[cp].name, dum1)); + neg_moles.insert(std::make_pair(icell, els)); + } + + if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) + { + dum2 += ct[icell].m_s[cp].tot1; + sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : 1e-16); + } + if (dum2 < -1e-12 && incr < 0) + { + els.insert(std::make_pair(ct[icell].m_s[cp].name, dum2)); + neg_moles.insert(std::make_pair(icell + 1, els)); + } + } + } + return; +} + + /* ---------------------------------------------------------------------- */ int Phreeqc:: multi_D(LDBLE DDt, int mobile_cell, int stagnant) @@ -2062,7 +2750,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant) /* * 1. obtain J_ij... */ - il_calcs = find_J(icell, jcell, mixf, DDt, stagnant); + il_calcs = (int) find_J(icell, jcell, mixf, DDt, stagnant); if (find_current) { @@ -2120,7 +2808,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant) } count_m_s = 0; } - fill_m_s(ct[icell].J_ij, ct[icell].J_ij_count_spec); + fill_m_s(ct[icell].J_ij, ct[icell].J_ij_count_spec, icell); /* * 3. find the solutions, add or subtract the moles... @@ -2305,13 +2993,14 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant) /* ---------------------------------------------------------------------- */ int Phreeqc:: -fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec) +fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell) /* ---------------------------------------------------------------------- */ { /* sum up the primary or secondary master_species from solute species * H and O go in tot1&2_h and tot1&2_o */ int j, k, l; + LDBLE fraction; char *ptr; for (j = 0; j < l_J_ij_count_spec; j++) @@ -2323,38 +3012,68 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec) get_elts_in_species(&ptr, 1); free_check_null(temp_name); } - for (k = 0; k < count_elts; k++) + if (implicit) { - if (strcmp(elt_list[k].elt->name, "X") == 0) - continue; - if (strcmp(elt_list[k].elt->name, "H") == 0) - { - tot1_h += elt_list[k].coef * l_J_ij[j].tot1; - tot2_h += elt_list[k].coef * l_J_ij[j].tot2; - } - else if (strcmp(elt_list[k].elt->name, "O") == 0) - { - tot1_o += elt_list[k].coef * l_J_ij[j].tot1; - tot2_o += elt_list[k].coef * l_J_ij[j].tot2; - } - else + for (k = 0; k < count_elts; k++) { for (l = 0; l < count_m_s; l++) { - if (strcmp(m_s[l].name, elt_list[k].elt->name) == 0) + if (strcmp(ct[icell].m_s[l].name, elt_list[k].elt->name) == 0) { - m_s[l].tot1 += elt_list[k].coef * l_J_ij[j].tot1; - m_s[l].tot2 += elt_list[k].coef * l_J_ij[j].tot2; + fraction = elt_list[k].coef * l_J_ij[j].tot1 + ct[icell].m_s[l].tot1; + if (fraction) + fraction = elt_list[k].coef * l_J_ij[j].tot1 / fraction; + else if (!ct[icell].m_s[l].tot1) + fraction = 1; + ct[icell].m_s[l].tot1 += elt_list[k].coef * l_J_ij[j].tot1; + if (!strcmp(elt_list[k].elt->name, "H") || !strcmp(elt_list[k].elt->name, "O")) + break; + ct[icell].m_s[l].charge *= (1 - fraction); + ct[icell].m_s[l].charge += fraction * l_J_ij[j].charge; break; } } - if (l == count_m_s) + } + } + else + { + for (k = 0; k < count_elts; k++) + { + if (strcmp(elt_list[k].elt->name, "X") == 0) + continue; + if (strcmp(elt_list[k].elt->name, "H") == 0) { - //m_s[l].name = string_hsave(elt_list[k].elt->name); - m_s[l].name = elt_list[k].elt->name; - m_s[l].tot1 = elt_list[k].coef * l_J_ij[j].tot1; - m_s[l].tot2 = elt_list[k].coef * l_J_ij[j].tot2; - count_m_s++; + tot1_h += elt_list[k].coef * l_J_ij[j].tot1; + tot2_h += elt_list[k].coef * l_J_ij[j].tot2; + } + else if (strcmp(elt_list[k].elt->name, "O") == 0) + { + tot1_o += elt_list[k].coef * l_J_ij[j].tot1; + tot2_o += elt_list[k].coef * l_J_ij[j].tot2; + } + else + { + for (l = 0; l < count_m_s; l++) + { + if (strcmp(m_s[l].name, elt_list[k].elt->name) == 0) + { + fraction = elt_list[k].coef * l_J_ij[j].tot1 / (elt_list[k].coef * l_J_ij[j].tot1 + m_s[l].tot1); + m_s[l].tot1 += elt_list[k].coef * l_J_ij[j].tot1; + m_s[l].tot2 += elt_list[k].coef * l_J_ij[j].tot2; + m_s[l].charge *= (1 - fraction); + m_s[l].charge += fraction * l_J_ij[j].charge; + break; + } + } + if (l == count_m_s) + { + //m_s[l].name = string_hsave(elt_list[k].elt->name); + m_s[l].name = elt_list[k].elt->name; + m_s[l].tot1 = elt_list[k].coef * l_J_ij[j].tot1; + m_s[l].tot2 = elt_list[k].coef * l_J_ij[j].tot2; + m_s[l].charge = l_J_ij[j].charge; + count_m_s++; + } } } } @@ -2362,7 +3081,7 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec) return (OK); } /* ---------------------------------------------------------------------- */ -int Phreeqc:: +LDBLE Phreeqc:: find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) /* ---------------------------------------------------------------------- */ { @@ -2407,11 +3126,12 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) int i, i_max, j, j_max, k, k_il, only_counter, il_calcs; int i1; LDBLE A1 = 0.0, A2 = 0.0, ddlm, aq1, aq2, t_aq1, t_aq2, f_free_i, f_free_j; - LDBLE dl_aq1, dl_aq2, c_dl, dum, dum1, dum2, tort1, tort2, b_i, b_j; + LDBLE dl_aq1, dl_aq2, dum, dum1, dum2, tort1, tort2, b_i, b_j; LDBLE Sum_zM, aq_il1, aq_il2; LDBLE por_il1, por_il2, por_il12 = 0.0; LDBLE cec1, cec2, cec12 = 0.0, rc1 = 0.0, rc2 = 0.0; LDBLE dV, c1, c2; + LDBLE max_b = 0.0; // appt il_calcs = (interlayer_Dflag ? 1 : 0); cxxSurface *s_ptr1, *s_ptr2; @@ -2435,9 +3155,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) { if (stagnant) { - if (cell_data[icell].por < multi_Dpor_lim - || cell_data[jcell].por < multi_Dpor_lim) - return (il_calcs); + if (cell_data[icell].por < multi_Dpor_lim || cell_data[jcell].por < multi_Dpor_lim) + { + if (!implicit) + return (il_calcs); + } } else { /* regular column... */ @@ -2454,14 +3176,15 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) sum_R += current_cells[icell].R; sum_Rd += current_cells[0].dif * current_cells[icell].R; } - return (il_calcs); + if (!implicit) + return (il_calcs); } } } /* do the calcs */ - aq1 = Utilities::Rxn_find(Rxn_solution_map, icell)->Get_mass_water(); - aq2 = Utilities::Rxn_find(Rxn_solution_map, jcell)->Get_mass_water(); + aq1 = ct[icell].kgw = Utilities::Rxn_find(Rxn_solution_map, icell)->Get_mass_water(); + aq2 = ct[jcell].kgw = Utilities::Rxn_find(Rxn_solution_map, jcell)->Get_mass_water(); /* * check if DL calculations must be made, find amounts of water... */ @@ -2686,13 +3409,29 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) */ k = sol_D[icell].count_spec + sol_D[jcell].count_spec; - ct[icell].J_ij = (struct J_ij *) free_check_null(ct[icell].J_ij); - ct[icell].J_ij = (struct J_ij *) PHRQ_malloc((size_t) k * sizeof(struct J_ij)); + if (ct[icell].J_ij == NULL) + { + ct[icell].J_ij = (struct J_ij *) PHRQ_malloc((size_t)(k) * sizeof(struct J_ij)); + ct[icell].J_ij_size = k; + } + else if (k > ct[icell].J_ij_size) + { + ct[icell].J_ij = (struct J_ij *) PHRQ_realloc(ct[icell].J_ij, (size_t)(k) * sizeof(struct J_ij)); + ct[icell].J_ij_size = k; + } if (ct[icell].J_ij == NULL) malloc_error(); - ct[icell].v_m = (struct V_M *) free_check_null(ct[icell].v_m); - ct[icell].v_m = (struct V_M *) PHRQ_malloc((size_t) k * sizeof(struct V_M)); + if (ct[icell].v_m == NULL) + { + ct[icell].v_m = (struct V_M *) PHRQ_malloc((size_t)(k) * sizeof(struct V_M)); + ct[icell].v_m_size = k; + } + else if (k > ct[icell].v_m_size) + { + ct[icell].v_m = (struct V_M *) PHRQ_realloc(ct[icell].v_m, (size_t)(k) * sizeof(struct V_M)); + ct[icell].v_m_size = k; + } if (ct[icell].v_m == NULL) malloc_error(); @@ -2708,7 +3447,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) //ct[icell].v_m[i].Dzc = 0.0; ct[icell].v_m[i].b_ij = 0.0; } - ct[icell].Dz2c = ct[icell].Dz2c_dl = ct[icell].Dz2c_il = 0.0; + ct[icell].Dz2c = ct[icell].Dz2c_il = 0.0; if (il_calcs) { @@ -2747,6 +3486,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) while (i < i_max || j < j_max) { + // appt debug + //ct[icell].J_ij[i].name = sol_D[icell].spec[i].name; + //ct[icell].J_ij[i+1].name = sol_D[jcell].spec[j].name; if (j == j_max || (i < i_max && strcmp(sol_D[icell].spec[i].name, sol_D[jcell].spec[j].name) < 0)) { @@ -3006,7 +3748,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) if (ct[icell].v_m[k].z) ct[icell].v_m[k].zc = ct[icell].v_m[k].z * ct[icell].v_m[k].c; - if (dV_dcell && ct[icell].v_m[k].z && !fix_current) + if (dV_dcell && ct[icell].v_m[k].z && !fix_current && !implicit) // appt { // compare diffuse and electromotive forces dum = ct[icell].v_m[k].grad; @@ -3069,6 +3811,9 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) else if (icell == count_cells) ct[icell].v_m[k].b_ij = b_i; } + if (ct[icell].v_m[k].b_ij > max_b) + max_b = ct[icell].v_m[k].b_ij; + if (ct[icell].v_m[k].z) ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z; @@ -3084,6 +3829,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) j++; } } + if (implicit) + return(max_b); // appt /* * fill in J_ij... */ @@ -3097,17 +3844,16 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) //if (transport_step >= 100) // debug... // icell = icell; current_cells[icell].ele = current_cells[icell].dif = 0; - dum = dV_dcell * F_Re3 / tk_x2; for (i = 0; i < ct[icell].J_ij_count_spec; i++) { if (!ct[icell].v_m[i].z) continue; current_cells[icell].ele -= ct[icell].v_m[i].b_ij * ct[icell].v_m[i].z * - ct[icell].v_m[i].zc * dum; + ct[icell].v_m[i].zc; current_cells[icell].dif -= ct[icell].v_m[i].b_ij * ct[icell].v_m[i].z * ct[icell].v_m[i].grad; } - current_cells[icell].R = dV_dcell / current_cells[icell].ele; + current_cells[icell].R = tk_x2 / (current_cells[icell].ele * F_Re3); sum_R += current_cells[icell].R; sum_Rd += (current_cells[0].dif - current_cells[icell].dif) * current_cells[icell].R; return(il_calcs); @@ -3118,7 +3864,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) dV_dcell2: ct[icell].J_ij_sum = 0; - Sum_zM = c_dl = 0.0; + Sum_zM = 0.0; for (i = 0; i < ct[icell].J_ij_count_spec; i++) { @@ -3128,6 +3874,7 @@ dV_dcell2: for (i = 0; i < ct[icell].J_ij_count_spec; i++) { ct[icell].J_ij[i].tot1 = -ct[icell].v_m[i].grad; + ct[icell].J_ij[i].charge = ct[icell].v_m[i].z; if (!dV_dcell && ct[icell].v_m[i].z && ct[icell].Dz2c > 0) ct[icell].J_ij[i].tot1 += Sum_zM * ct[icell].v_m[i].zc / ct[icell].Dz2c; if (stagnant) @@ -3143,6 +3890,7 @@ dV_dcell2: if (dV_dcell) { + ct[icell].J_ij_sum = 0; dV = cell_data[jcell].potV - cell_data[icell].potV; dum = dV * F_Re3 / tk_x2; // perhaps adapt dV for getting equal current... @@ -3167,9 +3915,15 @@ dV_dcell2: ct[icell].J_ij[i].tot1 -= ct[icell].v_m[i].b_ij * ct[icell].v_m[i].zc * dum * DDt; ct[icell].J_ij[i].tot2 = ct[icell].J_ij[i].tot1; + ct[icell].J_ij_sum += ct[icell].v_m[i].z * ct[icell].J_ij[i].tot1; } current_A = current_x * F_C_MOL; } + if (abs(ct[icell].J_ij_sum / DDt - current_x) > 1e-13 && (dV_dcell || (icell > 0 || jcell >= count_cells))) // appt + { + sprintf(token, "Moles charge added in cell %d is: %.3e.\n", icell, ct[icell].J_ij_sum / DDt - current_x); + screen_msg(token); + } /* * calculate interlayer mass transfer... */ @@ -3191,6 +3945,7 @@ dV_dcell2: ct[icell].J_ij_il[i].tot1 *= ct[icell].A_ij_il * DDt; ct[icell].J_ij_sum += ct[icell].v_m_il[i].z * ct[icell].J_ij_il[i].tot1; ct[icell].J_ij_il[i].tot2 = ct[icell].J_ij_il[i].tot1; + ct[icell].J_ij_il[i].charge = ct[icell].v_m_il[i].z; } /* express the transfer in elemental moles... */ @@ -3202,12 +3957,13 @@ dV_dcell2: malloc_error(); for (i1 = 0; i1 < count_elements; i1++) { + m_s[i1].charge = 0; m_s[i1].name = NULL; m_s[i1].tot1 = 0; m_s[i1].tot2 = 0; } count_m_s = 0; - fill_m_s(ct[icell].J_ij_il, k_il); + fill_m_s(ct[icell].J_ij_il, k_il, icell); /* do the mass transfer... */ if (icell > 0 || stagnant) From 55ea163dd0fab53e022370f78fa3364e5e9f29b7 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Thu, 30 May 2019 22:03:38 -0600 Subject: [PATCH 2/4] Implicit seems to be working with Tony's latest changes --- Phreeqc.h | 2 +- basicsubs.cpp | 3 +- class_main.cpp | 10 +- common/PHRQ_io.cpp | 4 +- pitzer.cpp | 3 - tidy.cpp | 2 +- transport.cpp | 483 ++++++++++++++++++++++++++++----------------- utilities.cpp | 3 +- 8 files changed, 313 insertions(+), 197 deletions(-) diff --git a/Phreeqc.h b/Phreeqc.h index f1f03b64..96df04da 100644 --- a/Phreeqc.h +++ b/Phreeqc.h @@ -1069,7 +1069,7 @@ public: void diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant); int fill_spec(int cell_no); void define_ct_structures(void); - int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec, int i); + int fill_m_s(struct J_ij *J_ij, int J_ij_count_spec, int i, int stagnant); static int sort_species_name(const void *ptr1, const void *ptr2); int disp_surf(LDBLE stagkin_time); int diff_stag_surf(int mobile_cell); diff --git a/basicsubs.cpp b/basicsubs.cpp index c48d56b7..17ae7a83 100644 --- a/basicsubs.cpp +++ b/basicsubs.cpp @@ -1470,7 +1470,8 @@ kinetics_moles(const char *kinetics_name) error_string = sformatf( "No data for rate %s in KINETICS keyword.", kinetics_name); - warning_msg(error_string); + //if (count_warnings >= 0) // appt debug cvode + // warning_msg(error_string); return (0); } /* ---------------------------------------------------------------------- */ diff --git a/class_main.cpp b/class_main.cpp index f6477e02..9ad59d8f 100644 --- a/class_main.cpp +++ b/class_main.cpp @@ -277,8 +277,8 @@ write_banner(void) /* version */ #ifdef NPP - //len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8 AmpŠre"); - len = sprintf(buffer, "* PHREEQC-%s *", "3.4.8"); + //len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1 AmpŠre"); + len = sprintf(buffer, "* PHREEQC-%s *", "3.5.1"); #else len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@"); #endif @@ -302,7 +302,7 @@ write_banner(void) /* date */ #ifdef NPP - len = sprintf(buffer, "%s", "January 17, 2019"); + len = sprintf(buffer, "%s", "May 29, 2019"); #else len = sprintf(buffer, "%s", "@VER_DATE@"); #endif @@ -492,12 +492,12 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie, screen_msg(sformatf("Database file: %s\n\n", token)); strcpy(db_file, token); #ifdef NPP - //output_msg(sformatf("Using PHREEQC: version 3.4.8 Ampère, compiled January 17, 2019\n")); + //output_msg(sformatf("Using PHREEQC: version 3.5.1 Ampère, compiled May 29, 2019\n")); #endif output_msg(sformatf(" Input file: %s\n", in_file)); output_msg(sformatf(" Output file: %s\n", out_file)); #ifdef NPP - output_msg(sformatf("Using PHREEQC: version 3.4.8, compiled January 17, 2019\n")); + output_msg(sformatf("Using PHREEQC: version 3.5.1, compiled May 29, 2019\n")); #endif output_msg(sformatf("Database file: %s\n\n", token)); #ifdef NPP diff --git a/common/PHRQ_io.cpp b/common/PHRQ_io.cpp index ab34f8c1..c57c4b7f 100644 --- a/common/PHRQ_io.cpp +++ b/common/PHRQ_io.cpp @@ -32,7 +32,7 @@ PHRQ_io(void) punch_on = true; error_on = true; dump_on = true; - echo_on = true; //**appt + echo_on = true; screen_on = true; echo_destination = ECHO_OUTPUT; @@ -800,7 +800,7 @@ get_line(void) this->push_istream(next_stream); std::ostringstream errstr; errstr << "\n\tReading data from " << file_name <<" ...\n"; - output_msg(errstr.str().c_str()); // **appt + output_msg(errstr.str().c_str()); } continue; } diff --git a/pitzer.cpp b/pitzer.cpp index 90ce1bd8..c911633e 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -2555,7 +2555,6 @@ gammas_pz(bool exch_a_f) case 9: /* activity water */ break; case 4: /* Exchange */ - /* * Find CEC * z contains valence of cation for exchange species, alk contains cec @@ -2580,7 +2579,6 @@ gammas_pz(bool exch_a_f) /* * All other species */ - /* modific 29 july 2005... */ if (s_x[i]->equiv != 0 && s_x[i]->alk > 0) { @@ -2595,7 +2593,6 @@ gammas_pz(bool exch_a_f) continue; coef = s_x[i]->rxn_x->token[j].coef; s_x[i]->lg += coef * s_x[i]->rxn_x->token[j].s->lg; - /*s_x[i]->dg += coef * s_x[i]->rxn_x->token[j].s->dg;*//* remove? */ } } if (s_x[i]->a_f && s_x[i]->primary == NULL && s_x[i]->moles) diff --git a/tidy.cpp b/tidy.cpp index f3b15be8..09502d18 100644 --- a/tidy.cpp +++ b/tidy.cpp @@ -1437,7 +1437,7 @@ tidy_inverse(void) s_eminus->primary->in = TRUE; /* Include electrons */ if (master_alk_ptr) { - master_alk_ptr->in = TRUE; /* Include alkalinity */ + master_alk_ptr->in = TRUE; /* Include alkalinity */ } else { diff --git a/transport.cpp b/transport.cpp index 2230b10a..75413c97 100644 --- a/transport.cpp +++ b/transport.cpp @@ -136,6 +136,7 @@ transport(void) moles_added = (struct MOLES_ADDED *) PHRQ_malloc((size_t) (count_elements)* sizeof(struct MOLES_ADDED)); if (moles_added == NULL) malloc_error(); + for (i = 0; i < count_elements; i++) { moles_added[i].name = NULL; @@ -283,9 +284,9 @@ transport(void) transport_step = 0; for (i = 0; i <= count_cells + 1; i++) { - //if ((bcon_first == 2 && i == 0) || - // (bcon_last == 2 && i == count_cells + 1)) - // continue; + if (!implicit && ((bcon_first == 2 && i == 0) || + (bcon_last == 2 && i == count_cells + 1))) + continue; set_initial_moles(i); cell_no = i; set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0); @@ -332,23 +333,18 @@ transport(void) * Initialize mixing factors, define kinetics times * for multicomponent diffusion, limit mixing by diffc_max (usually from H+) */ - if (multi_Dflag == TRUE) + if (multi_Dflag) diffc_tr = diffc_max; if ((stag_data->exch_f > 0) && (stag_data->count_stag == 1)) { /* multi_D calc's are OK if all cells have the same amount of water */ - /* if (multi_Dflag == TRUE) + if (multi_Dflag == TRUE) { - sprintf(token, "multi_D calc's and stagnant: define MIXing factors explicitly, or \n\t give all cells the same amount of water."); - warning_msg(token); + sprintf(token, "multi_D calc's and stagnant: define MIXing factors explicitly, or \n\t adapt the default Dw of -multi_d to match the mobile-immobile exchange factor."); + warning_msg(token); } - */ - + Rxn_mix_map.clear(); - - /* - * stagnant mix factors go in mix[0 .. count_cells] - */ } /* * mix[] is extended in init_mix(), to accommodate column mix factors @@ -484,7 +480,7 @@ transport(void) */ for (i = 0; i <= count_cells + 1; i++) { - if (!dV_dcell && (i == 0 || i == count_cells + 1)) + if (!dV_dcell && !implicit && (i == 0 || i == count_cells + 1)) continue; set_initial_moles(i); } @@ -555,8 +551,14 @@ transport(void) for (i = 0; i <= count_cells + 1; i++) { - if (!dV_dcell && (i == 0 || i == count_cells + 1 && !implicit)) + if (!dV_dcell && (i == 0 || i == count_cells + 1) && !implicit) + { + if (j == nmix && stag_data->count_stag == 0 && + (cell_data[0].print || cell_data[0].punch || + cell_data[count_cells + 1].print || cell_data[count_cells + 1].punch)) + print_punch(i, false); continue; + } if (overall_iterations > max_iter) max_iter = overall_iterations; cell_no = i; @@ -571,10 +573,9 @@ transport(void) transport_step, j, i, max_iter); status(0, token); - if (implicit || i == 0 || i == count_cells + 1) + if (i == 0 || i == count_cells + 1) { run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i - //print_punch(i, true); // appt debug } else run_reactions(i, kin_time, DISP, step_fraction); // nsaver = -2 @@ -584,7 +585,7 @@ transport(void) /* punch and output file */ if (ishift == 0 && j == nmix && stag_data->count_stag == 0) print_punch(i, true); - if (!implicit && i > 1) + if (i > 1) Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1); saver(); @@ -610,7 +611,7 @@ transport(void) } } - if (!dV_dcell && !implicit) + if (!dV_dcell) Utilities::Rxn_copy(Rxn_solution_map, -2, count_cells); /* Stagnant zone mixing after completion of each diffusive/dispersive step ... */ @@ -722,7 +723,7 @@ transport(void) /* * thermal diffusion when nmix = 0... */ - if ((nmix == 0) && (heat_nmix > 0)) // must add disp in implicit appt ?? + if (nmix == 0 && heat_nmix > 0 && !implicit) { heat_mix(heat_nmix); /* equilibrate again ... */ @@ -857,7 +858,7 @@ transport(void) /* for each cell in column */ for (i = 0; i <= count_cells + 1; i++) { - if (!dV_dcell && (i == 0 || i == count_cells + 1)) + if (!dV_dcell && (i == 0 || i == count_cells + 1) && !implicit) { if (j == nmix && stag_data->count_stag == 0 && (cell_data[0].print || cell_data[0].punch || @@ -879,7 +880,7 @@ transport(void) transport_step, j, i, max_iter); status(0, token); - if (implicit || i == 0 || i == count_cells + 1) + if (i == 0 || i == count_cells + 1) run_reactions(i, kin_time, NOMIX, step_fraction); else run_reactions(i, kin_time, DISP, step_fraction); @@ -887,7 +888,7 @@ transport(void) fill_spec(i); if (j == nmix && stag_data->count_stag == 0) print_punch(i, true); - if (!implicit && i > 1) + if (i > 1) Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1); saver(); @@ -1217,7 +1218,7 @@ init_mix(void) if (maxmix > max_mixf) l_nmix = 1 + (int)floor(maxmix / max_mixf); if (ishift != 0 && (bcon_first == 1 || bcon_last == 1)) - { // must be adjusted appt + { if (l_nmix < 2) l_nmix = 2; } @@ -1235,9 +1236,9 @@ init_mix(void) error_msg(token, STOP); } if (bcon_first == 1 || bcon_last == 1) - l_nmix = 1 + (int) floor(2.25 * maxmix); + l_nmix = 1 + (int)floor(2.25 * maxmix); else - l_nmix = 1 + (int) floor(1.5 * maxmix); + l_nmix = 1 + (int)floor(1.5 * maxmix); if (ishift != 0 && (bcon_first == 1 || bcon_last == 1)) { @@ -1245,24 +1246,24 @@ init_mix(void) l_nmix = 2; } if (mcd_substeps > 1) - l_nmix = (int) ceil(l_nmix * mcd_substeps); + l_nmix = (int)ceil(l_nmix * mcd_substeps); + } - for (i = 1; i <= count_cells; i++) - { - m[i] /= l_nmix; - m1[i] /= l_nmix; - /* - * Fill mix structure - */ - cxxMix temp_mix; - temp_mix.Set_n_user(i); - temp_mix.Set_n_user_end(i); + for (i = 1; i <= count_cells; i++) + { + m[i] /= l_nmix; + m1[i] /= l_nmix; + /* + * Fill mix structure + */ + cxxMix temp_mix; + temp_mix.Set_n_user(i); + temp_mix.Set_n_user_end(i); - temp_mix.Add(i - 1, m[i]); - temp_mix.Add(i + 1, m1[i]); - temp_mix.Add(i, 1.0 - m[i] - m1[i]); - Dispersion_mix_map[i] = temp_mix; - } + temp_mix.Add(i - 1, m[i]); + temp_mix.Add(i + 1, m1[i]); + temp_mix.Add(i, 1.0 - m[i] - m1[i]); + Dispersion_mix_map[i] = temp_mix; } m = (LDBLE *)free_check_null(m); m1 = (LDBLE *)free_check_null(m1); @@ -1609,13 +1610,12 @@ init_heat_mix(int l_nmix) corr_disp = 1.; if (correct_disp == TRUE && ishift != 0) { + mixf = (l_nmix > 1 ? l_nmix : 1); if (bcon_first == 3) - corr_disp += 1. / count_cells; + corr_disp += 1. / count_cells / mixf; if (bcon_last == 3) - corr_disp += 1. / count_cells; + corr_disp += 1. / count_cells / mixf; } - if (l_nmix > 0 && !implicit) - corr_disp /= l_nmix; maxmix = 0.0; for (i = 1; i < count_cells; i++) { @@ -1847,6 +1847,8 @@ fill_spec(int l_cell_no) // for implicit std::set ::iterator it; std::pair ::iterator, bool> name_ret; + if (implicit && !l_cell_no) + dif_spec_names.clear(); size_xt = 5; //sol_D[l_cell_no].spec = (struct spec *) free_check_null(sol_D[l_cell_no].spec); @@ -1869,8 +1871,8 @@ fill_spec(int l_cell_no) sol_D[l_cell_no].spec[i].aq_name = NULL; sol_D[l_cell_no].spec[i].type = -1; sol_D[l_cell_no].spec[i].a = 0.0; - sol_D[l_cell_no].spec[i].lm = 0.0; - sol_D[l_cell_no].spec[i].lg = 0.0; + sol_D[l_cell_no].spec[i].lm = -3.0; + sol_D[l_cell_no].spec[i].lg = -0.04; sol_D[l_cell_no].spec[i].c = 0.0; sol_D[l_cell_no].spec[i].z = 0.0; sol_D[l_cell_no].spec[i].Dwt = 0.0; @@ -2017,6 +2019,7 @@ fill_spec(int l_cell_no) count_exch_spec++; count_spec++; } + continue; } lm = s_ptr->lm; @@ -2067,9 +2070,8 @@ fill_spec(int l_cell_no) memmove(&sol_D[l_cell_no].spec[i1], &sol_D[l_cell_no - 1].spec[i1], sizeof(struct spec)); sol_D[l_cell_no].spec[i1].c = 0.0; sol_D[l_cell_no].spec[i1].a = 0.0; - sol_D[l_cell_no].spec[i1].lm = 0.0; - sol_D[l_cell_no].spec[i1].lg = 0.0; - // may have to redefine Dwt appt + sol_D[l_cell_no].spec[i1].lm = -3.0; + sol_D[l_cell_no].spec[i1].lg = -0.04; loc_spec_names.insert(sol_D[l_cell_no].spec[i1].name); count_spec++; @@ -2148,8 +2150,8 @@ fill_spec(int l_cell_no) memmove(&sol_D[i1].spec[i2], &sol_D[l_cell_no].spec[i2], sizeof(struct spec)); sol_D[i1].spec[i2].a = 0.0; - sol_D[i1].spec[i2].lm = 0.0; - sol_D[i1].spec[i2].lg = 0.0; + sol_D[i1].spec[i2].lm = -3.0; + sol_D[i1].spec[i2].lg = -0.04; sol_D[i1].spec[i2].c = 0.0; sol_D[i1].count_spec += 1; } @@ -2175,8 +2177,8 @@ fill_spec(int l_cell_no) memmove(&sol_D[l_cell_no].spec[i1], &sol_D[l_cell_no - 1].spec[i1], sizeof(struct spec)); sol_D[l_cell_no].spec[i1].c = 0.0; sol_D[l_cell_no].spec[i1].a = 0.0; - sol_D[l_cell_no].spec[i1].lm = 0.0; - sol_D[l_cell_no].spec[i1].lg = 0.0; + sol_D[l_cell_no].spec[i1].lm = -3.0; + sol_D[l_cell_no].spec[i1].lg = -0.04; sol_D[l_cell_no].count_spec += 1; loc_spec_names.insert(sol_D[l_cell_no].spec[i1].name); @@ -2310,97 +2312,93 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) else Ct2[i] = Ct1[i] = sol_D[i].spec[cp].c; } - // calculate diffuson... - //for (int ss = 0; ss < substeps; ss++) + // fill coefficient matrix A ... + // boundary cells ... + if (heat_nmix && cp == comp - 1) { - // fill coefficient matrix A ... - // boundary cells ... - if (heat_nmix && cp == comp - 1) + mfr = mixf[0][cp] = heat_mix_array[0]; + mfr1 = mixf[1][cp] = heat_mix_array[1]; + } + else + { + mfr = mixf[0][cp] = DDt * ct[0].v_m[cp].b_ij; // maybe add retardation appt + mfr1 = mixf[1][cp] = DDt * ct[1].v_m[cp].b_ij; + } + if (bcon_first == 2) + { + A[0][1] = 1; A[0][2] = 0; + A[1][0] = 0; A[1][1] = 1 + mfr1; A[1][2] = -mfr1; + } + else + { + if (dV_dcell) { - mfr = mixf[0][cp] = heat_mix_array[0]; - mfr1 = mixf[1][cp] = heat_mix_array[1]; + A[0][1] = 1 + mfr; A[0][2] = -mfr; } else - { - mfr = mixf[0][cp] = DDt * ct[0].v_m[cp].b_ij; // maybe add retardation appt - mfr1 = mixf[1][cp] = DDt * ct[1].v_m[cp].b_ij; - } - if (bcon_first == 2) { A[0][1] = 1; A[0][2] = 0; - A[1][0] = 0; A[1][1] = 1 + mfr1; A[1][2] = -mfr1; + } + A[1][0] = -mfr; A[1][1] = 1 + mfr + mfr1; A[1][2] = -mfr1; + } + if (heat_nmix && cp == comp - 1) + { + mfr = mixf[count_cells - 1][cp] = heat_mix_array[count_cells - 1]; + mfr1 = mixf[count_cells][cp] = heat_mix_array[count_cells]; + } + else + { + mfr = mixf[count_cells - 1][cp] = DDt * ct[count_cells - 1].v_m[cp].b_ij; + mfr1 = mixf[count_cells][cp] = DDt * ct[count_cells].v_m[cp].b_ij; + } + if (bcon_last == 2) + { + A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr; A[count_cells][2] = 0; + A[count_cells + 1][0] = 0; A[count_cells + 1][1] = 1; + } + else + { + A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr + mfr1; A[count_cells][2] = -mfr1; + if (dV_dcell) + { + A[count_cells + 1][0] = -mfr1; A[count_cells + 1][1] = 1 + mfr1; } else { - if (dV_dcell) - { - A[0][1] = 1 + mfr; A[0][2] = -mfr; - } - else - { - A[0][1] = 1; A[0][2] = 0; - } - A[1][0] = -mfr; A[1][1] = 1 + mfr + mfr1; A[1][2] = -mfr1; - } - if (heat_nmix && cp == comp - 1) - { - mfr = mixf[count_cells - 1][cp] = heat_mix_array[count_cells - 1]; - mfr1 = mixf[count_cells][cp] = heat_mix_array[count_cells]; - } - else - { - mfr = mixf[count_cells - 1][cp] = DDt * ct[count_cells - 1].v_m[cp].b_ij; - mfr1 = mixf[count_cells][cp] = DDt * ct[count_cells].v_m[cp].b_ij; - } - if (bcon_last == 2) - { - A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr; A[count_cells][2] = 0; A[count_cells + 1][0] = 0; A[count_cells + 1][1] = 1; } + } + // inner cells ... + for (i = 2; i < count_cells; i++) + { + if (heat_nmix && cp == comp - 1) + { + mfr = mixf[i - 1][cp] = heat_mix_array[i - 1]; + mfr1 = mixf[i][cp] = heat_mix_array[i]; + } else { - A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr + mfr1; A[count_cells][2] = -mfr1; - if (dV_dcell) - { - A[count_cells + 1][0] = -mfr1; A[count_cells + 1][1] = 1 + mfr1; - } - else - { - A[count_cells + 1][0] = 0; A[count_cells + 1][1] = 1; - } + mfr = mixf[i - 1][cp] = DDt * ct[i - 1].v_m[cp].b_ij; + mfr1 = mixf[i][cp] = DDt * ct[i].v_m[cp].b_ij; } - - // inner cells ... - for (i = 2; i < count_cells; i++) - { - if (heat_nmix && cp == comp - 1) - { - mfr = mixf[i - 1][cp] = heat_mix_array[i - 1]; - mfr1 = mixf[i][cp] = heat_mix_array[i]; - } - else - { - mfr = mixf[i - 1][cp] = DDt * ct[i - 1].v_m[cp].b_ij; - mfr1 = mixf[i][cp] = DDt * ct[i].v_m[cp].b_ij; - } - A[i][0] = -mfr; A[i][1] = 1 + mfr + mfr1; A[i][2] = -mfr1; - } - // decompose A in LU : store L in A[..][0..1] and U in A[..][2] ... - for (i = 1; i <= count_cells + 1; i++) - { - A[i - 1][2] = A[i - 1][2] / A[i - 1][1]; - A[i][1] = A[i][1] - A[i][0] * A[i - 1][2]; - } - // solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1 ... - // First, find y in L.y = Ct1 - y[0] = Ct2[0]; - for (i = 1; i <= count_cells + 1; i++) - y[i] = (Ct2[i] - A[i][0] * y[i - 1]) / A[i][1]; - // Now obtain Ct2 in U.Ct2 = y ... - Ct2[count_cells + 1] = y[count_cells + 1]; - for (i = count_cells; i >= 0; i--) - Ct2[i] = y[i] - A[i][2] * Ct2[i + 1]; + A[i][0] = -mfr; A[i][1] = 1 + mfr + mfr1; A[i][2] = -mfr1; } + // decompose A in LU : store L in A[..][0..1] and U in A[..][2] ... + for (i = 1; i <= count_cells + 1; i++) + { + A[i - 1][2] = A[i - 1][2] / A[i - 1][1]; + A[i][1] = A[i][1] - A[i][0] * A[i - 1][2]; + } + // solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1 ... + // First, find y in L.y = Ct1 + y[0] = Ct2[0]; + for (i = 1; i <= count_cells + 1; i++) + y[i] = (Ct2[i] - A[i][0] * y[i - 1]) / A[i][1]; + // Now obtain Ct2 in U.Ct2 = y ... + Ct2[count_cells + 1] = y[count_cells + 1]; + for (i = count_cells; i >= 0; i--) + Ct2[i] = y[i] - A[i][2] * Ct2[i + 1]; + // Moles transported by concentration gradient from cell [i] to [i + 1] go in tot1. // Correct for electro-neutrality... for (i = ifirst; i <= ilast + 1; i++) @@ -2431,7 +2429,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) else ct[i].v_m[cp].zc = ct[i].v_m[cp].z * (Ct2[i] + Ct2[i + 1]) / 2; - current_cells[i].dif -= ct[i].v_m[cp].z * mixf[i][cp] * grad; + current_cells[i].dif += ct[i].v_m[cp].z * ct[i].J_ij[cp].tot1; current_cells[i].ele -= ct[i].v_m[cp].z * mixf[i][cp] * ct[i].v_m[cp].zc; } } @@ -2492,6 +2490,10 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) for (i = ifirst + 1; i < ilast; i++) { dVc = current_cells[i].R * (current_x - current_cells[i].dif); + if ((dV_dcell && dVc * j_0e > 0 || + (dV_dcell > 0 && (cell_data[i].potV + dVc) > (cell_data[count_cells + 1].potV)) || + (dV_dcell < 0 && (cell_data[i].potV + dVc) < (cell_data[count_cells + 1].potV)))) + dVc = (cell_data[count_cells + 1].potV - cell_data[i].potV) / (count_cells + 1 - i); cell_data[i + 1].potV = cell_data[i].potV + dVc; } if (!dV_dcell || fix_current) @@ -2516,7 +2518,6 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) } current_A = current_x / DDt * F_C_MOL; - for (i = ifirst; i <= ilast + 1; i++) { // preserve the potentials... @@ -2531,13 +2532,13 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) ct[i].m_s = (struct M_S *) free_check_null(ct[i].m_s); if (ct[i].m_s == NULL) { - ct[i].m_s = (struct M_S *) PHRQ_malloc((size_t)(count_m_s) * sizeof(struct M_S)); - ct[i].m_s_size = count_m_s; + ct[i].m_s = (struct M_S *) PHRQ_malloc((size_t)(count_m_s + 5) * sizeof(struct M_S)); + ct[i].m_s_size = count_m_s + 5; } else if (count_m_s > ct[i].m_s_size) { - ct[i].m_s = (struct M_S *) PHRQ_realloc(ct[i].m_s, (size_t)(count_m_s) * sizeof(struct M_S)); - ct[i].m_s_size = count_m_s; + ct[i].m_s = (struct M_S *) PHRQ_realloc(ct[i].m_s, (size_t)(count_m_s + 5) * sizeof(struct M_S)); + ct[i].m_s_size = count_m_s + 5; } if (ct[i].m_s == NULL) malloc_error(); @@ -2549,15 +2550,16 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) ct[i].m_s[i1].name = (*it).c_str(); it++; } - fill_m_s(ct[i].J_ij, ct[i].J_ij_count_spec, i); + fill_m_s(ct[i].J_ij, ct[i].J_ij_count_spec, i, 0); } /* * 3. find the solutions, add or subtract the moles... */ std::map> ::iterator it1; std::map ::iterator it2; - int icell, if1, il1, incr; - LDBLE dum1, dum2; + cxxNameDouble::iterator kit; + int icell, if1, il1, incr, length, length2; + LDBLE dum1, dum2, min_mol; for (cp = 0; cp < count_m_s; cp++) { if (!dV_dcell || dV_dcell * ct[1].m_s[cp].charge <= 0) @@ -2570,6 +2572,9 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) } for (icell = if1; icell != il1; icell += incr) { + min_mol = 1e-12 * ct[i].kgw; + if (fabs(ct[icell].m_s[cp].tot1) < min_mol) + continue; dum1 = dum2 = 0; sptr1 = Utilities::Rxn_find(Rxn_solution_map, icell); sptr2 = Utilities::Rxn_find(Rxn_solution_map, icell + 1); @@ -2577,7 +2582,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) { if (dV_dcell || (icell > 0 && icell <= ilast)) sptr1->Set_total_h(sptr1->Get_total_h() - ct[icell].m_s[cp].tot1); - if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last ==2)) + if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) sptr2->Set_total_h(sptr2->Get_total_h() + ct[icell].m_s[cp].tot1); continue; } @@ -2589,7 +2594,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) sptr2->Set_total_o(sptr2->Get_total_o() + ct[icell].m_s[cp].tot1); continue; } - // see if icell - incr has negative moles, subtract it, so that it is not transported... + // see if icell - incr has negative moles, subtract it, so that it is not transported... if (icell > 0 && icell <= ilast && (it1 = neg_moles.find(icell - incr)) != neg_moles.end() && (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end()) @@ -2601,19 +2606,129 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) neg_moles.insert(std::make_pair(icell - incr, els)); } dum1 = sptr1->Get_totals()[ct[icell].m_s[cp].name]; - dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name]; + if (!dum1) { - if (ct[icell].m_s[cp].tot1 > dum1) - ct[icell].m_s[cp].tot1 = dum1; - else if (-ct[icell].m_s[cp].tot1 > dum2) - ct[icell].m_s[cp].tot1 = -dum2; + /* Get moles from redox states... */ + length = (int)strlen(ct[icell].m_s[cp].name); + for (kit = sptr1->Get_totals().begin(); kit != sptr1->Get_totals().end(); kit++) + { + length2 = (int)(size_t)strcspn(kit->first.c_str(), "("); + if (!strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) + { + dum1 += kit->second; kit->second = 0; + } + } + if (dum1) + sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1; } + dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name]; + if (!dum2) + { + length = (int)strlen(ct[icell].m_s[cp].name); + for (kit = sptr2->Get_totals().begin(); kit != sptr2->Get_totals().end(); kit++) + { + length2 = (int)(size_t)strcspn(kit->first.c_str(), "("); + if (!strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) + { + dum2 += kit->second; kit->second = 0; + } + } + if (dum2) + sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2; + } + + // check for negative moles, add moles from the donnan layer when necessary and available... + if (ct[icell].m_s[cp].tot1 > dum1 && + (dV_dcell || (icell > 0 && icell <= ilast))) + { + if (!ct[icell].dl_s) + ct[icell].m_s[cp].tot1 = dum1; + else + { + cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell); + if (s_ptr) + { + cxxSurfaceCharge * charge_ptr = NULL; + for (size_t j = 0; j < s_ptr->Get_surface_charges().size(); j++) + { + if (s_ptr->Get_dl_type() == cxxSurface::DONNAN_DL) + { + charge_ptr = &(s_ptr->Get_surface_charges()[j]); + for (kit = charge_ptr->Get_diffuse_layer_totals().begin(); kit != charge_ptr->Get_diffuse_layer_totals().end(); kit++) + { + //if (strcmp(kit->first.c_str(), "H") == 0 || strcmp(kit->first.c_str(), "O") == 0) + // continue; + if (strcmp(kit->first.c_str(), ct[icell].m_s[cp].name) == 0) + { + dummy = ct[icell].m_s[cp].tot1 - dum1 + min_mol; + if (kit->second > dummy) + { + dum1 += dummy; kit->second -= dummy; + } + else + { + dum1 += kit->second; kit->second = 0; + } + } + } + } + } + } + } + sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1; + } + // check for negative moles, in the other cell... + ct[icell].m_s[cp].tot2 = ct[icell].m_s[cp].tot1; + if (-ct[icell].m_s[cp].tot2 > dum2 && + (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))) + { + if (!ct[icell + 1].dl_s) + ct[icell].m_s[cp].tot2 = -dum2; + else + { + cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell + 1); + if (s_ptr) + { + cxxSurfaceCharge * charge_ptr = NULL; + for (size_t j = 0; j < s_ptr->Get_surface_charges().size(); j++) + { + if (s_ptr->Get_dl_type() == cxxSurface::DONNAN_DL) + { + charge_ptr = &(s_ptr->Get_surface_charges()[j]); + for (kit = charge_ptr->Get_diffuse_layer_totals().begin(); kit != charge_ptr->Get_diffuse_layer_totals().end(); kit++) + { + //if (strcmp(kit->first.c_str(), "H") == 0 || strcmp(kit->first.c_str(), "O") == 0) + // continue; + if (strcmp(kit->first.c_str(), ct[icell].m_s[cp].name) == 0) + { + dummy = ct[icell].m_s[cp].tot2 + dum2 - min_mol; + if (kit->second > -dummy) + { + dum2 -= dummy; kit->second += dummy; + } + else + { + dum2 += kit->second; kit->second = 0; + } + } + } + } + } + } + sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2; + if (-ct[icell].m_s[cp].tot2 > dum2) + ct[icell].m_s[cp].tot2 = -dum2; + } + } + if (fabs(ct[icell].m_s[cp].tot2) < fabs(ct[icell].m_s[cp].tot1)) + ct[icell].m_s[cp].tot1 = ct[icell].m_s[cp].tot2; + if (dV_dcell || (icell > 0 && icell <= ilast)) { dum1 -= ct[icell].m_s[cp].tot1; - sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : 1e-16); + sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : 0e-16); } - if (dum1 < -1e-12 && incr > 0) + if (dum1 < -min_mol && incr > 0) { els.insert(std::make_pair(ct[icell].m_s[cp].name, dum1)); neg_moles.insert(std::make_pair(icell, els)); @@ -2622,9 +2737,9 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) if (dV_dcell || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)) { dum2 += ct[icell].m_s[cp].tot1; - sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : 1e-16); + sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : 0e-16); } - if (dum2 < -1e-12 && incr < 0) + if (dum2 < -min_mol && incr < 0) { els.insert(std::make_pair(ct[icell].m_s[cp].name, dum2)); neg_moles.insert(std::make_pair(icell + 1, els)); @@ -2808,7 +2923,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant) } count_m_s = 0; } - fill_m_s(ct[icell].J_ij, ct[icell].J_ij_count_spec, icell); + fill_m_s(ct[icell].J_ij, ct[icell].J_ij_count_spec, icell, stagnant); /* * 3. find the solutions, add or subtract the moles... @@ -2993,7 +3108,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant) /* ---------------------------------------------------------------------- */ int Phreeqc:: -fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell) +fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell, int stagnant) /* ---------------------------------------------------------------------- */ { /* sum up the primary or secondary master_species from solute species @@ -3005,6 +3120,8 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell) for (j = 0; j < l_J_ij_count_spec; j++) { + if (abs(l_J_ij[j].tot1) < 1e-15) + continue; { char * temp_name = string_duplicate(l_J_ij[j].name); ptr = temp_name; @@ -3012,7 +3129,7 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell) get_elts_in_species(&ptr, 1); free_check_null(temp_name); } - if (implicit) + if (implicit && !stagnant) { for (k = 0; k < count_elts; k++) { @@ -3020,14 +3137,14 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell) { if (strcmp(ct[icell].m_s[l].name, elt_list[k].elt->name) == 0) { - fraction = elt_list[k].coef * l_J_ij[j].tot1 + ct[icell].m_s[l].tot1; + fraction = fabs((double)elt_list[k].coef * l_J_ij[j].tot1) + fabs(ct[icell].m_s[l].tot1); if (fraction) - fraction = elt_list[k].coef * l_J_ij[j].tot1 / fraction; - else if (!ct[icell].m_s[l].tot1) + fraction = fabs((double)elt_list[k].coef * l_J_ij[j].tot1) / fraction; + else fraction = 1; ct[icell].m_s[l].tot1 += elt_list[k].coef * l_J_ij[j].tot1; - if (!strcmp(elt_list[k].elt->name, "H") || !strcmp(elt_list[k].elt->name, "O")) - break; + //if (!strcmp(elt_list[k].elt->name, "H") || !strcmp(elt_list[k].elt->name, "O")) + // break; ct[icell].m_s[l].charge *= (1 - fraction); ct[icell].m_s[l].charge += fraction * l_J_ij[j].charge; break; @@ -3057,21 +3174,16 @@ fill_m_s(struct J_ij *l_J_ij, int l_J_ij_count_spec, int icell) { if (strcmp(m_s[l].name, elt_list[k].elt->name) == 0) { - fraction = elt_list[k].coef * l_J_ij[j].tot1 / (elt_list[k].coef * l_J_ij[j].tot1 + m_s[l].tot1); m_s[l].tot1 += elt_list[k].coef * l_J_ij[j].tot1; m_s[l].tot2 += elt_list[k].coef * l_J_ij[j].tot2; - m_s[l].charge *= (1 - fraction); - m_s[l].charge += fraction * l_J_ij[j].charge; break; } } if (l == count_m_s) { - //m_s[l].name = string_hsave(elt_list[k].elt->name); m_s[l].name = elt_list[k].elt->name; m_s[l].tot1 = elt_list[k].coef * l_J_ij[j].tot1; m_s[l].tot2 = elt_list[k].coef * l_J_ij[j].tot2; - m_s[l].charge = l_J_ij[j].charge; count_m_s++; } } @@ -3131,7 +3243,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) LDBLE por_il1, por_il2, por_il12 = 0.0; LDBLE cec1, cec2, cec12 = 0.0, rc1 = 0.0, rc2 = 0.0; LDBLE dV, c1, c2; - LDBLE max_b = 0.0; // appt + LDBLE max_b = 0.0; il_calcs = (interlayer_Dflag ? 1 : 0); cxxSurface *s_ptr1, *s_ptr2; @@ -3486,9 +3598,6 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) while (i < i_max || j < j_max) { - // appt debug - //ct[icell].J_ij[i].name = sol_D[icell].spec[i].name; - //ct[icell].J_ij[i+1].name = sol_D[jcell].spec[j].name; if (j == j_max || (i < i_max && strcmp(sol_D[icell].spec[i].name, sol_D[jcell].spec[j].name) < 0)) { @@ -3748,7 +3857,7 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) if (ct[icell].v_m[k].z) ct[icell].v_m[k].zc = ct[icell].v_m[k].z * ct[icell].v_m[k].c; - if (dV_dcell && ct[icell].v_m[k].z && !fix_current && !implicit) // appt + if (dV_dcell && ct[icell].v_m[k].z && !fix_current && !implicit) { // compare diffuse and electromotive forces dum = ct[icell].v_m[k].grad; @@ -3811,15 +3920,25 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) else if (icell == count_cells) ct[icell].v_m[k].b_ij = b_i; } - if (ct[icell].v_m[k].b_ij > max_b) - max_b = ct[icell].v_m[k].b_ij; if (ct[icell].v_m[k].z) ct[icell].Dz2c += ct[icell].v_m[k].b_ij * ct[icell].v_m[k].zc * ct[icell].v_m[k].z; + //ddlm = sol_D[jcell].spec[j].lm - sol_D[icell].spec[i].lm; // appt: this could give an incorrect large factor for implicit + //if (fabs(ddlm) > 1e-10) + // ct[icell].v_m[k].grad *= (1 + (sol_D[jcell].spec[j].lg - sol_D[icell].spec[i].lg) / ddlm); ddlm = sol_D[jcell].spec[j].lm - sol_D[icell].spec[i].lm; - if (fabs(ddlm) > 1e-10) - ct[icell].v_m[k].grad *= (1 + (sol_D[jcell].spec[j].lg - sol_D[icell].spec[i].lg) / ddlm); + dum = sol_D[jcell].spec[j].lg - sol_D[icell].spec[i].lg; + if (fabs(dum * sol_D[jcell].spec[j].lg) > 5e-3 && fabs(ddlm) > 0.025) + { + dum = (1 + dum / ddlm); + ct[icell].v_m[k].grad *= dum; + if (implicit) + ct[icell].v_m[k].b_ij *= dum; + } + + if (ct[icell].v_m[k].b_ij > max_b) + max_b = ct[icell].v_m[k].b_ij; k++; } @@ -3829,8 +3948,8 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant) j++; } } - if (implicit) - return(max_b); // appt + if (implicit && !stagnant) + return(max_b); /* * fill in J_ij... */ @@ -3919,11 +4038,11 @@ dV_dcell2: } current_A = current_x * F_C_MOL; } - if (abs(ct[icell].J_ij_sum / DDt - current_x) > 1e-13 && (dV_dcell || (icell > 0 || jcell >= count_cells))) // appt - { - sprintf(token, "Moles charge added in cell %d is: %.3e.\n", icell, ct[icell].J_ij_sum / DDt - current_x); - screen_msg(token); - } + //if (fabs(ct[icell].J_ij_sum / DDt - current_x) > 1e-13 && (dV_dcell || (icell > 0 || jcell >= count_cells))) // appt + //{ + // sprintf(token, "Moles charge added in cell %d is: %.3e.\n", icell, ct[icell].J_ij_sum / DDt - current_x); + // screen_msg(token); + //} /* * calculate interlayer mass transfer... */ @@ -3945,7 +4064,7 @@ dV_dcell2: ct[icell].J_ij_il[i].tot1 *= ct[icell].A_ij_il * DDt; ct[icell].J_ij_sum += ct[icell].v_m_il[i].z * ct[icell].J_ij_il[i].tot1; ct[icell].J_ij_il[i].tot2 = ct[icell].J_ij_il[i].tot1; - ct[icell].J_ij_il[i].charge = ct[icell].v_m_il[i].z; + // ct[icell].J_ij_il[i].charge = ct[icell].v_m_il[i].z; // appt not used now } /* express the transfer in elemental moles... */ @@ -3963,7 +4082,7 @@ dV_dcell2: m_s[i1].tot2 = 0; } count_m_s = 0; - fill_m_s(ct[icell].J_ij_il, k_il, icell); + fill_m_s(ct[icell].J_ij_il, k_il, icell, stagnant); /* do the mass transfer... */ if (icell > 0 || stagnant) diff --git a/utilities.cpp b/utilities.cpp index a3d98b75..1a93255e 100644 --- a/utilities.cpp +++ b/utilities.cpp @@ -6,7 +6,6 @@ #include "Solution.h" #include - /* ---------------------------------------------------------------------- */ int Phreeqc:: add_elt_list(struct elt_list *elt_list_ptr, LDBLE coef) @@ -1516,7 +1515,7 @@ status(int count, const char *str, bool rk_string) } else { - screen_string = sformatf("%-15s%-27s%1s%37s", sim_str, state_str, spin_str, stdstr.c_str()); + screen_string = sformatf("%-15s%-27s%1s%45s", sim_str, state_str, spin_str, stdstr.c_str()); status_string = screen_string; } } From c92911348b9f5ebf3777e278fea7225f5900b9b2 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Sat, 1 Jun 2019 10:32:24 -0600 Subject: [PATCH 3/4] Tony fix May 31 --- transport.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/transport.cpp b/transport.cpp index 75413c97..d0f442ec 100644 --- a/transport.cpp +++ b/transport.cpp @@ -2613,7 +2613,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) for (kit = sptr1->Get_totals().begin(); kit != sptr1->Get_totals().end(); kit++) { length2 = (int)(size_t)strcspn(kit->first.c_str(), "("); - if (!strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) + if (length == length2 && !strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) { dum1 += kit->second; kit->second = 0; } @@ -2628,7 +2628,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) for (kit = sptr2->Get_totals().begin(); kit != sptr2->Get_totals().end(); kit++) { length2 = (int)(size_t)strcspn(kit->first.c_str(), "("); - if (!strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) + if (length == length2 && !strncmp(ct[icell].m_s[cp].name, kit->first.c_str(), length2)) { dum2 += kit->second; kit->second = 0; } From b3bf691b4964b3eedbaf04d6ee9bd35f98afb113 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Mon, 3 Jun 2019 21:05:15 -0600 Subject: [PATCH 4/4] fixed > > in templates for gcc --- Makefile.old | 4 ++-- transport.cpp | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/Makefile.old b/Makefile.old index b780ea89..c0b66adb 100644 --- a/Makefile.old +++ b/Makefile.old @@ -971,8 +971,8 @@ tester: # cd ../mytest; make clean; make -k -j 1 $(SPOOL) make.out $(SPOOL2); make diff $(SPOOL) diff.out $(SPOOL2) cd ../mytest $(CONCAT) make clean $(CONCAT) make -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make diff $(SPOOL) diff.out $(SPOOL2) cd ../examples $(CONCAT) make -f Makefile.old clean $(CONCAT) make -f Makefile.old -k -j 1 $(SPOOL) make.out $(SPOOL2) $(CONCAT) make -f Makefile.old diff $(SPOOL) diff.out $(SPOOL2) - svn status -q ../mytest - svn status -q ../examples +# svn status -q ../mytest +# svn status -q ../examples #ld-option # Usage: ldflags += $(call ld-option, -Wl$(comma)--hash-style=sysv) diff --git a/transport.cpp b/transport.cpp index d0f442ec..214e2f17 100644 --- a/transport.cpp +++ b/transport.cpp @@ -18,7 +18,7 @@ char token[MAX_LENGTH]; // implicit... std::set dif_spec_names; std::set dif_els_names; -std::map> neg_moles; +std::map > neg_moles; std::map els; double *y, *Ct1, *Ct2, *l_tk_x2, **A, **mixf, mixf_comp_size = 0; @@ -952,7 +952,7 @@ transport(void) if (implicit) { - std::map> ::iterator it1; + std::map > ::iterator it1; std::map ::iterator it2; for (i = 0; i <= all_cells; i++) { @@ -2555,7 +2555,7 @@ diffuse_implicit(LDBLE max_mixf, LDBLE DDt, int stagnant) /* * 3. find the solutions, add or subtract the moles... */ - std::map> ::iterator it1; + std::map > ::iterator it1; std::map ::iterator it2; cxxNameDouble::iterator kit; int icell, if1, il1, incr, length, length2;