max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-15 16:18:22 +01:00
Code Issues Packages Projects Releases Wiki Activity
iphreeqc/phreeqcpp/transport.cpp
Darth Vader 1b2eb6ca0d Squashed 'src/' changes from 6e248c34..81f180a0 
81f180a0 Fixed bugs in inverse
6d575967 strcpy_safe and strcat_safe
6d98c4e1 trying different header files
72796f15 added <cstring>
e8481607 warnings, strcat, strcpy
988bdee0 Try using goto(s)
9b10ce3f Try updated logical expression
812061be Turn off optimizing on k_temp

git-subtree-dir: src
git-subtree-split: 81f180a069285bcb4d180c860664b4b2a193854d
2023-11-16 05:04:39 +00:00

6191 lines
181 KiB
C++
Raw Blame History

#include "Utils.h"
#include "Phreeqc.h"
#include "phqalloc.h"
#include "Exchange.h"
#include "GasPhase.h"
#include "PPassemblage.h"
#include "SSassemblage.h"
#include "cxxKinetics.h"
#include "Solution.h"
#include <limits.h>
LDBLE F_Re3 = F_C_MOL / (R_KJ_DEG_MOL * 1e3);
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 <std::string> dif_spec_names;
std::set <std::string> dif_els_names;
std::map<int, std::map<std::string, double> > neg_moles;
std::map<std::string, double> els;
double *Ct2, *l_tk_x2, **A, **LU, **mixf, **mixf_stag;
int mixf_comp_size = 0;
struct CURRENT_CELLS
{
LDBLE dif, ele, R; // diffusive and electric components, relative cell resistance
} *current_cells;
LDBLE sum_R, sum_Rd; // sum of R, sum of (current_cells[0].dif - current_cells[i].dif) * R
struct V_M // For calculating Vinograd and McBain's zero-charge, diffusive tranfer of individual solutes
{
LDBLE grad, D, z, c, zc, Dz, Dzc;
LDBLE b_ij; // harmonic mean of cell properties, with EDL enrichment
};
struct CT /* summed parts of V_M and mcd transfer in a timestep for all cells, for free + DL water */
{
LDBLE kgw, dl_s, Dz2c, Dz2c_stag, 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;
class J_ij *J_ij, *J_ij_il;
int count_m_s;
class M_S *m_s;
int v_m_size, J_ij_size, m_s_size;
} *ct = NULL;
std::map<int, std::map<std::string, J_ij_save> > cell_J_ij;
struct MOLES_ADDED /* total moles added to balance negative conc's */
{
char *name;
LDBLE moles;
} *moles_added;
int count_moles_added;
#if defined(PHREEQCI_GUI)
#ifdef _DEBUG
#define new DEBUG_NEW
#undef THIS_FILE
static char THIS_FILE[] = __FILE__;
#endif
#endif
/* ---------------------------------------------------------------------- */
int Phreeqc::
transport(void)
/* ---------------------------------------------------------------------- */
{
int i, j, k, n;
int j_imm;
LDBLE b, f, mix_f_m, mix_f_imm;
LDBLE water_m, water_imm;
int first_c, last_c, b_c;
int max_iter;
char token[MAX_LENGTH];
LDBLE kin_time, stagkin_time, kin_time_save;
int punch_boolean = 0;
LDBLE step_fraction;
state = TRANSPORT;
diffc_tr = diffc;
diffc_max = 0.0;
transp_surf = warn_fixed_Surf = warn_MCD_X = 0;
dV_dcell = current_A = 0.0;
current_cells = NULL;
/* mass_water_switch = TRUE; */
/*
* Check existence of solutions
*/
j = -1;
/* check column solutions */
for (i = 1; i <= count_cells; i++)
{
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, i));
if (use.Get_solution_ptr() == NULL)
{
input_error++;
error_string = sformatf(
"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)
{
sol_D = (class sol_D *) PHRQ_malloc((size_t) (all_cells)* sizeof(class sol_D));
if (sol_D == NULL)
malloc_error();
for (i = 0; i < all_cells; i++)
{
sol_D[i].count_spec = 0;
sol_D[i].count_exch_spec = 0;
sol_D[i].exch_total = 0;
sol_D[i].x_max = 0;
sol_D[i].viscos_f = 1.0;
sol_D[i].tk_x = 298.15;
sol_D[i].spec = NULL;
sol_D[i].spec_size = 0;
}
//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 = ct[i].Dz2c_stag = 0.0;
ct[i].visc1 = 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;
}
count_moles_added = (int)elements.size();
moles_added = (struct MOLES_ADDED *) PHRQ_malloc((size_t) (count_moles_added)* sizeof(struct MOLES_ADDED));
if (moles_added == NULL)
malloc_error();
for (i = 0; i < count_moles_added; i++)
{
moles_added[i].name = NULL;
moles_added[i].moles = 0;
}
}
try
{
/* check solution 0 */
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, 0));
if (!use.Get_solution_ptr())
{
if (ishift == 1)
{
input_error++;
error_string = sformatf(
"Solution 0 is needed for transport, but is not defined.");
error_msg(error_string, CONTINUE);
}
else
Utilities::Rxn_copy(Rxn_solution_map, 1, 0);
}
else
{
if ((cell_data[0].potV = use.Get_solution_ptr()->Get_potV()))
dV_dcell = 1;
}
/* check solution count_cells */
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, count_cells + 1));
if (!use.Get_solution_ptr())
{
if (ishift == -1)
{
input_error++;
error_string = sformatf(
"Solution %d is needed for transport, but is not defined.",
count_cells + 1);
error_msg(error_string, CONTINUE);
}
else
Utilities::Rxn_copy(Rxn_solution_map, count_cells, count_cells + 1);
}
else
{
if ((cell_data[count_cells + 1].potV = use.Get_solution_ptr()->Get_potV()))
dV_dcell = 1;
}
/*
* Initialize temperature in stagnant cells ...
*/
for (n = 1; n <= stag_data.count_stag; n++)
{
for (i = 1; i <= count_cells; i++)
{
k = i + 1 + n * count_cells;
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, k));
if (use.Get_solution_ptr() != NULL)
cell_data[k].temp = use.Get_solution_ptr()->Get_tc();
if (n == 1 && implicit && use.Get_solution_ptr() == NULL)
{
input_error++;
error_string = sformatf(
"Stagnant solution %d not found for implicit diffusion with 1 stagnant layer.\n Please define it, or set -implicit false.", k);
error_msg(error_string, CONTINUE);
}
}
}
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;
dV_dcell = 1;
}
if (dV_dcell)
{
if (ishift)
{
input_error++;
error_string = sformatf(
"Electro-diffusion cannot be combined with advective transport.");
error_msg(error_string, CONTINUE);
}
if (!multi_Dflag)
{
input_error++;
error_string = sformatf(
"Electrical Field (potential) was defined, but needs -multi_D.");
error_msg(error_string, CONTINUE);
}
else
{
if (bcon_first != 1)
{
bcon_first = 1;
error_string = sformatf(
"Electrical Field (potential) was defined, assuming constant boundary condition for first cell.");
warning_msg(error_string);
}
if (bcon_last != 1)
{
bcon_last = 1;
error_string = sformatf(
"Electrical Field (potential) was defined, assuming constant boundary condition for last cell.");
warning_msg(error_string);
}
current_cells = (struct CURRENT_CELLS *) PHRQ_malloc(
(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;
}
}
}
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 && !timest)
{
error_string = sformatf(
"Time step not defined for diffusion, using -time_step 1 # seconds");
warning_msg(error_string);
timest = 1;
}
//if (implicit && stag_data.count_stag > 1)
//{
// error_string = sformatf(
// "Sorry, implicit diffusion can handle only 1 stagnant layer for now. Please remove -implicit true, or set -stagnant 1.");
// //error_msg(error_string, CONTINUE);
//}
if (implicit && current_cells == NULL)
{
current_cells = (struct CURRENT_CELLS *) PHRQ_malloc(
(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
*/
dup_print("Equilibrating initial solutions", TRUE);
transport_step = 0;
for (i = 0; i <= count_cells + 1; i++)
{
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);
if (use.Get_surface_ptr() != NULL && use.Get_surface_ptr()->Get_transport())
transp_surf = TRUE;
if (transp_surf && !multi_Dflag)
{
error_string = sformatf(
"-multi_d must be defined for surface transport");
error_msg(error_string, CONTINUE);
}
if (multi_Dflag == TRUE)
{
fill_spec(cell_no, 0);
}
print_punch(i, true);
/* if (i > 0 && i <= count_cells)*/
saver();
}
/*
* Also stagnant cells
*/
for (n = 1; n <= stag_data.count_stag; n++)
{
for (i = 1; i <= count_cells; i++)
{
k = i + 1 + n * count_cells;
cell_no = k;
if (Utilities::Rxn_find(Rxn_solution_map, k) != 0)
{
set_initial_moles(k);
set_and_run_wrapper(k, NOMIX, FALSE, k, 0.0);
if (multi_Dflag == TRUE)
{
fill_spec(cell_no, i);
}
print_punch(k, true);
saver();
}
}
}
/*
* Initialize mixing factors, define kinetics times
* for multicomponent diffusion, limit mixing by diffc_max (usually from H+)
*/
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)
{
snprintf(token, sizeof(token), "multi_D calc's and stagnant: define MIXing factors explicitly, or \n\t give in -multi_D the Dw used for calculating the mobile-immobile exchange factor.");
warning_msg(token);
}
Rxn_mix_map.clear();
}
/*
* mix[] is extended in init_mix(), to accommodate column mix factors
*/
nmix = init_mix();
heat_nmix = init_heat_mix(nmix);
if (nmix < 2)
stagkin_time = timest;
else
stagkin_time = timest / (double)nmix;
if (ishift != 0)
kin_time = timest / (1 + (double)nmix);
else
kin_time = stagkin_time;
kin_time_save = kin_time;
/* Reaction defined for a shift... */
if (!ishift)
{
if (nmix < 2)
step_fraction = 1.0;
else
step_fraction = 1.0 / nmix;
}
else
step_fraction = 1.0 / (1.0 + nmix);
/*
* Set boundary conditions, transport direction
*/
last_model.force_prep = true;
if ((ishift == 0) || (bcon_first == 1) || (bcon_last == 1))
b_c = 1;
else
b_c = 0;
if (ishift >= 0)
{
last_c = count_cells;
first_c = 1;
}
else
{
last_c = 1;
first_c = count_cells;
}
/*
* Define stagnant/mobile mix structure, if not read explicitly.
*
* With count_stag = 1, mix factors are calculated from exchange factor alpha
* (= exch_f), mobile th_m and immobile th_im porosity.
* These variables are read under keyword TRANSPORT, after stagnant, in
* structure stag_data.
* MIX 'cell_no' in input file can be an alternative for the calculation here.
*/
if ((stag_data.exch_f > 0) && (stag_data.count_stag == 1))
{
b = stag_data.th_m / (stag_data.th_m + stag_data.th_im);
f = exp(-stag_data.exch_f * stagkin_time / (b * stag_data.th_im));
mix_f_imm = b - b * f;
mix_f_m = mix_f_imm * stag_data.th_im / stag_data.th_m;
for (j = 1; j <= count_cells; j++)
{
j_imm = j + (1 + count_cells);
if (Utilities::Rxn_find(Rxn_solution_map, j) == NULL)
error_msg
("Could not find mobile cell solution in TRANSPORT.",
STOP);
if (Utilities::Rxn_find(Rxn_solution_map, j_imm) == NULL)
//error_msg
//("Could not find immobile cell solution in TRANSPORT.",
//STOP);
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();
/*
* Define C_m = (1 - mix_f_m) * C_m0 + mix_f_m) * C_im0
*/
{
cxxMix temp_mix;
temp_mix.Set_n_user(j);
temp_mix.Set_n_user_end(j);
temp_mix.Add(j, 1 - mix_f_m);
temp_mix.Add(j_imm, mix_f_m * water_m / water_imm);
Rxn_mix_map[j] = temp_mix;
}
/*
* Define C_im = mix_f_imm * C_m0 + (1 - mix_f_imm) * C_im0, or...
*/
{
cxxMix temp_mix;
temp_mix.Set_n_user(j_imm);
temp_mix.Set_n_user_end(j_imm);
temp_mix.Add(j_imm, 1 - mix_f_imm);
temp_mix.Add(j, mix_f_imm * water_imm / water_m);
Rxn_mix_map[j_imm] = temp_mix;
}
}
if (heat_nmix > 0)
{
/*
* Assumption: D_e used for calculating exch_f in input file equals diffc
*/
f = stag_data.exch_f * (heat_diffc - diffc) / diffc / tempr;
f = exp(-f * stagkin_time / (b * stag_data.th_im));
heat_mix_f_imm = b - b * f;
heat_mix_f_m =
heat_mix_f_imm * stag_data.th_im / stag_data.th_m;
}
}
/*
* Stop if error
*/
if (get_input_errors() > 0)
{
error_msg("Program terminating due to input errors.", STOP);
}
/*
* Now transport
*/
if (implicit)
snprintf(token, sizeof(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
snprintf(token, sizeof(token), "\nCalculating transport: %d (mobile) cells, %d shifts, %d mixruns...\n\n",
count_cells, count_shifts - transport_start + 1, nmix);
warning_msg(token);
max_iter = 0;
for (transport_step = transport_start; transport_step <= count_shifts;
transport_step++)
{
/*
* Set initial moles of phases
*/
for (i = 0; i <= count_cells + 1; i++)
{
if (!dV_dcell && !implicit && (i == 0 || i == count_cells + 1))
continue;
set_initial_moles(i);
}
/*
* Also stagnant cells
*/
for (n = 1; n <= stag_data.count_stag; n++)
{
for (i = 1; i <= count_cells; i++)
{
k = i + 1 + n * count_cells;
cell_no = k;
if (Utilities::Rxn_find(Rxn_solution_map, k) != 0)
set_initial_moles(k);
}
}
/*
* Start diffusing if boundary cond = 1, (fixed c, or closed)
*/
if (b_c == 1)
{
/* For half of mixing steps */
for (j = 1; j <= floor((LDBLE)nmix / 2); j++)
{
rate_sim_time_start = ((double)transport_step - 1) *
timest + ((double)j - 1) * kin_time;
rate_sim_time = rate_sim_time_start + kin_time;
mixrun = j;
if (multi_Dflag && j == floor((LDBLE)nmix / 2))
{
//snprintf(token, sizeof(token),
// "Transport step %3d. Multicomponent diffusion run %3d.",
// transport_step, j);
//dup_print(token, FALSE);
}
else if (!multi_Dflag)
{
snprintf(token, sizeof(token), "Transport step %3d. Mixrun %3d.",
transport_step, j);
dup_print(token, FALSE);
}
if (heat_nmix > 0 && !implicit)
{
heat_mix(heat_nmix);
/* equilibrate again ... */
for (i = 1; i <= count_cells; i++)
{
cell_no = i;
set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0);
if (multi_Dflag)
fill_spec(i, i - 1);
saver();
}
}
/* Go through cells */
if (transp_surf)
{
if (disp_surf(stagkin_time) == ERROR)
error_msg("Error in surface transport, stopping.",
STOP);
}
if (implicit)
diffuse_implicit(stagkin_time, stag_data.count_stag);
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) && !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;
if (multi_Dflag)
snprintf(token, sizeof(token),
"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
transport_step, j, i, max_iter);
else
snprintf(token, sizeof(token),
"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 (dV_dcell)
run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i
else
run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
}
else
{
run_reactions(i, kin_time, DISP, step_fraction); // n_saver = -2
}
if (multi_Dflag)
fill_spec(i, 0);
/* punch and output file */
if (ishift == 0 && j == nmix && (stag_data.count_stag == 0 || (implicit && stag_data.count_stag == 1)))
print_punch(i, true);
if (i > 1)
Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1);
saver();
}
if (!dV_dcell)
Utilities::Rxn_copy(Rxn_solution_map, -2, count_cells);
/* Stagnant zone mixing after completion of each
diffusive/dispersive step ... */
rate_sim_time_start = ((double)transport_step - 1) *
timest + ((double)j - 1) * stagkin_time;
rate_sim_time = rate_sim_time_start + stagkin_time;
if (stag_data.count_stag > 0)
{
if (ishift == 0 && j == nmix)
punch_boolean = TRUE;
else
punch_boolean = FALSE;
for (i = 0; i <= count_cells + 1; i++) // allow for stagnant cell mixing with boundary cells
{
mix_stag(i, stagkin_time, punch_boolean, step_fraction);
}
}
if (ishift == 0 && j == nmix && stag_data.count_stag > 0)
{
for (n = 1; n <= stag_data.count_stag; n++)
{
for (i = 1; i <= count_cells; i++)
{
k = i + 1 + n * count_cells;
if (Utilities::Rxn_find(Rxn_solution_map, k) == NULL)
continue;
print_punch(k, false);
}
}
}
}
}
/*
* Advective transport
*/
if (ishift != 0)
{
snprintf(token, sizeof(token), "Transport step %3d.", transport_step);
dup_print(token, FALSE);
if (b_c == 1)
rate_sim_time_start = ((double)transport_step - 1) *
timest + ((double)j - 1) * kin_time;
else
rate_sim_time_start = ((double)transport_step - 1) * timest;
rate_sim_time = rate_sim_time_start + kin_time;
/* halftime kinetics for resident water in first cell ... */
if (Utilities::Rxn_find(Rxn_kinetics_map, first_c) != NULL && count_cells > 1)
{
cell_no = first_c;
kin_time = kin_time_save / 2;
run_reactions(first_c, kin_time, NOMIX, 0.0);
saver();
kin_time = kin_time_save;
}
/* for each cell in column */
for (i = last_c; i != (first_c - ishift); i -= ishift)
Utilities::Rxn_copy(Rxn_solution_map, i - ishift, i);
/* if boundary_solutions must be flushed by the flux from the column...
if (ishift == 1 && bcon_last == 3)
solution_duplicate (last_c, last_c + 1);
else if (ishift == -1 && bcon_first == 3)
solution_duplicate (last_c, last_c - 1);
*/
if (transp_surf)
{
for (i = last_c + ishift; i != (first_c - ishift);
i -= ishift)
{
if ((ishift == 1 && i == last_c + 1 && bcon_last != 3) ||
(ishift == -1 && i == last_c - 1 && bcon_first != 3))
continue;
cxxSurface * surface_ptr = Utilities::Rxn_find(Rxn_surface_map, i - ishift);
if (surface_ptr == NULL)
{
if ((Utilities::Rxn_find(Rxn_surface_map, i) != NULL) &&
((i == 0 && bcon_first == 3)
|| (i == count_cells + 1 && bcon_last == 3)))
{
Rxn_surface_map.erase(i);
}
continue;
}
if (surface_ptr->Get_transport())
{
cxxSurface * surface_ptr1 = Utilities::Rxn_find(Rxn_surface_map, i);
if (surface_ptr1 == NULL)
{
cxxSurface surf;
surf.Set_n_user(i);
surf.Set_n_user_end(i);
Rxn_surface_map[i] = surf;
}
if (i == first_c)
{
Rxn_surface_map[i] = mobile_surface_copy(surface_ptr, i, false);
}
else
{
Rxn_surface_map[i] = mobile_surface_copy(surface_ptr, i, true);
}
}
}
}
/*
* thermal diffusion when nmix = 0...
*/
if (nmix == 0 && heat_nmix > 0 && !implicit)
{
heat_mix(heat_nmix);
/* equilibrate again ... */
for (i = 1; i <= count_cells; i++)
{
cell_no = i;
set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0);
if (multi_Dflag)
fill_spec(i, i - 1);
saver();
}
}
for (i = 1; i <= count_cells; i++)
{
if (i == first_c && count_cells > 1)
kin_time /= 2;
cell_no = i;
if (multi_Dflag)
snprintf(token, sizeof(token),
"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
transport_step, 0, i, max_iter);
else
snprintf(token, sizeof(token),
"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)
fill_spec(i, i - 1);
if (overall_iterations > max_iter)
max_iter = overall_iterations;
if (nmix == 0 && stag_data.count_stag == 0)
print_punch(i, true);
if (i == first_c && count_cells > 1)
kin_time = kin_time_save;
saver();
/* If nmix is zero, stagnant zone mixing after advective step ... */
if ((nmix == 0) && (stag_data.count_stag > 0))
{
mix_stag(i, stagkin_time, TRUE, step_fraction);
}
}
if (nmix == 0 && stag_data.count_stag > 0)
{
for (n = 1; n <= stag_data.count_stag; n++)
{
for (i = 1; i <= count_cells; i++)
{
k = i + 1 + n * count_cells;
if (Utilities::Rxn_find(Rxn_solution_map, k) == NULL)
continue;
print_punch(k, false);
}
}
}
}
/*
* Further dispersive and diffusive transport
*/
if (b_c != 1)
j = 1;
for (; j <= nmix; j++) // loop on j
{
mixrun = j;
if (multi_Dflag && j == nmix && (transport_step % print_modulus == 0))
{
snprintf(token, sizeof(token),
"Transport step %3d. Multicomponent diffusion run %3d.",
transport_step, j);
dup_print(token, FALSE);
}
else if (!multi_Dflag)
{
snprintf(token, sizeof(token), "Transport step %3d. Mixrun %3d.",
transport_step, j);
dup_print(token, FALSE);
}
rate_sim_time_start = ((double)transport_step - 1) *
timest + ((double)j - 1) * kin_time;
if (ishift != 0)
rate_sim_time_start += kin_time;
rate_sim_time = rate_sim_time_start + kin_time;
if (heat_nmix > 0 && !implicit)
{
heat_mix(heat_nmix);
/* equilibrate again ... */
for (i = 1; i <= count_cells; i++)
{
cell_no = i;
set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0);
if (multi_Dflag)
fill_spec(i, i - 1);
saver();
}
}
if (transp_surf)
{
if (disp_surf(stagkin_time) == ERROR)
error_msg("Error in surface transport, stopping.", STOP);
}
if (implicit)
diffuse_implicit(stagkin_time, stag_data.count_stag);
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) && !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;
if (multi_Dflag)
snprintf(token, sizeof(token),
"Transport step %3d. MCDrun %3d. Cell %3d. (Max. iter %3d)",
transport_step, j, i, max_iter);
else
snprintf(token, sizeof(token),
"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 (dV_dcell)
run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i
else
run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i
}
else
{
run_reactions(i, kin_time, DISP, step_fraction);
}
if (multi_Dflag == TRUE)
fill_spec(i, 0);
if (j == nmix && (stag_data.count_stag == 0 || (implicit && stag_data.count_stag == 1)))
print_punch(i, true);
if (i > 1)
Utilities::Rxn_copy(Rxn_solution_map, -2, i - 1);
saver();
}
if (!dV_dcell)
Utilities::Rxn_copy(Rxn_solution_map, -2, count_cells);
/* Stagnant zone mixing after completion of each diffusive/dispersive step ... */
rate_sim_time_start = ((double)transport_step - 1) *
timest + ((double)j - 1) * stagkin_time;
rate_sim_time = rate_sim_time_start + stagkin_time;
if (stag_data.count_stag > 0)
{
if (j == nmix)
punch_boolean = TRUE;
else
punch_boolean = FALSE;
for (i = 0; i <= count_cells + 1; i++) // allow for stagnant cell mixing with boundary cells
{
mix_stag(i, stagkin_time, punch_boolean, step_fraction);
}
}
if (j == nmix && ((stag_data.count_stag > 0/* && !implicit) || (implicit && stag_data.count_stag == 1*/)))
{
for (n = 1; n <= stag_data.count_stag; n++)
{
for (i = 1; i <= count_cells; i++)
{
k = i + 1 + n * count_cells;
if (Utilities::Rxn_find(Rxn_solution_map, k) == NULL)
continue;
cell_no = k;
print_punch(k, false);
}
}
}
}
if (dump_modulus != 0 && (transport_step % dump_modulus) == 0)
dump();
}
screen_msg("\n");
if (implicit && neg_moles.size())
{
std::map<int, std::map<std::string, double> > ::iterator it1;
std::map<std::string, double> ::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_moles_added; 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)
{
snprintf(token, sizeof(token),
"\nFor balancing negative concentrations in MCD, added in total to the system:");
warning_msg(token);
for (i = 0; i < count_moles_added; i++)
{
if (!moles_added[i].moles)
break;
snprintf(token, sizeof(token),
"\t %.4e moles %s.",
(double)moles_added[i].moles, moles_added[i].name);
warning_msg(token);
}
}
}
catch (...)
{
transport_cleanup();
throw;
}
transport_cleanup();
initial_total_time += rate_sim_time;
rate_sim_time = 0;
mass_water_switch = FALSE;
return (OK);
}
/* ---------------------------------------------------------------------- */
void Phreeqc::
transport_cleanup(void)
/* ---------------------------------------------------------------------- */
{
int i;
/*
* free mix structures
*/
Dispersion_mix_map.clear();
if ((stag_data.exch_f > 0) && (stag_data.count_stag == 1))
{
Rxn_mix_map.clear();
}
if (heat_nmix > 0)
{
heat_mix_array = (LDBLE *)free_check_null(heat_mix_array);
temp1 = (LDBLE *)free_check_null(temp1);
temp2 = (LDBLE *)free_check_null(temp2);
}
if (multi_Dflag)
{
for (i = 0; i < all_cells; i++)
{
sol_D[i].spec = (class spec *) free_check_null(sol_D[i].spec);
}
sol_D = (class sol_D *) free_check_null(sol_D);
for (int i = 0; i < all_cells; i++)
{
ct[i].v_m = (struct V_M *) free_check_null(ct[i].v_m);
ct[i].v_m_il = (struct V_M *) free_check_null(ct[i].v_m_il);
ct[i].J_ij = (class J_ij *) free_check_null(ct[i].J_ij);
ct[i].J_ij_il = (class J_ij *) free_check_null(ct[i].J_ij_il);
ct[i].m_s = (class M_S *) free_check_null(ct[i].m_s);
}
ct = (struct CT *) free_check_null(ct);
for (int i = 0; i < count_moles_added; i++)
{
moles_added[i].name = (char *)free_check_null(moles_added[i].name);
}
moles_added = (struct MOLES_ADDED *) free_check_null(moles_added);
}
if (implicit)
{
int l_stag = (stag_data.count_stag < 2 ? stag_data.count_stag : 0);
Ct2 = (LDBLE *)free_check_null(Ct2);
l_tk_x2 = (LDBLE *)free_check_null(l_tk_x2);
if (A)
{
for (i = 0; i < count_cells + 2 + l_stag * count_cells; i++)
{
A[i] = (LDBLE *)free_check_null(A[i]);
LU[i] = (LDBLE *)free_check_null(LU[i]);
}
}
if (mixf)
{
for (i = 0; i < count_cells + 2; i++)
{
mixf[i] = (LDBLE *)free_check_null(mixf[i]);
if (l_stag)
mixf_stag[i] = (LDBLE *)free_check_null(mixf_stag[i]);
if (!dV_dcell && !fix_current)
{
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);
LU = (LDBLE **)free_check_null(LU);
mixf = (LDBLE **)free_check_null(mixf);
mixf_stag = (LDBLE **)free_check_null(mixf_stag);
dif_spec_names.clear();
mixf_comp_size = 0;
}
current_cells = (struct CURRENT_CELLS *) free_check_null(current_cells);
}
/* ---------------------------------------------------------------------- */
void Phreeqc::
print_punch(int i, boolean active)
/* ---------------------------------------------------------------------- */
{
if ((!(cell_data[i].punch && (transport_step % punch_modulus == 0)) &&
!(cell_data[i].print && (transport_step % print_modulus == 0))) ||
(bcon_first == 2 && i == 0) ||
(bcon_last == 2 && i == count_cells + 1))
return;
if (!active)
run_reactions(i, 0, NOMIX, 0);
cell_no = i;
if (dV_dcell || fix_current)
{
use.Set_n_solution_user(i);
use.Get_solution_ptr()->Set_potV(cell_data[i].potV);
potV_x = cell_data[i].potV;
}
use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, i));
if (use.Get_kinetics_ptr() != NULL)
{
use.Set_n_kinetics_user(i);
use.Set_kinetics_in(true);
}
if (cell_data[i].punch && (transport_step % punch_modulus == 0))
punch_all();
if (cell_data[i].print && (transport_step % print_modulus == 0))
print_all();
/* maybe sorb a surface component... */
if (change_surf_count > 0)
{
for (int k = 0; k < change_surf_count; k++)
{
if (change_surf[k].cell_no != i)
break;
reformat_surf(change_surf[k].comp_name,
change_surf[k].fraction,
change_surf[k].new_comp_name,
change_surf[k].new_Dw,
change_surf[k].cell_no);
change_surf[k].cell_no = -99;
}
change_surf_count = 0;
save.n_surface_user = save.n_solution_user;
save.n_surface_user_end = save.n_solution_user_end;
}
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
init_mix(void)
/* ---------------------------------------------------------------------- */
{
LDBLE lav, dav = 0.0, mf12, maxmix, corr_disp, diffc_here, mD;
bool warning = false;
int i, l_nmix;
LDBLE *m, *m1;
m = (LDBLE *)PHRQ_malloc((count_cells + 1) * sizeof(LDBLE));
if (m == NULL)
malloc_error();
m1 = (LDBLE *)PHRQ_malloc((count_cells + 1) * sizeof(LDBLE));
if (m1 == NULL)
malloc_error();
for (i = 0; i < count_cells + 1; i++)
{
m[i] = 0;
m1[i] = 0;
}
corr_disp = 1.;
if (correct_disp == TRUE && ishift != 0)
{
if (bcon_first == 3)
corr_disp += 1. / count_cells;
if (bcon_last == 3)
corr_disp += 1. / count_cells;
}
maxmix = 0.0;
if (multi_Dflag)
{
if (dV_dcell)
dV_dcell = (cell_data[count_cells + 1].potV - cell_data[0].potV) / count_cells;
for (i = 1; i <= count_cells; i++)
{
// estimate the number of mixes
if (i < count_cells)
{
lav = (cell_data[i + 1].length + cell_data[i].length) / 2;
mD = diffc_max * timest / (lav * lav);
if (mD > maxmix)
maxmix = mD;
}
if (ishift != 0)
{
/* define dispersive mixing (diffusion is done in multi_D),
with harmonic mean of cell properties (manual eqn 49):
M1 = -a1 * (v1 * timest / 2) * A1 * (c12 - c1) / (h1 / 2)
a1 is dispersivity in cell 1, A1 = V1 / h1, v1 * timest / 2 = h1 / 2.
All cells have equal V in advective transport, delta(c) = M1 / V1.
*/
if (i < count_cells)
{
if (cell_data[i].disp)
dav = cell_data[i].length / cell_data[i].disp;
if (cell_data[i + 1].disp)
dav += cell_data[i + 1].length / cell_data[i + 1].disp;
if (dav)
m1[i] = 2 * corr_disp / dav; /* m1[i] has mixf with higher cell */
}
if (i > 1)
{
if (cell_data[i].disp)
dav = cell_data[i].length / cell_data[i].disp;
if (cell_data[i - 1].disp)
dav += cell_data[i - 1].length / cell_data[i - 1].disp;
if (dav)
m[i] = 2 * corr_disp / dav; /* m[i] has mixf with lower cell */
}
mf12 = m[i] + m1[i];
if (mf12 > maxmix)
maxmix = mf12;
}
}
/*
* Also for boundary cells
*/
if (bcon_first == 1)
{
mD = 2 * diffc_max * timest / (cell_data[1].length * cell_data[1].length);
if (mD > maxmix)
maxmix = mD;
if (ishift != 0)
{
m[1] = 2 * cell_data[1].disp / cell_data[1].length * corr_disp;
mf12 = m[1] + m1[1];
if (mf12 > maxmix)
maxmix = mf12;
}
}
if (bcon_last == 1)
{
mD = 2 * diffc_max * timest / (cell_data[count_cells].length * cell_data[count_cells].length);
if (mD > maxmix)
maxmix = mD;
if (ishift != 0)
{
m1[count_cells] = 2 * cell_data[count_cells].disp / cell_data[count_cells].length * corr_disp;
mf12 = m[count_cells] + m1[count_cells];
if (mf12 > maxmix)
maxmix = mf12;
}
}
/*
* Find number of mixes
*/
if (maxmix == 0)
{
l_nmix = 0;
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))
{
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)
{
m = (LDBLE *)free_check_null(m);
m1 = (LDBLE *)free_check_null(m1);
char token[MAX_LENGTH];
snprintf(token, sizeof(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)
l_nmix = 1 + (int)floor(2.25 * maxmix);
else
l_nmix = 1 + (int)floor(1.5 * maxmix);
if (ishift != 0 && (bcon_first == 1 || bcon_last == 1))
{
if (l_nmix < 2)
l_nmix = 2;
}
if (mcd_substeps > 1)
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);
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);
return l_nmix;
}
else // multi_D false
{
diffc_here = 2 * diffc_tr * timest;
/*
* Define mixing factors among inner cells
*/
for (i = 1; i <= count_cells; i++)
{
if (i < count_cells)
{
// find mix with higher numbered cell...
if (ishift != 0)
{
if (cell_data[i].disp)
dav = cell_data[i].length / cell_data[i].disp;
if (cell_data[i + 1].disp)
dav += cell_data[i + 1].length / cell_data[i + 1].disp;
if (dav)
m1[i] = 2 / dav;
}
// add diffusive mixf...
m1[i] += diffc_here /
(cell_data[i].length * cell_data[i].length + cell_data[i].length * cell_data[i + 1].length);
m1[i] *= corr_disp; /* m1[i] has mixf with higher cell */
}
if (i > 1)
{
// and with lower numbered cell...
if (ishift != 0)
{
if (cell_data[i].disp)
dav = cell_data[i].length / cell_data[i].disp;
if (cell_data[i - 1].disp)
dav += cell_data[i - 1].length / cell_data[i - 1].disp;
if (dav)
m[i] = 2 / dav;
}
// add diffusive mixf...
m[i] += diffc_here /
(cell_data[i].length * cell_data[i].length + cell_data[i].length * cell_data[i - 1].length);
m[i] *= corr_disp; /* m[i] has mixf with lower numbered cell */
if (m[i] != m1[i - 1] && !warning && (!dav || m[i] / (2 / dav) > 1.00001))
{
warning_msg("Unequal cell-lengths may give mass-balance error, consider using -multi_D");
warning = true;
}
}
mf12 = m[i] + m1[i];
if (mf12 > maxmix)
maxmix = mf12;
}
/*
* Also for boundary cells
*/
if (bcon_first == 1)
{
m[1] = diffc_here / (cell_data[1].length * cell_data[1].length);
if (ishift != 0)
m[1] += cell_data[1].disp / cell_data[1].length;
mf12 = m[1] + m1[1];
if (mf12 > maxmix)
maxmix = mf12;
}
if (bcon_last == 1)
{
m1[count_cells] = diffc_here / (cell_data[count_cells].length * cell_data[count_cells].length);
if (ishift != 0)
m1[count_cells] += cell_data[count_cells].disp / cell_data[count_cells].length;
mf12 = m[count_cells] + m1[count_cells];
if (mf12 > maxmix)
maxmix = mf12;
}
/*
* Find number of mixes
*/
if (maxmix == 0)
l_nmix = 0;
else
{
if (1.5 * maxmix > (double)INT_MAX)
{
m = (LDBLE *)free_check_null(m);
m1 = (LDBLE *)free_check_null(m1);
char token[MAX_LENGTH];
snprintf(token, sizeof(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);
if ((ishift != 0) && ((bcon_first == 1) || (bcon_last == 1)))
{
if (l_nmix < 2)
l_nmix = 2;
}
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;
}
}
m = (LDBLE *)free_check_null(m);
m1 = (LDBLE *)free_check_null(m1);
return (l_nmix);
}
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction)
/* ---------------------------------------------------------------------- */
{
int n, k;
LDBLE t_imm;
cxxSolution *ptr_imm, *ptr_m;
k = -1000; // compiler says k may be undefined
ptr_imm = NULL;
boolean done_mixing = false;
/*
* Kinetics in transport cell is done while transporting
*/
for (n = 1; n <= stag_data.count_stag; n++)
{
if (i == 0 || i == count_cells + 1)
{
use.Set_mix_ptr(NULL);
use.Set_mix_in(false);
use.Set_mix_ptr(Utilities::Rxn_find(Rxn_mix_map, i));
if (use.Get_mix_ptr())
{
for (std::map < int, LDBLE >::const_iterator it = use.Get_mix_ptr()->Get_mixComps().begin();
it != use.Get_mix_ptr()->Get_mixComps().end(); it++)
{
if (it->first > i && it->first < all_cells && it->first != count_cells + 1)
{
k = it->first;
ptr_imm = Utilities::Rxn_find(Rxn_solution_map, k);
break;
}
}
}
}
else
{
k = i + 1 + n * count_cells;
if (k < all_cells)
ptr_imm = Utilities::Rxn_find(Rxn_solution_map, k);
}
if (ptr_imm != NULL)
{
if (n == 1)
{
if (heat_nmix > 0 && (!implicit || (implicit && stag_data.count_stag > 1)))
{
ptr_m = Utilities::Rxn_find(Rxn_solution_map, i);
t_imm =
heat_mix_f_imm * ptr_m->Get_tc() + (1 - heat_mix_f_imm) * ptr_imm->Get_tc();
ptr_m->Set_tc(heat_mix_f_m * ptr_imm->Get_tc() + (1 - heat_mix_f_m) * ptr_m->Get_tc());
cell_data[i].temp = ptr_m->Get_tc();
cell_data[k].temp = t_imm = ptr_imm->Get_tc();
/* equilibrate again ... */
cell_no = i;
set_and_run_wrapper(i, NOMIX, FALSE, i, 0.0);
if (multi_Dflag == TRUE)
fill_spec(cell_no, 0);
saver();
cell_no = k;
set_and_run_wrapper(k, NOMIX, FALSE, k, 0.0);
if (multi_Dflag == TRUE)
fill_spec(cell_no, i);
saver();
}
/*
* Mobile cell, kinetics already done ...
*/
cell_no = i;
if (transp_surf)
{
if (diff_stag_surf(i) == ERROR)
error_msg("Error in surface transport, stopping.",
STOP);
}
if (!implicit || (implicit && stag_data.count_stag > 1))
{
if (multi_Dflag == TRUE)
multi_D(1.0, i, 2);
set_and_run_wrapper(i, STAG, FALSE, -2, 0.0);
if (multi_Dflag == TRUE)
fill_spec(cell_no, 0);
if (l_punch)
print_punch(i, true);
saver(); // save solution i in -2, original can be used in other stagnant mixes
}
}
cell_no = k;
if (implicit)
run_reactions(k, kin_time, NOMIX, step_fraction);
else
run_reactions(k, kin_time, STAG, step_fraction);
if (multi_Dflag == TRUE)
fill_spec(cell_no, i);
saver(); // save solution k in -2 - k, original k can be used in other stagnant mixes
done_mixing = true;
}
else if (n == 1 && l_punch && !implicit)
print_punch(i, false);
}
if (done_mixing) // after all mixing is done, the temporal solution becomes the original for the next timestep
{
for (n = 1; n <= stag_data.count_stag; n++)
{
k = i + 1 + n * count_cells;
if (Utilities::Rxn_find(Rxn_solution_map, k) != 0)
{
Utilities::Rxn_copy(Rxn_solution_map, -2 - k, k);
if (n == 1 && !implicit)
Utilities::Rxn_copy(Rxn_solution_map, -2, i);
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
init_heat_mix(int l_nmix)
/* ---------------------------------------------------------------------- */
{
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 && !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++)
{
if (fabs(cell_data[i].temp - t0) > 1.0)
{
l_heat_nmix = 1;
break;
}
}
if (l_heat_nmix == 0)
{
if (fabs(Utilities::Rxn_find(Rxn_solution_map, count_cells + 1)->Get_tc() - t0) > 1.0)
l_heat_nmix = 1;
for (n = 1; n <= stag_data.count_stag; n++)
{
for (i = 1; i < count_cells; i++)
{
k = i + 1 + n * count_cells;
if (Utilities::Rxn_find(Rxn_solution_map, k) != 0)
{
if (fabs(cell_data[k].temp - t0) > 1.0)
{
l_heat_nmix = 1;
break;
}
}
}
}
}
if (l_heat_nmix == 0)
return (0);
/*
* Initialize arrays...
*/
heat_mix_array = (LDBLE *)PHRQ_malloc((count_cells + 2) * sizeof(LDBLE));
if (heat_mix_array == NULL)
malloc_error();
temp1 = (LDBLE *)PHRQ_malloc((count_cells + 2) * sizeof(LDBLE));
if (temp1 == NULL)
malloc_error();
temp2 = (LDBLE *)PHRQ_malloc((count_cells + 2) * sizeof(LDBLE));
if (temp2 == NULL)
malloc_error();
/*
* Define mixing factors among inner cells...
*/
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 / mixf;
if (bcon_last == 3)
corr_disp += 1. / count_cells / mixf;
}
maxmix = 0.0;
for (i = 1; i < count_cells; i++)
{
lav = (cell_data[i + 1].length + cell_data[i].length) / 2;
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 */
}
/*
* Also for boundary cells
*/
if (bcon_first == 1)
{
lav = cell_data[1].length;
mixf = (l_diffc) * timest * corr_disp / tempr / (lav * lav);
if (2 * mixf > maxmix)
maxmix = 2 * mixf;
heat_mix_array[1] = 2 * mixf;
}
else
heat_mix_array[1] = 0;
if (bcon_last == 1)
{
lav = cell_data[count_cells].length;
mixf = (l_diffc) * timest * corr_disp / tempr / (lav * lav);
if (2 * mixf > maxmix)
maxmix = 2 * mixf;
heat_mix_array[count_cells + 1] = 2 * mixf;
}
else
heat_mix_array[count_cells + 1] = 0;
/*
* Find number of mixes
*/
if (maxmix == 0)
l_heat_nmix = 0;
else
{
if (implicit)
{
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);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
heat_mix(int l_heat_nmix)
/* ---------------------------------------------------------------------- */
{
int i, j;
for (i = 1; i <= count_cells; i++)
temp1[i] = Utilities::Rxn_find(Rxn_solution_map, i)->Get_tc();
temp1[0] = Utilities::Rxn_find(Rxn_solution_map, 0)->Get_tc();
temp1[count_cells + 1] =
Utilities::Rxn_find(Rxn_solution_map, (count_cells + 1))->Get_tc();
for (i = 1; i <= l_heat_nmix; i++)
{
for (j = 1; j <= count_cells; j++)
temp2[j] =
heat_mix_array[j] * temp1[j - 1] + heat_mix_array[j + 1] * temp1[j + 1] +
(1 - heat_mix_array[j] - heat_mix_array[j + 1]) * temp1[j];
for (j = 1; j <= count_cells; j++)
temp1[j] = temp2[j];
}
for (i = 1; i <= count_cells; i++)
{
cell_data[i].temp = temp1[i];
Utilities::Rxn_find(Rxn_solution_map, i)->Set_tc(temp1[i]);
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
set_initial_moles(int i)
/* ---------------------------------------------------------------------- */
{
cxxKinetics *kinetics_ptr;
char token[MAX_LENGTH];
const char* cptr;
int j, k, l;
/*
* Pure phase assemblage
*/
{
cxxPPassemblage * pp_assemblage_ptr = Utilities::Rxn_find(Rxn_pp_assemblage_map, i);
if (pp_assemblage_ptr != NULL)
{
std::map<std::string, cxxPPassemblageComp>::iterator it;
it = pp_assemblage_ptr->Get_pp_assemblage_comps().begin();
for (; it != pp_assemblage_ptr->Get_pp_assemblage_comps().end(); it++)
{
it->second.Set_initial_moles(it->second.Get_moles());
if (it->second.Get_initial_moles() < 0)
it->second.Set_initial_moles(0.0);
}
}
}
/*
* Gas phase
*/
{
cxxGasPhase * gas_phase_ptr = Utilities::Rxn_find(Rxn_gas_phase_map, i);
if (gas_phase_ptr != NULL)
{
std::vector<cxxGasComp> gc = gas_phase_ptr->Get_gas_comps();
for (size_t l = 0; l < gc.size(); l++)
{
gc[l].Set_initial_moles(gc[l].Get_moles());
}
gas_phase_ptr->Set_gas_comps(gc);
}
}
/*
* Kinetics
*/
kinetics_ptr = Utilities::Rxn_find(Rxn_kinetics_map, i);
if (kinetics_ptr != NULL)
{
for (j = 0; j < (int) kinetics_ptr->Get_kinetics_comps().size(); j++)
{
cxxKineticsComp *kinetics_comp_ptr = &(kinetics_ptr->Get_kinetics_comps()[j]);
kinetics_comp_ptr->Set_initial_moles(kinetics_comp_ptr->Get_m());
}
}
/*
* Solid solutions
*/
{
cxxSSassemblage *ss_assemblage_ptr = Utilities::Rxn_find(Rxn_ss_assemblage_map, i);
if (ss_assemblage_ptr != NULL)
{
std::vector<cxxSS *> ss_ptrs = ss_assemblage_ptr->Vectorize();
for (k = 0; k < (int) ss_ptrs.size(); k++)
{
cxxSS * ss_ptr = ss_ptrs[k];
for (j = 0; j < (int) ss_ptr->Get_ss_comps().size(); j++)
{
cxxSScomp *comp_ptr = &(ss_ptr->Get_ss_comps()[j]);
comp_ptr->Set_init_moles(comp_ptr->Get_moles());
}
}
}
}
/*
* For interlayer diffusion: add tiny bit of exchanger if absent
*/
cxxExchange * exchange_ptr = Utilities::Rxn_find(Rxn_exchange_map, i);
if (interlayer_Dflag && exchange_ptr == NULL)
{
cxxExchange temp_exchange;
temp_exchange.Set_n_user_both(i);
temp_exchange.Set_description("Interlayer diffusion: added 2e-10 moles X-");
use.Set_exchange_in(true);
use.Set_n_exchange_user(i);
temp_exchange.Set_new_def(true);
temp_exchange.Set_solution_equilibria(true);
temp_exchange.Set_n_solution(i);
cxxExchComp comp;
count_elts = 0;
paren_count = 0;
Utilities::strcpy_safe(token, MAX_LENGTH, "X");
cptr = token;
get_elts_in_species(&cptr, 2e-10);
cptr = token;
LDBLE z;
{
std::string token1;
get_token(&cptr, token1, &z, &l);
comp.Set_formula(token1.c_str());
}
comp.Set_formula_z(z);
comp.Set_totals(elt_list_NameDouble());
comp.Set_charge_balance(0.0);
temp_exchange.Get_exchange_comps().push_back(comp);
Rxn_exchange_map[i] = temp_exchange;
state = INITIAL_EXCHANGE;
initial_exchangers(TRUE);
state = TRANSPORT;
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
fill_spec(int l_cell_no, int ref_cell)
/* ---------------------------------------------------------------------- */
{
/* copy species activities into sol_D.spec... */
int i, i1, i2, i3, count_spec, count_exch_spec, size_xt;
char token[MAX_LENGTH];
const char * name;
class species *s_ptr, *s_ptr2;
class master *master_ptr;
LDBLE dum, dum2;
LDBLE lm;
LDBLE por, por_il, viscos_f, viscos_il_f, viscos;
bool x_max_done = false;
std::set <std::string> loc_spec_names;
s_ptr2 = NULL;
// for implicit
std::set <std::string> ::iterator it;
std::pair <std::set<std::string> ::iterator, bool> name_ret;
if (implicit && !l_cell_no)
dif_spec_names.clear();
size_xt = 5;
//sol_D[l_cell_no].spec = (class spec *) free_check_null(sol_D[l_cell_no].spec);
if (sol_D[l_cell_no].spec == NULL)
{
sol_D[l_cell_no].spec = (class spec *) PHRQ_malloc((species_list.size() + (size_t)size_xt) * sizeof(class spec));
sol_D[l_cell_no].spec_size = (int)species_list.size() + size_xt;
}
else if ((int)species_list.size() + size_xt > sol_D[l_cell_no].spec_size)
{
sol_D[l_cell_no].spec = (class spec *) PHRQ_realloc(sol_D[l_cell_no].spec, (species_list.size() + (size_t)size_xt) * sizeof(class spec));
sol_D[l_cell_no].spec_size = (int)species_list.size() + size_xt;
}
if (sol_D[l_cell_no].spec == NULL)
malloc_error();
for (i = 0; i < sol_D[l_cell_no].spec_size; 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 = min_dif_LM;
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;
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;
if (l_cell_no == 0)
{
por = cell_data[1].por;
por_il = cell_data[1].por_il;
}
else if (l_cell_no == count_cells + 1)
{
por = cell_data[count_cells].por;
por_il = cell_data[count_cells].por_il;
}
else
{
por = cell_data[l_cell_no].por;
por_il = cell_data[l_cell_no].por_il;
}
if (por < multi_Dpor_lim)
por = viscos_f = 0.0;
if (por_il < interlayer_Dpor_lim)
por_il = viscos_il_f = 0.0;
/*
* correct diffusion coefficient for temperature and viscosity, D_T = D_298 * viscos_298 / viscos
* modify viscosity effect: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15), SC data from Robinson and Stokes, 1959
*/
viscos = viscos_0;
/*
* put temperature factor in por_factor which corrects for porous medium...
*/
dum = viscos_0_25 / viscos;
viscos_f *= dum;
viscos_il_f *= dum;
sol_D[l_cell_no].viscos_f = dum;
count_spec = count_exch_spec = 0;
/*
* sort species by name...
*/
if (species_list.size() > 1)
{
qsort(&species_list[0], species_list.size(),
sizeof(class species_list), sort_species_name);
}
if (correct_Dw)
calc_SC();
for (i = 0; i < (int)species_list.size(); i++)
{
/*
* copy species data
*/
s_ptr = species_list[i].s;
if (s_ptr->type > HPLUS &&
!(s_ptr->type == EX && interlayer_Dflag))
continue;
//if (s_ptr->type == EX && !interlayer_Dflag)
// continue;
//if (s_ptr->type == SURF)
// continue;
if (i > 0 && strcmp(s_ptr->name, species_list[(size_t)i - 1].s->name) == 0)
continue;
//if (s_ptr == s_h2o)
// continue;
if (s_ptr->type == EX)
{
if (s_ptr->lm > min_dif_LM)
{
/* find exchanger's name, use only master exchanger 'X' */
if (species_list[i].master_s->secondary != NULL)
master_ptr = species_list[i].master_s->secondary;
else
master_ptr = species_list[i].master_s->primary;
if (s_ptr->equiv != 0.0)
dum = fabs(s_ptr->equiv) / master_ptr->total;
else
{
if (species_list[i].master_s->z == 0)
dum = 1 / master_ptr->total;
else
dum = 1;
}
name = master_ptr->elt->name;
if (strcmp(name, "X") != 0)
{
if (!warn_MCD_X)
{
snprintf(token, sizeof(token),
"MCD found more than 1 exchanger, uses X for interlayer diffusion.");
warning_msg(token);
warn_MCD_X = 1;
}
continue;
}
dum2 = s_ptr->moles * dum; /* equivalent fraction */
sol_D[l_cell_no].spec[count_spec].name = s_ptr->name;
//string_hsave(s_ptr->name);
sol_D[l_cell_no].spec[count_spec].type = EX;
sol_D[l_cell_no].spec[count_spec].c = dum2;
sol_D[l_cell_no].spec[count_spec].lg = s_ptr->lg - log10(dum);
sol_D[l_cell_no].spec[count_spec].a = dum2 * pow(10, sol_D[l_cell_no].spec[count_spec].lg);
sol_D[l_cell_no].exch_total = master_ptr->total;
if (transport_step == 0 && !x_max_done)
{
x_max_done = true;
dum = master_ptr->total / Utilities::Rxn_find(Rxn_solution_map, l_cell_no)->Get_mass_water();
if (dum > sol_D[1].x_max)
sol_D[1].x_max = dum;
}
/* find the aqueous species in the exchange reaction... */
for (i2 = 0; (s_ptr->rxn.token[i2].s != NULL); i2++)
{
if ((s_ptr2 = s_ptr->rxn.token[i2].s)->type == AQ)
break;
}
/* copy its name and Dw and charge... */
sol_D[l_cell_no].spec[count_spec].aq_name = s_ptr2->name;
//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;
else
{
if (s_ptr2->dw_t)
{
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr2->dw *
exp(s_ptr2->dw_t / 298.15 - s_ptr2->dw_t / tk_x) * viscos_il_f;
sol_D[l_cell_no].spec[count_spec].dw_t = s_ptr2->dw_t;
}
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_dif_LM)
{
if (implicit/* && l_cell_no < count_cells + 2 */)
{
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;
if (i3 + count_spec + 1 > sol_D[l_cell_no].spec_size)
{
sol_D[l_cell_no].spec = (class spec *) PHRQ_realloc(sol_D[l_cell_no].spec,
((size_t)i3 + count_spec + 1 + (size_t)size_xt) * sizeof(class spec));
if (sol_D[l_cell_no].spec == NULL)
malloc_error();
sol_D[l_cell_no].spec_size = i3 + count_spec + 1 + size_xt;
}
for (; i1 < i2; i1++) // i1 is loop variable
{
// memmove(&sol_D[l_cell_no].spec[i1], &sol_D[ref_cell].spec[i1], sizeof(class spec));
sol_D[l_cell_no].spec[i1] = sol_D[ref_cell].spec[i1];
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 = min_dif_LM;
sol_D[l_cell_no].spec[i1].lg = -0.04;
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 = (class spec *) PHRQ_realloc(sol_D[l_cell_no].spec,
(count_spec + (size_t)size_xt) * sizeof(class 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].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;
sol_D[l_cell_no].spec[count_spec].z = s_ptr->z;
if (s_ptr->dw == 0)
sol_D[l_cell_no].spec[count_spec].Dwt = default_Dw * viscos_f;
else
{
if (s_ptr->dw_t)
{
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw *
exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15) * viscos_f;
sol_D[l_cell_no].spec[count_spec].dw_t = s_ptr->dw_t;
}
else
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw * viscos_f;
}
if (correct_Dw)
{
//calc_SC(); // removed that neutral species are corrected as if z = 1, but is viscosity-dependent
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw_corr * viscos_f;
}
if (l_cell_no <= count_cells + 1 && sol_D[l_cell_no].spec[count_spec].Dwt * pow(por, multi_Dn) > diffc_max)
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/* && l_cell_no < count_cells + 2 */)
{
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 = (class spec *) PHRQ_realloc(sol_D[i1].spec,
((size_t)i2 + (size_t)size_xt) * sizeof(class spec));
if (sol_D[i1].spec == NULL)
malloc_error();
sol_D[i1].spec_size = i2 + size_xt;
}
i2--;
for (; i2 > i3; i2--) // i2 is loop variable
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(class spec));
sol_D[i1].spec[i2] = sol_D[l_cell_no].spec[i2];
sol_D[i1].spec[i2].a = 0.0;
sol_D[i1].spec[i2].lm = min_dif_LM;
sol_D[i1].spec[i2].lg = -0.04;
sol_D[i1].spec[i2].c = 0.0;
sol_D[i1].count_spec += 1;
}
}
}
count_spec++;
}
}
sol_D[l_cell_no].count_spec = count_spec;
sol_D[l_cell_no].count_exch_spec = count_exch_spec;
if (implicit/* && l_cell_no < count_cells + 2 */ && 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 = (class spec *) PHRQ_realloc(sol_D[l_cell_no].spec,
((size_t)i3 + (size_t)size_xt) * sizeof(class 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[ref_cell].spec[i1], sizeof(class spec));
sol_D[l_cell_no].spec[i1] = sol_D[ref_cell].spec[i1];
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 = min_dif_LM;
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);
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/* && l_cell_no < count_cells + 2 */)
{
for (i = 0; i < count_spec; i++)
{
//name = sol_D[l_cell_no].spec[i].name;
if (strcmp(sol_D[l_cell_no].spec[i].name, sol_D[ref_cell].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[ref_cell].spec[i].name);
error_msg(error_string, CONTINUE);
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
sort_species_name(const void *ptr1, const void *ptr2)
/* ---------------------------------------------------------------------- */
{
const class species_list *nptr1, *nptr2;
nptr1 = (const class species_list *) ptr1;
nptr2 = (const class species_list *) ptr2;
return (strcmp(nptr1->s->name, nptr2->s->name));
}
/* ---------------------------------------------------------------------- */
void Phreeqc::
diffuse_implicit(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.
// Version 3.5.2 includes diffusion in 1 stagnant layer in the implicit calc'n.
// 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, icell, cp, comp;
// ifirst = (bcon_first == 2 ? 1 : 0); ilast = (bcon_last == 2 ? count_cells - 1 : count_cells);
int ifirst, ilast;
int i_1, i0, i1, i2 = 0;
//double mfr, mfr1, max_b = 0, b, grad, dVc, j_0e, min_dif_M = pow(10, min_dif_LM);
double mfr, mfr1, grad, dVc, j_0e, min_dif_M = pow(10, min_dif_LM);
LDBLE dum1, dum2, dum_stag = 0.0, min_mol;
LDBLE dum = 0;
//cxxSurfaceCharge * charge_ptr = NULL;
if (stagnant > 1)
stagnant = 0;
cxxSolution *sptr1, *sptr2, *sptr_stag;
std::vector<int> n_solution;
std::vector<LDBLE> fraction;
std::map<int, std::map<std::string, double> > ::iterator it1;
std::map<std::string, double> ::iterator it2;
// c= count_cells, c1=(c+1)= end boundary-cell, c_1=(c-1), c2=(c+2)= first stagnant cell, cc1=(c+c+1)= last stagnant cell
int c = count_cells, c1 = c + 1, c2 = c + 2, c_1 = c - 1, cc = c + stagnant * c, cc1 = cc + 1;
comp = sol_D[1].count_spec - sol_D[1].count_exch_spec;
cell_J_ij.clear();
for (i = 0; i <= count_cells; i++)
{
std::map<std::string, class J_ij_save> J_map;
for (cp = 0; cp < comp; cp++)
{
class J_ij_save J_save;
J_map[sol_D[1].spec[cp].name] = J_save;
}
cell_J_ij[i] = J_map;
}
if (heat_nmix)
comp += 1;
if (Ct2 == NULL)
Ct2 = (LDBLE *) PHRQ_malloc((count_cells + 2 + (size_t)stagnant * count_cells) * sizeof(LDBLE));
if (Ct2 == NULL) malloc_error();
if (l_tk_x2 == NULL)
l_tk_x2 = (LDBLE *) PHRQ_malloc((count_cells + 2 + (size_t)stagnant * count_cells) * sizeof(LDBLE));
if (l_tk_x2 == NULL) malloc_error();
if (A == NULL)
{
A = (LDBLE **)PHRQ_malloc((count_cells + 2 + (size_t)stagnant * count_cells) * sizeof(LDBLE *));
if (A == NULL) malloc_error();
for (i = 0; i < count_cells + 2 + stagnant * count_cells; i++)
{
if (stagnant)
A[i] = (LDBLE *)PHRQ_calloc((2 * count_cells + 2), sizeof(LDBLE));
else
A[i] = (LDBLE *)PHRQ_malloc(3 * sizeof(LDBLE));
if (A[i] == NULL) malloc_error();
}
}
if (LU == NULL)
{
LU = (LDBLE **)PHRQ_malloc((count_cells + 2 + (size_t)stagnant * count_cells) * sizeof(LDBLE *));
if (LU == NULL) malloc_error();
for (i = 0; i < count_cells + 2 + stagnant * count_cells; i++)
{
if (stagnant)
LU[i] = (LDBLE *)PHRQ_calloc((2 * count_cells + 2), sizeof(LDBLE));
else
LU[i] = (LDBLE *)PHRQ_malloc(3 * sizeof(LDBLE));
if (LU[i] == NULL) malloc_error();
}
}
if (mixf == NULL)
{
mixf = (LDBLE **)PHRQ_malloc((count_cells + 2) * sizeof(LDBLE *));
if (mixf == NULL) malloc_error();
for (i = 0; i < count_cells + 2; i++)
{
mixf[i] = NULL;
}
}
if (stagnant)
{
if (mixf_stag == NULL)
{
mixf_stag = (LDBLE **)PHRQ_malloc((count_cells + 2) * sizeof(LDBLE *));
if (mixf_stag == NULL) malloc_error();
for (i = 0; i < count_cells + 2; i++)
{
mixf_stag[i] = NULL;
}
}
}
if (comp + 2 > mixf_comp_size)
{
for (i = 0; i < count_cells + 2; i++)
{
if (mixf[i] == NULL)
mixf[i] = (LDBLE *)PHRQ_malloc(((size_t)comp + 2) * sizeof(LDBLE));
else
mixf[i] = (LDBLE *)PHRQ_realloc(mixf[i], ((size_t)comp + 2) * sizeof(LDBLE));
if (mixf[i] == NULL) malloc_error();
if (stagnant)
{
if (mixf_stag[i] == NULL)
mixf_stag[i] = (LDBLE *)PHRQ_malloc(((size_t)comp + 2) * sizeof(LDBLE));
else
mixf_stag[i] = (LDBLE *)PHRQ_realloc(mixf_stag[i], ((size_t)comp + 2) * sizeof(LDBLE));
if (mixf_stag[i] == NULL) malloc_error();
for (cp = 0; cp < comp; cp++)
mixf_stag[i][cp] = 0;
}
}
mixf_comp_size = comp + 2;
}
if (dV_dcell)
find_current = 1;
// obtain b_ij...
dummy = default_Dw * pow(multi_Dpor, multi_Dn) * multi_Dpor;
for (i = 0; i <= count_cells + 1; i++)
{
if (!heat_nmix)
{
if (i <= count_cells)
l_tk_x2[i] = (sol_D[i].tk_x + sol_D[i + 1].tk_x) / 2;
//if (stagnant && i > 0 && i <= c)
// l_tk_x2[i + c1] = (sol_D[i].tk_x + sol_D[i + c1].tk_x) / 2;
}
if (stagnant)
{
use.Set_mix_ptr(Utilities::Rxn_find(Rxn_mix_map, i));
if (use.Get_mix_ptr() == NULL)
{
i1 = 0;
if (i == 0 && (sptr1 = Utilities::Rxn_find(Rxn_solution_map, c2)) != NULL)
{
// the 1st boundary solution diffuses into 1st stagnant cell (=c2) by the ratio immobile_kgw in cell c2 / mobile_kgw in cell 1
find_J(0, 1, 1, 1, 0);
ct[c2].kgw = sptr1->Get_mass_water();
for (cp = 0; cp < comp; cp++)
{
if (!heat_nmix || cp < comp - 1)
mixf_stag[i][cp] = DDt * ct[i].v_m[cp].b_ij * ct[c2].kgw / ct[1].kgw;
else if (heat_nmix && cp == comp - 1)
mixf_stag[i][cp] = mixf_stag[i][0] * heat_diffc / tempr / sol_D[i].spec[0].Dwt;
}
i1 = c2;
}
else if (i == c1 && (sptr1 = Utilities::Rxn_find(Rxn_solution_map, cc1)) != NULL)
{
// the end boundary solution diffuses into stagnant cell cc1
find_J(c, c1, 1, 1, 0);
ct[cc1].kgw = sptr1->Get_mass_water();
for (cp = 0; cp < comp; cp++)
{
if (!heat_nmix || cp < comp - 1)
mixf_stag[i][cp] = DDt * ct[c].v_m[cp].b_ij * ct[cc1].kgw / ct[c].kgw;
else if (heat_nmix && cp == comp - 1)
mixf_stag[i][cp] = mixf_stag[i][0] * heat_diffc / tempr / sol_D[c].spec[0].Dwt;
}
i1 = cc1;
}
if (i1)
{
find_J(i, i1, 1, 1, 1);
}
}
else
{
use.Get_mix_ptr()->Vectorize(n_solution, fraction);
for (i1 = 0; i1 < (int) use.Get_mix_ptr()->Get_mixComps().size(); i1++)
{
if (n_solution[i1] > count_cells && n_solution[i1] <= cc1)
{
mfr = fraction[i1] / nmix;
find_J(i, n_solution[i1], mfr, 1, 1);
mfr *= (2.0 / dummy);
for (cp = 0; cp < comp; cp++)
{
if (!heat_nmix || cp < comp - 1)
mixf_stag[i][cp] = ct[i].v_m[cp].b_ij * mfr;
else if (heat_nmix && cp == comp - 1)
mixf_stag[i][cp] = mixf_stag[i][0] * heat_diffc / tempr / sol_D[i].spec[0].Dwt;
}
}
}
}
}
if (i <= count_cells)
find_J(i, i + 1, 1, 1, 0);
}
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 <= cc1; i++)
{
if (heat_nmix && cp == comp - 1)
Ct2[i] = sol_D[i].tk_x;
else
Ct2[i] = sol_D[i].spec[cp].c;
}
// 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;
mfr1 = mixf[1][cp] = DDt * ct[1].v_m[cp].b_ij;
}
if (bcon_first == 2)
{
A[1][0] = 0; A[1][1] = 1 + mfr1; A[1][2] = -mfr1;
if (stagnant)
{
A[0][0] = A[c2][c2] = 1;
if (mixf_stag[1][cp])
{
A[1][1] += mixf_stag[1][cp];
A[1][c2] = A[c2][1] = -mixf_stag[1][cp];
A[c2][c2] += mixf_stag[1][cp];
}
}
else
{
A[0][1] = 1; A[0][2] = 0;
}
}
else
{
A[1][0] = -mfr; A[1][1] = 1 + mfr + mfr1; A[1][2] = -mfr1;
if (stagnant)
{
if (dV_dcell)
{
A[0][0] = 1 + mfr; A[0][1] = -mfr;
}
else
A[0][0] = 1;
A[c2][c2] = 1;
if (mixf_stag[0][cp])
{
if (dV_dcell)
{
A[0][0] += mixf_stag[0][cp]; A[0][c2] = -mixf_stag[0][cp];
}
A[c2][0] = -mixf_stag[0][cp];
A[c2][c2] += mixf_stag[0][cp];
}
if (mixf_stag[1][cp])
{
A[1][1] += mixf_stag[1][cp];
A[1][c2] = A[c2][1] = -mixf_stag[1][cp];
A[c2][c2] += mixf_stag[1][cp];
}
}
else
{
if (dV_dcell)
{
A[0][1] = 1 + mfr; A[0][2] = -mfr;
}
else
{
A[0][1] = 1; A[0][2] = 0;
}
}
}
if (heat_nmix && cp == comp - 1)
{
mfr = mixf[c_1][cp] = heat_mix_array[c_1];
mfr1 = mixf[count_cells][cp] = heat_mix_array[count_cells];
}
else
{
mfr = mixf[c_1][cp] = DDt * ct[c_1].v_m[cp].b_ij;
mfr1 = mixf[count_cells][cp] = DDt * ct[count_cells].v_m[cp].b_ij;
}
if (bcon_last == 2)
{
if (stagnant)
{
A[count_cells][c_1] = -mfr; A[count_cells][count_cells] = 1 + mfr;
A[c1][c1] = A[cc1][cc1] = 1;
if (mixf_stag[count_cells][cp])
{
A[count_cells][count_cells] += mixf_stag[count_cells][cp];
A[count_cells][cc1] = A[cc1][count_cells] = -mixf_stag[count_cells][cp];
A[cc1][cc1] += mixf_stag[count_cells][cp];
}
}
else
{
A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr; A[count_cells][2] = 0;
A[c1][0] = 0; A[c1][1] = 1;
}
}
else
{
if (stagnant)
{
A[count_cells][c_1] = -mfr; A[count_cells][count_cells] = 1 + mfr + mfr1;
A[count_cells][c1] = -mfr1;
if (dV_dcell)
{
A[c1][count_cells] = -mfr1; A[c1][c1] = 1 + mfr1;
}
else
A[c1][c1] = 1;
A[cc1][cc1] = 1;
if (mixf_stag[count_cells][cp])
{
A[count_cells][count_cells] += mixf_stag[count_cells][cp];
A[count_cells][cc1] = A[cc1][count_cells] = -mixf_stag[count_cells][cp];
A[cc1][cc1] += mixf_stag[count_cells][cp];
}
if (mixf_stag[c1][cp])
{
if (dV_dcell)
{
A[c1][c1] += mixf_stag[c1][cp]; A[c1][cc1] = -mixf_stag[c1][cp];
}
A[cc1][c1] = -mixf_stag[c1][cp];
A[cc1][cc1] += mixf_stag[c1][cp];
}
}
else
{
A[count_cells][0] = -mfr; A[count_cells][1] = 1 + mfr + mfr1; A[count_cells][2] = -mfr1;
if (dV_dcell)
{
A[c1][0] = -mfr1; A[c1][1] = 1 + mfr1;
}
else
{
A[c1][0] = 0; A[c1][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;
}
if (stagnant && mixf_stag[i][cp])
{
A[i][i - 1] = -mfr; A[i][i] = 1 + mfr + mfr1 + mixf_stag[i][cp]; A[i][i + 1] = -mfr1;
A[i + c1][i + c1] = 1 + mixf_stag[i][cp];
A[i][i + c1] = A[i + c1][i] = -mixf_stag[i][cp];
}
else
{
A[i][0] = -mfr; A[i][1] = 1 + mfr + mfr1; A[i][2] = -mfr1;
}
}
if (stagnant)
// decompose full A in LU, bypass A[][] = 0...
{
LDBLE s;
LU[0][0] = A[0][0]; LU[1][0] = A[1][0]; LU[c2][0] = A[c2][0];
LU[0][1] = A[0][1] / LU[0][0]; LU[0][c2] = A[0][c2] / LU[0][0];
for (i = 1; i <= cc1; i++)
{
i_1 = i - 1;
if (i < c2)
{
LU[i][i] = A[i][i] - LU[i][i_1] * LU[i_1][i];
LU[i + 1][i] = A[i + 1][i];
LU[i][i + 1] = A[i][i + 1] / LU[i][i];
if (i == 1)
{
LU[c2][1] = A[c2][1] - LU[c2][0] * LU[0][1];
LU[1][c2] = (A[1][c2] - LU[1][0] * LU[0][c2]) / LU[1][1];
}
else
{
for (i0 = 0; i0 <= i_1; i0++)
{
if (i0 == i_1 && i0 < c)
{
LU[c2 + i0][i] = A[c2 + i0][i];
LU[i][c2 + i0] = A[i][c2 + i0] / LU[i][i];
}
else if (i0 < c)
{
LU[c2 + i0][i] = -LU[c2 + i0][i_1] * LU[i_1][i];
LU[i][c2 + i0] = -LU[i][i_1] * LU[i_1][c2 + i0] / LU[i][i];
}
if (i > c && i0 == c_1)
{
LU[c2 + i0][i] += A[c2 + i0][i];
LU[i][c2 + i0] += A[i][c2 + i0] / LU[i][i];
}
}
}
}
else
{
// the i = count_cells + 2 and higher rows are filled from i0 onwards, need to subtract the L*U terms for the ith column of L.
for (i0 = i; i0 <= cc1; i0++)
{
s = 0;
if (i0 == i)
{
i2 = (i == c2 ? 0 : i - c1);
for (i1 = i2; i1 <= i_1; i1++)
s += LU[i0][i1] * LU[i1][i];
LU[i][i] = A[i][i] - s;
i2 += (i2 == 0 ? 2 : 1);
}
else
{
for (i1 = i2; i1 <= i_1; i1++)
s += LU[i0][i1] * LU[i1][i];
LU[i0][i] = -s;
}
}
// and the ith row of U...
if (i == cc1)
break;
for (i0 = i + 1; i0 <= cc1; i0++)
{
s = 0;
if (i0 == i + 1)
i2 = i - c;
for (i1 = i2; i1 <= i_1; i1++)
s += LU[i][i1] * LU[i1][i0];
LU[i][i0] = -s / LU[i][i];
i2++;
}
}
}
Ct2[0] /= LU[0][0];
for (i = 1; i <= cc1; i++)
{
i_1 = i - 1;
if (i < c2)
s = LU[i][i_1] * Ct2[i_1];
else
{
i2 = 0;
s = 0;
for (i1 = i2; i1 <= i_1; i1++)
{
s += LU[i][i1] * Ct2[i1];
i2 += (i2 == 0 ? 2 : 1);
}
}
Ct2[i] = (Ct2[i] - s) / LU[i][i];
}
for (i = cc; i >= 0; i--)
{
s = 0;
for (i1 = i + 1; i1 <= cc1; i1++)
s += LU[i][i1] * Ct2[i1];
Ct2[i] -= s;
}
}
else
{
// 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][1];
A[i][1] -= A[i][0] * A[i - 1][2];
}
// solve Ct2 in A.Ct2 = L.U.Ct2 = Ct1, Ct1 was put in Ct2 ...
// First, find y in L.y = Ct1, put y in Ct2
Ct2[0] /= A[0][1];
for (i = 1; i <= count_cells + 1; i++)
Ct2[i] = (Ct2[i] - A[i][0] * Ct2[i - 1]) / A[i][1];
// Now obtain Ct2 in U.Ct2 = 'y' ...
for (i = count_cells; i >= 0; i--)
Ct2[i] -= A[i][2] * Ct2[i + 1];
}
// Moles transported by concentration gradient from cell [i] to [i + 1] go in tot1,
// moles by stagnant exchange from cell [i] to [i + c1] go in tot_stag
// 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);
if (stagnant && i > 0 && i < c1 && mixf_stag[i][cp])
{
i1 = i + c1;
cell_data[i1].temp = Ct2[i1] - 273.15;
sptr2 = Utilities::Rxn_find(Rxn_solution_map, i1);
sptr2->Set_tc(Ct2[i1] - 273.15);
}
continue;
}
else if (i == ilast + 1 && (!stagnant || !mixf_stag[i][cp]))
continue;
if (!cp)
{
ct[i].J_ij_sum = 0;
ct[i].Dz2c_stag = 0;
if (i < ilast + 1)
current_cells[i].ele = current_cells[i].dif = 0;
}
for (i0 = stagnant; i0 >= 0; i0--)
{
i1 = (i0 == 0 ? i + 1 : i == 0 ? c2 : i == c1 ? cc1 : i + c1);
if (Utilities::Rxn_find(Rxn_solution_map, i1) == NULL)
continue;
ct[i].J_ij[cp].name = sol_D[i].spec[cp].name;
ct[i].J_ij[cp].charge = ct[i].v_m[cp].z;
grad = (Ct2[i1] - Ct2[i])/* * ct[i].v_m[cp].D*/; // .D has the d{lg(gamma)} / d{lg(molal)} correction
if (!i0)
{
ct[i].J_ij[cp].tot1 = -mixf[i][cp] * grad;
std::map<int, std::map<std::string, class J_ij_save> >::iterator
cell_iter = cell_J_ij.find(i);
if (cell_iter == cell_J_ij.end())
continue;
std::map< std::string, class J_ij_save >::iterator
s_iter = cell_iter->second.find(sol_D[1].spec[cp].name);
s_iter->second.flux_c = ct[i].J_ij[cp].tot1;
s_iter->second.flux_t = ct[i].J_ij[cp].tot1;
}
else
ct[i].J_ij[cp].tot_stag = -mixf_stag[i][cp] * grad;
if (ct[i].v_m[cp].z)
{
if (i == 0)
ct[i].v_m[cp].zc = ct[i].v_m[cp].z * Ct2[i1];
else if (i == ilast)
{
if (!i0)
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[i1];
}
else
ct[i].v_m[cp].zc = ct[i].v_m[cp].z * (Ct2[i] + Ct2[i1]) / 2;
if (i0)
{
mixf_stag[i][cp] *= ct[i].v_m[cp].zc;
ct[i].J_ij_sum += ct[i].v_m[cp].z * ct[i].J_ij[cp].tot_stag;
ct[i].Dz2c_stag -= ct[i].v_m[cp].z * mixf_stag[i][cp];
}
else if (i < ilast + 1)
{
mixf[i][cp] *= ct[i].v_m[cp].zc;
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];
}
}
}
}
// 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);
const char* cptr = temp_name;
count_elts = 0;
get_elts_in_species(&cptr, 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 || fix_current)
{
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 = 0e-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);
//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))))
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) / ((double)count_cells + 1 - (double)i);
}
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++)
{
if (!ct[ifirst].v_m[cp].z)
continue;
for (i = ifirst; i <= ilast + stagnant; i++)
{
if (i < ilast + 1)
{
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] * dVc;
std::map<int, std::map<std::string, class J_ij_save> >::iterator
cell_iter = cell_J_ij.find(i);
std::map< std::string, class J_ij_save >::iterator
s_iter = cell_iter->second.find(sol_D[1].spec[cp].name);
s_iter->second.flux_t = ct[i].J_ij[cp].tot1;
}
if (stagnant && ct[i].Dz2c_stag)
{
ct[i].J_ij[cp].tot_stag += (mixf_stag[i][cp] * ct[i].J_ij_sum / ct[i].Dz2c_stag);
}
}
}
current_A = current_x / DDt * F_C_MOL;
for (i = ifirst; i <= ilast + stagnant + (bcon_last == 2 ? 1 : 0); i++)
{
if (i <= ilast + 1)
{
// preserve the potentials...
sptr1 = Utilities::Rxn_find(Rxn_solution_map, i);
sptr1->Set_potV(cell_data[i].potV);
}
// 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 = (class M_S *) free_check_null(ct[i].m_s);
if (ct[i].m_s == NULL)
{
ct[i].m_s = (class M_S *) PHRQ_malloc((count_m_s + 5) * sizeof(class 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 = (class M_S *) PHRQ_realloc(ct[i].m_s, (count_m_s + 5) * sizeof(class M_S));
ct[i].m_s_size = count_m_s + 5;
}
if (ct[i].m_s == NULL)
malloc_error();
std::set <std::string> ::iterator it = dif_els_names.begin();
for (i1 = 0; i1 < count_m_s; i1++)
{
ct[i].m_s[i1].tot1 = ct[i].m_s[i1].tot2 = ct[i].m_s[i1].tot_stag = 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, stagnant);
}
/*
* 3. find the solutions in the column, add or subtract the moles...
*/
cxxNameDouble::iterator kit;
int if1, il1, incr;
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)
{
min_mol = min_dif_M * ct[icell].kgw;
if (min_mol < 1e-13)
min_mol = 1e-13;
dum1 = dum2 = 0;
sptr1 = Utilities::Rxn_find(Rxn_solution_map, icell);
sptr2 = Utilities::Rxn_find(Rxn_solution_map, icell + 1);
if (stagnant && mixf_stag[icell][cp])
{
i1 = (icell == 0 ? c1 + 1 : icell == c1 ? cc1 : icell + c1);
sptr_stag = Utilities::Rxn_find(Rxn_solution_map, i1);
}
else
sptr_stag = NULL;
//if (!cp)
//{
// ct[icell].J_ij_sum = ct[icell + 1].J_ij_sum = 0.0;
// if (sptr_stag)
// ct[i1].J_ij_sum = 0.0;
//}
if (!strcmp(ct[icell].m_s[cp].name, "H"))
{
dummy = ct[icell].m_s[cp].tot1;
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_h(sptr1->Get_total_h() - dummy);
if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
sptr2->Set_total_h(sptr2->Get_total_h() + dummy);
if (sptr_stag)
{
dummy = ct[icell].m_s[cp].tot_stag;
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_h(sptr1->Get_total_h() - dummy);
if (icell == c)
{
// mix the boundary solution with the stagnant column end-cell
dummy += ct[icell + 1].m_s[cp].tot_stag;
if (dV_dcell || fix_current)
sptr2->Set_total_h(sptr2->Get_total_h() - ct[icell + 1].m_s[cp].tot_stag);
}
sptr_stag->Set_total_h(sptr_stag->Get_total_h() + dummy);
}
continue;
}
if (!strcmp(ct[icell].m_s[cp].name, "O"))
{
dummy = ct[icell].m_s[cp].tot1;
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_o(sptr1->Get_total_o() - dummy);
if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
sptr2->Set_total_o(sptr2->Get_total_o() + dummy);
if (sptr_stag)
{
dummy = ct[icell].m_s[cp].tot_stag;
if (dV_dcell || fix_current || (icell > 0 && icell <= ilast))
sptr1->Set_total_o(sptr1->Get_total_o() - dummy);
if (icell == c)
{
dummy += ct[icell + 1].m_s[cp].tot_stag;
if (dV_dcell || fix_current)
sptr2->Set_total_o(sptr2->Get_total_o() - ct[icell + 1].m_s[cp].tot_stag);
}
sptr_stag->Set_total_o(sptr_stag->Get_total_o() + dummy);
}
//if (cp == count_m_s - 1) // transport the charge imbalance
//{
// sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum);
// sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum);
// if (sptr_stag)
// sptr_stag->Set_cb(sptr_stag->Get_cb() + ct[i1].J_ij_sum);
//}
continue;
}
dum1 = sptr1->Get_totals()[ct[icell].m_s[cp].name];
dum2 = sptr2->Get_totals()[ct[icell].m_s[cp].name];
if (sptr_stag)
dum_stag = sptr_stag->Get_totals()[ct[icell].m_s[cp].name];
// check for negative moles, add moles from other redox states and the donnan layer when necessary and available...
if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol &&
(dV_dcell || fix_current || (icell > 0 && icell <= ilast)))
{
dum1 = moles_from_redox_states(sptr1, ct[icell].m_s[cp].name);
if (ct[icell].dl_s > 1e-8)
{
cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell);
if (s_ptr)
{
dum1 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, ct[icell].m_s[cp].tot1 + ct[icell].m_s[cp].tot_stag - dum1 + min_mol);
}
}
sptr1->Get_totals()[ct[icell].m_s[cp].name] = dum1;
if (dum1 - ct[icell].m_s[cp].tot1 - ct[icell].m_s[cp].tot_stag < min_mol)
ct[icell].m_s[cp].tot1 = dum1 - ct[icell].m_s[cp].tot_stag - min_mol;
}
// check for negative moles, in the other cell...
ct[icell].m_s[cp].tot2 = ct[icell].m_s[cp].tot1;
dum = 0;
if (icell == c && sptr_stag && ct[c1].m_s[cp].tot_stag)
dum = ct[c1].m_s[cp].tot_stag;
if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol &&
(dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2)))
{
dum2 = moles_from_redox_states(sptr2, ct[icell].m_s[cp].name);
if (ct[icell + 1].dl_s > 1e-8)
{
cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, icell + 1);
if (s_ptr)
{
dum2 += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(ct[icell].m_s[cp].tot2 + dum2 - min_mol));
}
}
sptr2->Get_totals()[ct[icell].m_s[cp].name] = dum2;
if (dum2 + ct[icell].m_s[cp].tot2 - dum < min_mol)
ct[icell].m_s[cp].tot2 = -dum2 + min_mol + dum;
}
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 || fix_current || (icell > 0 && icell <= ilast))
{
dum = ct[icell].m_s[cp].tot1;
if (stagnant)
dum += ct[icell].m_s[cp].tot_stag;
dum1 -= dum;
sptr1->Get_totals()[ct[icell].m_s[cp].name] = (dum1 > 0 ? dum1 : min_mol);
if (dum1 < 0)
{
dum += dum1 - min_mol;
if ((it1 = neg_moles.find(icell)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
dum1 += it2->second;
els.clear();
els.insert(std::make_pair(ct[icell].m_s[cp].name, dum1));
neg_moles.erase(icell);
neg_moles.insert(std::make_pair(icell, els));
}
//ct[icell].J_ij_sum -= dum * ct[icell].m_s[cp].charge;
}
if (dV_dcell || fix_current || (icell >= 0 && icell < ilast) || (icell == ilast && bcon_last == 2))
{
dum = ct[icell].m_s[cp].tot1;
if (stagnant && icell == c && (dV_dcell || fix_current))
dum -= ct[c1].m_s[cp].tot_stag;
dum2 += dum;
sptr2->Get_totals()[ct[icell].m_s[cp].name] = (dum2 > 0 ? dum2 : min_mol);
if (dum2 < 0)
{
dum -= dum2 - min_mol;
if ((it1 = neg_moles.find(icell + 1)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
dum2 += it2->second;
els.clear();
els.insert(std::make_pair(ct[icell].m_s[cp].name, dum2));
neg_moles.erase(icell + 1);
neg_moles.insert(std::make_pair(icell + 1, els));
}
//ct[icell + 1].J_ij_sum += dum * ct[icell].m_s[cp].charge;
}
if (sptr_stag)
{
dum = ct[icell].m_s[cp].tot_stag;
if (icell == c)
dum += ct[c1].m_s[cp].tot_stag;
if (dum_stag + dum < 0)
{
dum_stag = moles_from_redox_states(sptr_stag, ct[icell].m_s[cp].name);
if (ct[i1].dl_s)
{
cxxSurface * s_ptr = Utilities::Rxn_find(Rxn_surface_map, i1);
if (s_ptr)
{
dum_stag += moles_from_donnan_layer(s_ptr, ct[icell].m_s[cp].name, -(dum + dum_stag - min_mol));
}
sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = dum_stag;
}
}
dum_stag += dum;
sptr_stag->Get_totals()[ct[icell].m_s[cp].name] = (dum_stag > 0 ? dum_stag : min_mol);
if (dum_stag < 0)
{
dum -= dum_stag - min_mol;
if ((it1 = neg_moles.find(i1)) != neg_moles.end()
&& (it2 = (els = it1->second).find(ct[icell].m_s[cp].name)) != els.end())
dum_stag += it2->second;
els.clear();
els.insert(std::make_pair(ct[icell].m_s[cp].name, dum_stag));
neg_moles.erase(i1);
neg_moles.insert(std::make_pair(i1, els));
}
//ct[i1].J_ij_sum += dum * ct[icell].m_s[cp].charge;
}
// reduce oscillations in the column-boundary cells, but not for H and O, and current_A is not adjusted...
if ((dV_dcell/* || fix_current*/) && icell == il1 - incr && dV_dcell * ct[0].m_s[cp].charge < 0 && strcmp(ct[0].m_s[cp].name, "H") && strcmp(ct[0].m_s[cp].name, "O") && c > 3 && mixrun > 1)
{
dummy = Utilities::Rxn_find(Rxn_solution_map, 0)->Get_totals()[ct[0].m_s[cp].name] / ct[0].kgw * (1 - ct[0].dl_s);
if (dummy > 1e-6)
{
sptr1 = Utilities::Rxn_find(Rxn_solution_map, 1);
sptr2 = Utilities::Rxn_find(Rxn_solution_map, 2);
dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[1].kgw * (1 - ct[1].dl_s) - dummy;
dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[2].kgw * (1 - ct[2].dl_s) - dummy;
if (dum1 / dum2 < 0 || dum1 / dum2 > 1)
{
dum = cell_data[1].mid_cell_x / cell_data[2].mid_cell_x;
//ct[1].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[1].m_s[cp].charge;
dum1 = (dummy + dum * dum2) * ct[1].kgw / (1 - ct[1].dl_s);
sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1;
//ct[1].J_ij_sum += dum1 * ct[1].m_s[cp].charge;
}
}
dummy = Utilities::Rxn_find(Rxn_solution_map, c1)->Get_totals()[ct[0].m_s[cp].name] / ct[c1].kgw * (1 - ct[c1].dl_s);
if (dummy > 1e-6)
{
sptr1 = Utilities::Rxn_find(Rxn_solution_map, c);
sptr2 = Utilities::Rxn_find(Rxn_solution_map, c_1);
dum1 = sptr1->Get_totals()[ct[0].m_s[cp].name] / ct[c].kgw * (1 - ct[c].dl_s) - dummy;
dum2 = sptr2->Get_totals()[ct[0].m_s[cp].name] / ct[c_1].kgw * (1 - ct[c_1].dl_s) - dummy;
if (dum1 / dum2 < 0 || dum1 / dum2 > 1)
{
dum = (cell_data[c].mid_cell_x - cell_data[c_1].mid_cell_x) /
(cell_data[c1].mid_cell_x - cell_data[c_1].mid_cell_x);
//ct[c].J_ij_sum -= sptr1->Get_totals()[ct[0].m_s[cp].name] * ct[c].m_s[cp].charge;
dum1 = (dummy + (1 - dum) * dum2) * ct[c].kgw / (1 - ct[c].dl_s);
sptr1->Get_totals()[ct[0].m_s[cp].name] = dum1;
//ct[c].J_ij_sum += dum1 * ct[c].m_s[cp].charge;
}
}
}
//if (cp == count_m_s - 1)
//{
// sptr1->Set_cb(sptr1->Get_cb() + ct[icell].J_ij_sum);
// sptr2->Set_cb(sptr2->Get_cb() + ct[icell + 1].J_ij_sum);
// if (sptr_stag)
// sptr_stag->Set_cb(sptr_stag->Get_cb() + ct[i1].J_ij_sum);
//}
}
}
return;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
moles_from_redox_states(cxxSolution *sptr, const char *name)
/* ---------------------------------------------------------------------- */
{
int length, length2;
cxxNameDouble::iterator kit;
LDBLE dum = 0;
length = (int)strlen(name);
for (kit = sptr->Get_totals().begin(); kit != sptr->Get_totals().end(); kit++)
{
length2 = (int)(size_t)strcspn(kit->first.c_str(), "(");
if (length == length2 && !strncmp(name, kit->first.c_str(), length2))
{
dum += kit->second; kit->second = 0;
}
}
return(dum);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
add_MCD_moles(LDBLE moles, LDBLE min_mol, int i, cxxSolution *sptr, const char *name)
/* ---------------------------------------------------------------------- */
{
cxxNameDouble::iterator kit;
std::map<int, std::map<std::string, double> > ::iterator it1;
std::map<std::string, double> ::iterator it2;
LDBLE dum = sptr->Get_totals()[name];
if (!dum)
dum = moles_from_redox_states(sptr, name);
if ((it1 = neg_moles.find(i)) != neg_moles.end()
&& (it2 = (els = it1->second).find(name)) != els.end())
{
dum += it2->second;
neg_moles.erase(it1);
els.erase(it2);
neg_moles.insert(std::make_pair(i, els));
}
dum += moles;
if (dum < -min_mol)
{
if (ct[i].dl_s)
{
cxxSurface *su_ptr = Utilities::Rxn_find(Rxn_surface_map, i);
if (su_ptr != NULL)
dum += moles_from_donnan_layer(su_ptr, name, -dum + min_mol);
}
}
sptr->Get_totals()[name] = (dum > 0 ? dum : 0e-16);
if (dum < -min_mol)
{
els.insert(std::make_pair(name, dum));
if ((it1 = neg_moles.find(i)) != neg_moles.end())
neg_moles.erase(it1);
neg_moles.insert(std::make_pair(i, els));
}
return(dum);
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
moles_from_donnan_layer(cxxSurface *sptr, const char *name, LDBLE moles_needed)
/* ---------------------------------------------------------------------- */
{
cxxNameDouble::iterator kit;
LDBLE dum = 0;
cxxSurfaceCharge * charge_ptr = NULL;
for (size_t j = 0; j < sptr->Get_surface_charges().size(); j++)
{
if (sptr->Get_dl_type() == cxxSurface::DONNAN_DL)
{
charge_ptr = &(sptr->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(), name) == 0)
{
if (kit->second > moles_needed)
{
dum += moles_needed; kit->second -= moles_needed;
}
else
{
dum += kit->second; kit->second = 0;
}
}
}
}
}
return(dum);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
multi_D(LDBLE DDt, int mobile_cell, int stagnant)
/* ---------------------------------------------------------------------- */
{
/*
* 1. determine mole transfer (mol/s) of solute species for the interface between 2 cells.
* 2. sum up as mole transfer of master_species
* 3. add moles of master_species to the 2 cells
* NOTE. Define the water content of stagnant cells relative to the
* mobile cell (with, for example, 1 kg water)
* Define properties of each interface only 1 time with MIX.
* If an electrical field is applied (dV_dcell != 0), the currents j_i = current_cells[i].ele + dif (C * s)
are calculated for all cells. Then the ele part from cell 0 -> 1 is calculated:
current_x = (j_0d + j_0e) = j_0 = j_1 = ... = j_i
j_0e * (R0 + R1 + ...) + (j_0d - j_1d) * R1 + ... + (j_0d - j_id) * Ri = Vtot
or
j_0e * Sum_R + Sum_Rd = Vtot.
Ri = dV_dcell / j_ie, the relative cell resistance.
Solve j_0e, find (V1 - V0) = j_0e * R0. j_1e = current_x - j_1d, find (V2 - V1) = j_1e * R1, etc.
*/
int icell, jcell, i, l, n, length, length2, il_calcs;
int i1, loop_f_c;
int first_c, last_c, last_c2 = 0;
char token[MAX_LENGTH];
LDBLE mixf, temp;
LDBLE dVtemp = 0.0;
if (dV_dcell && stagnant)
{
dVtemp = dV_dcell;
dV_dcell = 0;
}
cell_J_ij.clear();
icell = jcell = -1;
first_c = last_c = -1;
il_calcs = -1;
current_x = sum_R = sum_Rd = 0.0;
if (dV_dcell)
find_current = loop_f_c = 1; // calculate J_ij once for all cells, loop to dV_dcell2.
else
find_current = loop_f_c = 0;
for (int f_c = 0; f_c <= loop_f_c; f_c++)
{
for (n = 0; n <= (stagnant ? stag_data.count_stag : 0); n++) // allow for stagnant cell mixing with higher cells in the layer
{
icell = mobile_cell + 1 + n * count_cells;
if (stagnant)
{
if (n == 0)
icell -= 1;
else if (mobile_cell == 0 || mobile_cell == count_cells + 1)
continue;
/*
* find the mix ptr for icell and go along the cells that mix with it
*/
use.Set_mix_ptr(Utilities::Rxn_find(Rxn_mix_map, icell));
if (use.Get_mix_ptr() == NULL)
{
if (first_c < 0)
first_c = icell;
if (last_c2 < icell)
last_c2 = icell;
continue;
}
first_c = 0;
last_c = (int) (use.Get_mix_ptr()->Get_mixComps().size() - 1);
}
else
{ /* regular column... */
if (bcon_first == 1)
first_c = 0;
else
first_c = 1;
if (bcon_last == 1)
last_c = count_cells;
else
last_c = count_cells - 1;
last_c2 = (dV_dcell ? count_cells + 1: count_cells);
}
for (i = first_c; i <= last_c; i++)
{
if (stagnant)
{
std::vector<int> n_solution;
std::vector<LDBLE> fraction;
(use.Get_mix_ptr())->Vectorize(n_solution, fraction);
jcell = n_solution[i];
if (jcell <= icell)
continue;
if (jcell >= all_cells || jcell < 0)
{
error_string = sformatf(
"Multi_D asked for diffusion from cell %d to %d, %d is beyond the number of cells", icell, jcell, jcell);
error_msg(error_string, CONTINUE);
}
if (jcell > last_c2)
last_c2 = jcell;
mixf = (nmix ? fraction[i] / nmix : fraction[i]);
}
else
{ /* regular column... */
icell = i;
jcell = i + 1;
mixf = 1.0;
}
if (dV_dcell)
{
tk_x2 = (sol_D[icell].tk_x + sol_D[jcell].tk_x) / 2;
}
/*
* 1. obtain J_ij...
*/
il_calcs = (int) find_J(icell, jcell, mixf, DDt, stagnant);
if (find_current)
{
if (i < last_c)
continue;
else
{
LDBLE dVc, j_0e;
// distribute dV_dcell according to relative resistance, calculate current_x if not fixed
j_0e = (dV_dcell * count_cells - sum_Rd) / sum_R;
current_x = j_0e + current_cells[0].dif;
if (fix_current)
{
int sign = (current_x >= 0 ? 1 : -1);
current_x = sign * fix_current / F_C_MOL;
j_0e = current_x - current_cells[0].dif;
}
dVc = j_0e * current_cells[0].R;
cell_data[1].potV = cell_data[0].potV + dVc;
for (i1 = 1; i1 < count_cells; i1++)
{
dVc = current_cells[i1].R * (current_x - current_cells[i1].dif);
cell_data[i1 + 1].potV = cell_data[i1].potV + dVc;
}
if (fix_current)
{
dVc = current_cells[i1].R * (current_x - current_cells[i1].dif);
cell_data[i1 + 1].potV = cell_data[i1].potV + dVc;
}
find_current = 0;
continue;
}
}
if (!ct[icell].J_ij_count_spec)
continue;
/*
* 2. sum up the primary or secondary master_species
*/
if (!il_calcs)
{
tot1_h = tot1_o = tot2_h = tot2_o = 0.0;
m_s = (class M_S *) free_check_null(m_s);
count_m_s = (ct[icell].J_ij_count_spec < count_moles_added ?
ct[icell].J_ij_count_spec : count_moles_added);
m_s = (class M_S *) PHRQ_malloc((size_t) count_m_s *
sizeof(class M_S));
if (m_s == NULL)
malloc_error();
for (i1 = 0; i1 < count_m_s; i1++)
{
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, ct[icell].J_ij_count_spec, icell, stagnant);
/*
* 3. find the solutions, add or subtract the moles...
*/
if (dV_dcell || (icell != 0 && icell != count_cells + 1))
{
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, icell));
use.Get_solution_ptr()->Set_total_h(use.Get_solution_ptr()->Get_total_h() - tot1_h);
use.Get_solution_ptr()->Set_total_o(use.Get_solution_ptr()->Get_total_o() - tot1_o);
if (dV_dcell && (icell > 0 || fix_current))
{
use.Get_solution_ptr()->Set_potV(cell_data[icell].potV);
}
for (l = 0; l < count_m_s; l++)
{
length = (int) strlen(m_s[l].name);
cxxNameDouble::iterator it;
for (it = use.Get_solution_ptr()->Get_totals().begin();
it != use.Get_solution_ptr()->Get_totals().end(); it++)
{
length2 =
(int) (size_t) strcspn(it->first.c_str(), "(");
if (strncmp(m_s[l].name, it->first.c_str(), length) == 0 && length == length2)
{
it->second -= m_s[l].tot1;
break;
}
}
if (it == use.Get_solution_ptr()->Get_totals().end())
{
use.Get_solution_ptr()->Get_totals()[m_s[l].name] = -m_s[l].tot1;
}
}
}
if (dV_dcell || jcell != count_cells + 1)
{
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, jcell));
dummy = use.Get_solution_ptr()->Get_total_h();
use.Get_solution_ptr()->Set_total_h(dummy + tot2_h);
dummy = use.Get_solution_ptr()->Get_total_o();
use.Get_solution_ptr()->Set_total_o(dummy + tot2_o);
if (icell == count_cells && fix_current && !stagnant)
{
use.Get_solution_ptr()->Set_potV(cell_data[jcell].potV);
}
for (l = 0; l < count_m_s; l++)
{
length = (int) strlen(m_s[l].name);
cxxNameDouble::iterator it;
for (it = use.Get_solution_ptr()->Get_totals().begin();
it != use.Get_solution_ptr()->Get_totals().end(); it++)
{
length2 = (int) (size_t) strcspn(it->first.c_str(), "(");
if (strncmp(m_s[l].name, it->first.c_str(), length) == 0 && length == length2)
{
it->second += m_s[l].tot2;
break;
}
}
if (it == use.Get_solution_ptr()->Get_totals().end())
{
use.Get_solution_ptr()->Get_totals()[m_s[l].name] = m_s[l].tot2;
}
}
}
}
}
}
// check for negative conc's...
//if (stagnant)
// first_c = mobile_cell; // allow for stagnant cell mixing with boundary cell 0
for (i = first_c; i <= last_c2; i++)
{
//if (stagnant && i > first_c && i <= count_cells + first_c)
// continue;
use.Set_solution_ptr(Utilities::Rxn_find(Rxn_solution_map, i));
if (!use.Get_solution_ptr())
continue;
cxxNameDouble::iterator it;
for (it = use.Get_solution_ptr()->Get_totals().begin();
it != use.Get_solution_ptr()->Get_totals().end(); it++)
{
if (strcmp(it->first.c_str(), "H(0)") == 0)
continue;
if (strcmp(it->first.c_str(), "O(0)") == 0)
continue;
LDBLE moles = it->second;
if (moles < 0 && ct[i].dl_s)
{
use.Set_surface_ptr(Utilities::Rxn_find(Rxn_surface_map, i));
cxxSurface * s_ptr = use.Get_surface_ptr();
cxxSurfaceCharge * charge_ptr = NULL;
cxxNameDouble::iterator jit;
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 (jit = charge_ptr->Get_diffuse_layer_totals().begin(); jit != charge_ptr->Get_diffuse_layer_totals().end(); jit++)
{
if (strcmp(jit->first.c_str(), "H") == 0 || strcmp(jit->first.c_str(), "O") == 0)
continue;
if (strcmp(jit->first.c_str(), it->first.c_str()) == 0)
{
moles += jit->second;
it->second += jit->second;
jit->second = 0;
}
}
}
}
}
if (moles < 0)
{
temp = moles;
it->second = 0;
/* see if other redox states have more moles... */
length = (int) strlen(it->first.c_str());
cxxNameDouble::iterator kit;
for (kit = use.Get_solution_ptr()->Get_totals().begin();
kit != use.Get_solution_ptr()->Get_totals().end(); kit++)
{
length2 = (int) (size_t) strcspn(kit->first.c_str(), "(");
if (!strncmp(it->first.c_str(), kit->first.c_str(), length2))
{
temp += kit->second;
if (temp < 0)
{
kit->second = 0;
}
else
{
kit->second = temp;
break;
}
}
}
if (temp < -1e-12)
{
snprintf(token, sizeof(token),
"Negative concentration in MCD: added %.4e moles %s in cell %d",
(double)-temp, it->first.c_str(), i);
warning_msg(token);
for (i1 = 0; i1 < count_moles_added; i1++)
{
if (moles_added[i1].name && !strcmp(moles_added[i1].name, it->first.c_str()))
{
moles_added[i1].moles -= temp;
break;
}
else if (!moles_added[i1].moles)
{
moles_added[i1].name = string_duplicate(it->first.c_str());
moles_added[i1].moles -= temp;
break;
}
}
}
}
}
}
m_s = (class M_S *) free_check_null(m_s);
for (i = first_c; i < last_c2; i++)
{
if (stagnant && i > first_c && i <= count_cells + first_c)
continue;
ct[i].J_ij = (class J_ij *) free_check_null(ct[i].J_ij);
if (il_calcs)
ct[i].J_ij_il = (class J_ij *) free_check_null(ct[i].J_ij_il);
ct[i].v_m = (struct V_M *) free_check_null(ct[i].v_m);
}
if (dVtemp && stagnant)
{
dV_dcell = dVtemp;
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
fill_m_s(class 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
* H and O go in tot1&2_h and tot1&2_o
*/
int j, k, l;
LDBLE fraction;
const char* cptr;
for (j = 0; j < l_J_ij_count_spec; j++)
{
{
char * temp_name = string_duplicate(l_J_ij[j].name);
cptr = temp_name;
count_elts = 0;
get_elts_in_species(&cptr, 1);
free_check_null(temp_name);
}
if (implicit && stagnant < 2)
{
for (k = 0; k < count_elts; k++)
{
for (l = 0; l < count_m_s; l++)
{
if (strcmp(ct[icell].m_s[l].name, elt_list[k].elt->name) == 0)
{
fraction = fabs((double)elt_list[k].coef * l_J_ij[j].tot1) + fabs(ct[icell].m_s[l].tot1);
if (fraction)
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;
ct[icell].m_s[l].charge *= (1 - fraction);
ct[icell].m_s[l].charge += fraction * l_J_ij[j].charge;
if (stagnant)
ct[icell].m_s[l].tot_stag += elt_list[k].coef * l_J_ij[j].tot_stag;
break;
}
}
}
}
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)
{
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)
{
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;
break;
}
}
if (l == count_m_s)
{
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++;
}
}
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
void Phreeqc::
calc_b_ij(int icell, int jcell, int k, LDBLE b_i, LDBLE b_j, LDBLE g_i, LDBLE g_j, LDBLE free_i, LDBLE free_j, int stagnant)
/* ---------------------------------------------------------------------- */
{
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) * b_j * (free_j + g_j) / (b_i * (free_i + g_i) + b_j * (free_j + g_j));
// At filterends, concentrations of ions change step-wise to the DL.
// We take the harmonic mean for f_free, the average for the DL.
if (ct[icell].v_m[k].z)
{
if (!g_i && g_j)
{
ct[icell].v_m[k].b_ij = free_j * b_i * b_j / (b_i + b_j) +
b_i * (1 - free_j) / 4 + b_j * g_j / 4;
}
else if (g_i && !g_j)
ct[icell].v_m[k].b_ij = free_i * b_i * b_j / (b_i + b_j) +
b_j * (1 - free_i) / 4 + b_i * g_i / 4;
}
// for boundary cells...
if (stagnant > 1)
{ /* for a diffusion experiment with well-mixed reservoir in cell 3 and the last stagnant cell,
and with the mixf * 2 for the boundary cells in the input... */
if (icell == 3 && !g_i && g_j)
ct[icell].v_m[k].b_ij = b_j * (free_j + g_j) / 2;
else if (jcell == all_cells - 1 && !g_j && g_i)
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i) / 2;
}
else
{
if (icell == 0 || (icell == count_cells + 1 && jcell == count_cells + count_cells + 1))
ct[icell].v_m[k].b_ij = b_j * (free_j + g_j);
else if (icell == count_cells && jcell == count_cells + 1)
ct[icell].v_m[k].b_ij = b_i * (free_i + g_i);
}
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;
return;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
/* ---------------------------------------------------------------------- */
{
/* mole transfer of the individual master_species:
* Eqn 1:
* J_ij = DDt * (A_ij / lav) * (-D_i*grad(c) + D_i*z_i*c_i * SUM(D_i*z_i*grad(c)) / SUM(D_i*(z_i)^2*c_i))
* regular column, stagnant FALSE:
* D_i = temperature-corrected Dw
* A_ij = A_icell * A_jcell
* A_icell = (L porewater in i_cell / length_icell) / tort_f_icell /
* (length_icell / 2)
* lav = A_icell + A_jcell
* grad(c) is concentration difference in icell and jcell (dx is in lav),
for activity corrections see Appelo & Wersin, 2007.
** dec. 28, 2015**
included aq_dl in the harmonic mean:
J_ij = - b_i * b_j / (b_i + b_j) * (c_j - c_i) in mol/s (see ex 21 in the manual 3).
b_i = A1 / (G_i * h_i / 2) * Dw for a pore without EDL. A1 = aq1 / h_i (m^2).
with EDL (no aq_il_i in A1, for now):
t_aq1 = aq1 + aq_dl_i. A1 = t_aq1 / h_i. f_free_i = aq1 / t_aq1.
b_i_cat = A1 / (G_i * h_i / 2) * Dw * {f_free + (1 - f_free) * Bm}. Bm = Boltzmann enrichment in EDL = g_dl.
b_i_ani = A1 / (G_i * h_i / 2) * Dw * {f_free + (1 - f_free) / Bm)}.
22/2/18: now calculates diffusion through EDL's of multiple, differently charged surfaces
* stagnant TRUE:
* same eqn for J_ij, but multplies with 2 * mixf. (times 2, because mixf = A / (G_i * h_i))
* mixf_ij = mixf / (Dw / init_tort_f) / new_tort_f * new_por / init_por
* mixf is defined in MIX; Dw is default multicomponent diffusion coefficient;
* init_tort_f equals multi_Dpor^(-multi_Dn); new_pf = new tortuosity factor.
* Interlayer diffusion (IL) takes the gradient in the equivalent concentrations on X-.
surface area A for IL:
stagnant: ct[icell].mixf_il is mixf * por_il / por.
por_il = interlayer porosity, from -interlayer_D true 'por_il'.
por = free + DL porewater porosity, from -multi_D true 'multi_Dpor'.
in regular column, A is calc'd from (free + DL porewater) and cell-length.
for IL: A * por_il / por.
por_il should be entered for the cell with the maximal cec.
IL water is related to X-, thus the cec (eq/L IL water) is the same for all cells if X is difined.
IL-water = (free + DL porewater) * por_il / por.
for IL: A * aq_il / t_aq.
*/
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, 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;
il_calcs = (interlayer_Dflag ? 1 : 0);
cxxSurface *s_ptr1, *s_ptr2;
LDBLE g_i, g_j;
std::vector<cxxSurfaceCharge> s_charge_p;
std::vector<cxxSurfaceCharge> s_charge_p1;
std::vector<cxxSurfaceCharge> s_charge_p2;
std::vector<cxxSurfaceCharge>::iterator it_sc;
std::vector<cxxSurfaceComp> s_com_p;
ct[icell].dl_s = ct[jcell].dl_s = dl_aq1 = dl_aq2 = 0.0;
if (dV_dcell && !find_current)
goto dV_dcell2;
/* check for immediate return and interlayer diffusion calcs... */
//ct[icell].J_ij_sum = 0.0;
ct[icell].J_ij_count_spec = 0;
if (!il_calcs)
{
if (stagnant)
{
if (cell_data[icell].por < multi_Dpor_lim || cell_data[jcell].por < multi_Dpor_lim)
{
if (!implicit)
return (il_calcs);
}
}
else
{ /* regular column... */
if ((icell == 0 && cell_data[1].por < multi_Dpor_lim)
|| (icell == count_cells && cell_data[count_cells].por < multi_Dpor_lim)
|| (icell != 0 && icell != count_cells && (cell_data[icell].por < multi_Dpor_lim
|| cell_data[jcell].por < multi_Dpor_lim)))
{
if (dV_dcell)
{
current_cells[icell].R = -1e15;
current_cells[icell].ele = dV_dcell / current_cells[icell].R;
current_cells[icell].dif = 0;
sum_R += current_cells[icell].R;
sum_Rd += current_cells[0].dif * current_cells[icell].R;
}
if (!implicit)
return (il_calcs);
}
}
}
/* do the calcs */
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...
*/
s_ptr1 = s_ptr2 = NULL;
ct[icell].visc1 = ct[icell].visc2 = 1.0;
only_counter = FALSE;
s_ptr1 = Utilities::Rxn_find(Rxn_surface_map, icell);
if (s_ptr1 != NULL)
{
if (s_ptr1->Get_dl_type() != cxxSurface::NO_DL)
{
s_charge_p.assign(s_ptr1->Get_surface_charges().begin(), s_ptr1->Get_surface_charges().end());
s_com_p.assign(s_ptr1->Get_surface_comps().begin(), s_ptr1->Get_surface_comps().end());
if (s_ptr1->Get_only_counter_ions())
only_counter = TRUE;
ct[icell].visc1 = s_ptr1->Get_DDL_viscosity();
/* find the immobile surface charges with DL... */
for (i = 0; i < (int)s_charge_p.size(); i++)
{
for (i1 = 0; i1 < (int)s_com_p.size(); i1++)
{
if (!(s_charge_p[i].Get_name().compare(s_com_p[i1].Get_charge_name())) && !s_com_p[i1].Get_Dw())
{
dl_aq1 += s_charge_p[i].Get_mass_water();
s_charge_p1.push_back(s_charge_p[i]);
break;
}
}
}
}
}
s_ptr2 = Utilities::Rxn_find(Rxn_surface_map, jcell);
if (s_ptr2 != NULL)
{
if (s_ptr2->Get_dl_type() != cxxSurface::NO_DL)
{
s_charge_p.assign(s_ptr2->Get_surface_charges().begin(), s_ptr2->Get_surface_charges().end());
s_com_p.assign(s_ptr2->Get_surface_comps().begin(), s_ptr2->Get_surface_comps().end());
if (s_ptr2->Get_only_counter_ions())
only_counter = TRUE;
ct[icell].visc2 = s_ptr2->Get_DDL_viscosity();
for (i = 0; i < (int)s_charge_p.size(); i++)
{
for (i1 = 0; i1 < (int)s_com_p.size(); i1++)
{
if (!(s_charge_p[i].Get_name().compare(s_com_p[i1].Get_charge_name())) && !s_com_p[i1].Get_Dw())
{
dl_aq2 += s_charge_p[i].Get_mass_water();
s_charge_p2.push_back(s_charge_p[i]);
break;
}
}
}
}
}
if (!stagnant)
{
if (icell == 0)
ct[icell].visc1 = ct[icell].visc2;
else if (icell == count_cells)
ct[icell].visc2 = ct[icell].visc1;
}
/* in each cell: DL surface = mass_water_DL / (cell_length)
free pore surface = mass_water_free / (cell_length)
determine DL surface as a fraction of the total pore surface... */
t_aq1 = aq1 + dl_aq1;
t_aq2 = aq2 + dl_aq2;
f_free_i = aq1 / t_aq1;
f_free_j = aq2 / t_aq2;
if (dl_aq1 > 0)
ct[icell].dl_s = dl_aq1 / t_aq1;
if (dl_aq2 > 0)
{
ct[jcell].dl_s = dl_aq2 / t_aq2;
if (!ct[icell].dl_s)
ct[icell].dl_s = 1e-8; // used in implicit appt
}
if (il_calcs)
{
/* find interlayer porosity por_il,
make it relative to exchange capacity (mol X/L), highest X in sol_D[1].x_max (mol X / L).
Find amounts of IL water and cec. */
por_il1 = por_il2 = por_il12 = 0.0;
cec1 = cec2 = cec12 = rc1 = rc2 = 0.0;
if (icell == 0)
{
por_il1 = sol_D[0].exch_total / aq1 / sol_D[1].x_max *
cell_data[1].por_il;
por_il2 = sol_D[1].exch_total / aq2 / sol_D[1].x_max *
cell_data[1].por_il;
if (sol_D[0].exch_total > 3e-10 && sol_D[1].exch_total > 3e-10)
/* take the harmonic mean... */
por_il12 = 2 * por_il1 * por_il2 / (por_il1 + por_il2);
else
/* at column ends, take the clay... */
por_il12 = (por_il1 >= por_il2 ? por_il1 : por_il2);
aq_il2 = t_aq2 * por_il12 / cell_data[1].por;
aq_il1 = t_aq2 * por_il12 / cell_data[1].por;
}
else if (icell == count_cells)
{
por_il1 = sol_D[count_cells].exch_total / aq1 / sol_D[1].x_max *
cell_data[count_cells].por_il;
por_il2 = sol_D[count_cells + 1].exch_total / aq2 / sol_D[1].x_max *
cell_data[count_cells].por_il;
if (sol_D[count_cells].exch_total > 3e-10 && sol_D[count_cells + 1].exch_total > 3e-10)
por_il12 = 2 * por_il1 * por_il2 / (por_il1 + por_il2);
else
por_il12 = (por_il1 >= por_il2 ? por_il1 : por_il2);
aq_il1 = t_aq1 * por_il12 / cell_data[count_cells].por;
aq_il2 = t_aq2 * por_il12 / cell_data[count_cells].por;
}
else
{
por_il1 = sol_D[icell].exch_total / aq1 / sol_D[1].x_max *
cell_data[icell].por_il;
por_il2 = sol_D[jcell].exch_total / aq2 / sol_D[1].x_max *
cell_data[jcell].por_il;
if (sol_D[icell].exch_total > 3e-10 && sol_D[jcell].exch_total > 3e-10)
por_il12 = 2 * por_il1 * por_il2 / (por_il1 + por_il2);
else
por_il12 = (por_il1 >= por_il2 ? por_il1 : por_il2);
aq_il1 = t_aq1 * por_il12 / cell_data[icell].por;
aq_il2 = t_aq2 * por_il12 / cell_data[jcell].por;
}
if (por_il12 < interlayer_Dpor_lim)
il_calcs = 0;
else
{
dum = sol_D[icell].exch_total;
dum2 = sol_D[jcell].exch_total;
// the rc's are for distributing the mole transfer (later on) over X and solution.
rc1 = (dum2 > dum ? dum / dum2 : 1);
rc2 = (dum > dum2 ? dum2 / dum : 1);
if (sol_D[icell].exch_total > 3e-10)
cec1 = dum / aq_il1;
else
cec1 = 2e-10;
if (sol_D[jcell].exch_total > 3e-10)
cec2 = dum2 / aq_il2;
else
cec2 = 2e-10;
/* and the largest for calculating the mass transfer since it is usually at the boundary... */
cec12 = (cec1 > cec2 ? cec1 : cec2);
}
}
/* Find ct[icell].mixf_il for IL diffusion.
In stagnant calc's, correct mixf by default values, A = por / tort.
In regular column, A = surface area / (0.5 * x * tort)*/
tort1 = tort2 = 1.0;
ct[icell].A_ij_il = ct[icell].mixf_il = 0.0;
if (stagnant)
{
mixf /= (default_Dw * pow(multi_Dpor, multi_Dn) * multi_Dpor);
if (il_calcs)
ct[icell].mixf_il = mixf * por_il12 / interlayer_tortf;
}
if (icell == 0)
{
tort1 = tort2 = pow(cell_data[1].por, -multi_Dn);
if (stagnant)
A2 = cell_data[1].por / tort2;
else
{
A2 = t_aq2 / (cell_data[1].length * 0.5 * cell_data[1].length);
if (il_calcs && !stagnant)
ct[icell].A_ij_il = A2 * por_il12 / (cell_data[1].por * interlayer_tortf);
A2 /= tort2;
}
A1 = A2;
}
else if (icell == count_cells)
{
tort1 = tort2 = pow(cell_data[count_cells].por, -multi_Dn);
if (stagnant)
A1 = cell_data[count_cells].por / tort1;
else
{
A1 = t_aq1 / (cell_data[count_cells].length * 0.5 * cell_data[count_cells].length);
if (il_calcs && !stagnant)
ct[icell].A_ij_il = A1 * por_il12 / (cell_data[count_cells].por * interlayer_tortf);
A1 /= tort1;
}
A2 = A1;
}
else
{
tort1 = pow(cell_data[icell].por, -multi_Dn);
tort2 = pow(cell_data[jcell].por, -multi_Dn);
if (stagnant)
{
A1 = cell_data[icell].por / tort1;
A2 = cell_data[jcell].por / tort2;
}
else
{
A1 = t_aq1 / (cell_data[icell].length * 0.5 * cell_data[icell].length);
A2 = t_aq2 / (cell_data[jcell].length * 0.5 * cell_data[jcell].length);
if (il_calcs && !stagnant)
{
dum = A1 * por_il12 / (cell_data[icell].por * interlayer_tortf);
dum2 = A2 * por_il12 / (cell_data[jcell].por * interlayer_tortf);
ct[icell].A_ij_il = dum * dum2 / (dum + dum2);
}
A1 /= tort1;
A2 /= tort2;
}
}
/* diffuse... */
/*
* malloc sufficient space...
*/
k = sol_D[icell].count_spec + sol_D[jcell].count_spec;
if (ct[icell].J_ij == NULL)
{
ct[icell].J_ij = (class J_ij *) PHRQ_malloc((size_t)(k) * sizeof(class J_ij));
ct[icell].J_ij_size = k;
}
else if (k > ct[icell].J_ij_size)
{
ct[icell].J_ij = (class J_ij *) PHRQ_realloc(ct[icell].J_ij, (size_t)(k) * sizeof(class J_ij));
ct[icell].J_ij_size = k;
}
if (ct[icell].J_ij == NULL)
malloc_error();
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();
for (i = 0; i < k; i++)
{
ct[icell].J_ij[i].tot1 = ct[icell].J_ij[i].tot_stag = 0.0;
ct[icell].v_m[i].grad = 0.0;
ct[icell].v_m[i].D = 1.0; // used for gamma correction
ct[icell].v_m[i].z = 0.0;
ct[icell].v_m[i].c = 0.0;
ct[icell].v_m[i].zc = 0.0;
//ct[icell].v_m[i].Dz = 0.0;
//ct[icell].v_m[i].Dzc = 0.0;
ct[icell].v_m[i].b_ij = 0.0;
}
ct[icell].Dz2c = ct[icell].Dz2c_il = 0.0;
if (il_calcs)
{
/* also for interlayer cations */
k = sol_D[icell].count_exch_spec + sol_D[jcell].count_exch_spec;
ct[icell].J_ij_il = (class J_ij *) free_check_null(ct[icell].J_ij_il);
ct[icell].J_ij_il = (class J_ij *) PHRQ_malloc((size_t) k * sizeof(class J_ij));
if (ct[icell].J_ij_il == NULL)
malloc_error();
ct[icell].v_m_il = (struct V_M *) free_check_null(ct[icell].v_m_il);
ct[icell].v_m_il = (struct V_M *) PHRQ_malloc((size_t) k * sizeof(struct V_M));
if (ct[icell].v_m_il == NULL)
malloc_error();
for (i = 0; i < k; i++)
{
ct[icell].J_ij_il[i].tot1 = 0.0;
ct[icell].v_m_il[i].grad = 0.0;
ct[icell].v_m_il[i].D = 0.0;
ct[icell].v_m_il[i].z = 0.0;
ct[icell].v_m_il[i].c = 0.0;
ct[icell].v_m_il[i].zc = 0.0;
ct[icell].v_m_il[i].Dz = 0.0;
ct[icell].v_m_il[i].Dzc = 0.0;
ct[icell].v_m_il[i].b_ij = 0.0;
}
}
/*
* coefficients in Eqn (1)...
*/
i = j = k = k_il = 0;
i_max = sol_D[icell].count_spec;
j_max = sol_D[jcell].count_spec;
while (i < i_max || j < j_max)
{
if (j == j_max
|| (i < i_max && strcmp(sol_D[icell].spec[i].name, sol_D[jcell].spec[j].name) < 0))
{
/* species 'name' is only in icell */
if (il_calcs && sol_D[icell].spec[i].type == EX)
{
ct[icell].J_ij_il[k_il].name = sol_D[icell].spec[i].name;
ct[icell].v_m_il[k_il].D = sol_D[icell].spec[i].Dwt;
ct[icell].v_m_il[k_il].z = sol_D[icell].spec[i].z;
ct[icell].v_m_il[k_il].Dz = ct[icell].v_m_il[k_il].D * ct[icell].v_m_il[k_il].z;
dum = sol_D[icell].spec[i].c * cec12 / (2 * ct[icell].v_m_il[k_il].z);
ct[icell].v_m_il[k_il].Dzc = ct[icell].v_m_il[k_il].Dz * dum;
ct[icell].Dz2c_il += ct[icell].v_m_il[k_il].Dzc * ct[icell].v_m_il[k_il].z;
ct[icell].v_m_il[k_il].grad = -sol_D[icell].spec[i].c * cec12 / ct[icell].v_m_il[k_il].z; /* use equivalent fraction */
k_il++;
}
else
{
ct[icell].J_ij[k].name = sol_D[icell].spec[i].name;
ct[icell].v_m[k].z = sol_D[icell].spec[i].z;
ct[icell].v_m[k].grad = -sol_D[icell].spec[i].c; /* assume d log(gamma) / d log(c) = 0 */
c1 = sol_D[icell].spec[i].c / 2;
ct[icell].v_m[k].c = c1;
if (ct[icell].v_m[k].z)
ct[icell].v_m[k].zc = ct[icell].v_m[k].z * c1;
if (dV_dcell && ct[icell].v_m[k].z && !fix_current)
{
// compare diffusive and electromotive forces
dum = ct[icell].v_m[k].grad;
if (icell == 0)
dum2 = (cell_data[1].potV - cell_data[0].potV) / (cell_data[1].length / 2);
else if (icell == count_cells)
dum2 = (cell_data[count_cells + 1].potV - cell_data[count_cells].potV) / (cell_data[count_cells].length / 2);
else
dum2 = (cell_data[jcell].potV - cell_data[icell].potV) / ((cell_data[jcell].length + cell_data[icell].length) / 2);
dum2 *= F_Re3 / tk_x2 * ct[icell].v_m[k].z * c1;
if (dum + dum2 > 0)
{
// step out: no transport against the dV_dcell gradient if c = 0 in jcell...
if (i < i_max)
i++;
continue;
}
}
g_i = g_j = 0;
if (ct[icell].dl_s > 0)
{
if (dl_aq1)
{
for (it_sc = s_charge_p1.begin(); it_sc != s_charge_p1.end(); it_sc++)
{
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_i *= sol_D[icell].spec[i].erm_ddl;
}
if (dl_aq2)
{
for (it_sc = s_charge_p2.begin(); it_sc != s_charge_p2.end(); it_sc++)
{
if (ct[icell].v_m[k].z == 0 || only_counter)
{
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
else
{
if (abs(ct[icell].v_m[k].z) == 1)
// there is always H+ and OH-...
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
else
{
dum1 = it_sc->Get_mass_water() / t_aq2;
dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
g_j += pow(dum2, ct[icell].v_m[k].z) * dum1;
}
}
}
g_j *= sol_D[icell].spec[i].erm_ddl;
}
}
b_i = A1 * sol_D[icell].spec[i].Dwt;
b_j = A2;
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
b_j *= sol_D[icell].spec[i].Dwt;
else
{
dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
dum2 *= exp(sol_D[icell].spec[i].dw_t / sol_D[jcell].tk_x - sol_D[icell].spec[i].dw_t / sol_D[icell].tk_x);
dum2 *= sol_D[jcell].viscos_f;
b_j *= dum2;
}
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
k++;
}
if (i < i_max)
i++;
}
else if (i == i_max ||
(j < j_max && strcmp(sol_D[icell].spec[i].name, sol_D[jcell].spec[j].name) > 0))
{
/* species 'name' is only in jcell */
if (il_calcs && sol_D[jcell].spec[j].type == EX)
{
ct[icell].J_ij_il[k_il].name = sol_D[jcell].spec[j].name;
ct[icell].v_m_il[k_il].D = sol_D[jcell].spec[j].Dwt;
ct[icell].v_m_il[k_il].z = sol_D[jcell].spec[j].z;
ct[icell].v_m_il[k_il].Dz = ct[icell].v_m_il[k_il].D * ct[icell].v_m_il[k_il].z;
ct[icell].v_m_il[k_il].Dzc = ct[icell].v_m_il[k_il].Dz * sol_D[jcell].spec[j].c *
cec12 / (2 * ct[icell].v_m_il[k_il].z);
ct[icell].Dz2c_il += ct[icell].v_m_il[k_il].Dzc * ct[icell].v_m_il[k_il].z;
ct[icell].v_m_il[k_il].grad = sol_D[jcell].spec[j].c * cec12 / ct[icell].v_m_il[k_il].z; /* use equivalent fraction */
k_il++;
}
else
{
ct[icell].J_ij[k].name = sol_D[jcell].spec[j].name;
ct[icell].v_m[k].z = sol_D[jcell].spec[j].z;
ct[icell].v_m[k].grad = sol_D[jcell].spec[j].c; /* assume d log(gamma) / d log(c) = 0 */
c2 = sol_D[jcell].spec[j].c / 2;
ct[icell].v_m[k].c = c2;
if (ct[icell].v_m[k].z)
ct[icell].v_m[k].zc = ct[icell].v_m[k].z * c2;
if (dV_dcell && ct[icell].v_m[k].z && !fix_current)
{
// compare diffuse and electromotive forces
dum = ct[icell].v_m[k].grad;
if (icell == 0)
dum2 = (cell_data[1].potV - cell_data[0].potV) / (cell_data[1].length / 2);
else if (icell == count_cells)
dum2 = (cell_data[count_cells + 1].potV - cell_data[count_cells].potV) / (cell_data[count_cells].length / 2);
else
dum2 = (cell_data[jcell].potV - cell_data[icell].potV) / ((cell_data[jcell].length + cell_data[icell].length) / 2);
dum2 *= F_Re3 / tk_x2 * ct[icell].v_m[k].z * c2;
// don't transport unavailable moles against the gradient
if (dum + dum2 < 0)
{
// step out...
if (j < j_max)
j++;
continue;
}
}
g_i = g_j = 0;
if (ct[icell].dl_s > 0)
{
if (dl_aq1)
{
for (it_sc = s_charge_p1.begin(); it_sc != s_charge_p1.end(); it_sc++)
{
if (ct[icell].v_m[k].z == 0 || only_counter)
{
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
else
{
if (abs(ct[icell].v_m[k].z) == 1)
// there is always H+ and OH-...
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
else
{
dum1 = it_sc->Get_mass_water() / t_aq1;
dum2 = it_sc->Get_z_gMCD_map()[1] / dum1;
g_i += pow(dum2, ct[icell].v_m[k].z) * dum1;
}
}
}
g_i *= sol_D[jcell].spec[j].erm_ddl;
}
if (dl_aq2)
{
for (it_sc = s_charge_p2.begin(); it_sc != s_charge_p2.end(); it_sc++)
{
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_j *= sol_D[jcell].spec[j].erm_ddl;
}
}
b_i = A1;
b_j = A2 * sol_D[jcell].spec[j].Dwt;
if (sol_D[icell].tk_x == sol_D[jcell].tk_x)
b_i *= sol_D[jcell].spec[j].Dwt;
else
{
dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
dum2 *= exp(sol_D[jcell].spec[j].dw_t / sol_D[icell].tk_x - sol_D[jcell].spec[j].dw_t / sol_D[jcell].tk_x);
dum2 *= sol_D[icell].viscos_f;
b_i *= dum2;
}
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
k++;
}
if (j < j_max)
j++;
}
else if (strcmp(sol_D[icell].spec[i].name, sol_D[jcell].spec[j].name) == 0)
{
/* species 'name' is in both cells */
if (il_calcs && sol_D[icell].spec[i].type == EX)
{
ct[icell].J_ij_il[k_il].name = sol_D[icell].spec[i].name;
if (sol_D[icell].spec[i].Dwt == 0 || sol_D[jcell].spec[j].Dwt == 0)
ct[icell].v_m_il[k_il].D = 0.0;
else
ct[icell].v_m_il[k_il].D =
(sol_D[icell].spec[i].Dwt + sol_D[jcell].spec[j].Dwt) / 2;
ct[icell].v_m_il[k_il].z = sol_D[icell].spec[i].z;
ct[icell].v_m_il[k_il].Dz = ct[icell].v_m_il[k_il].D * ct[icell].v_m_il[k_il].z;
ct[icell].v_m_il[k_il].Dzc = ct[icell].v_m_il[k_il].Dz * (sol_D[icell].spec[i].c +
sol_D[jcell].spec[j].c) * cec12 / ct[icell].v_m_il[k_il].z;
ct[icell].Dz2c_il += ct[icell].v_m_il[k_il].Dzc * ct[icell].v_m_il[k_il].z;
ct[icell].v_m_il[k_il].grad = (sol_D[jcell].spec[j].c - sol_D[icell].spec[i].c) *
cec12 / ct[icell].v_m_il[k_il].z; /* use equivalent fraction */
k_il++;
}
else
{
ct[icell].J_ij[k].name = sol_D[icell].spec[i].name;
ct[icell].v_m[k].z = sol_D[icell].spec[i].z;
ct[icell].v_m[k].grad = (sol_D[jcell].spec[j].c - sol_D[icell].spec[i].c);
c1 = sol_D[icell].spec[i].c / 2;
c2 = sol_D[jcell].spec[j].c / 2;
ct[icell].v_m[k].c = c1 + c2;
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)
{
// compare diffuse and electromotive forces
dum = ct[icell].v_m[k].grad;
if (icell == 0)
dum2 = (cell_data[1].potV - cell_data[0].potV) / (cell_data[1].length / 2);
else if (icell == count_cells)
dum2 = (cell_data[count_cells + 1].potV - cell_data[count_cells].potV) / (cell_data[count_cells].length / 2);
else
dum2 = (cell_data[jcell].potV - cell_data[icell].potV) / ((cell_data[jcell].length + cell_data[icell].length) / 2);
dum2 *= F_Re3 / tk_x2 * ct[icell].v_m[k].z * (c1 + c2);
// don't transport unavailable moles against the gradient
if (fabs(dum) < fabs(dum2) &&
((dum2 >= 0 && sol_D[jcell].spec[j].c * aq2 < 1e-12) ||
(dum2 <= 0 && sol_D[icell].spec[i].c * aq1 < 1e-12)))
{
// step out:
if (i < i_max)
i++;
if (j < j_max)
j++;
continue;
}
}
g_i = g_j = 0;
if (ct[icell].dl_s > 0)
{
if (dl_aq1)
{
for (it_sc = s_charge_p1.begin(); it_sc != s_charge_p1.end(); it_sc++)
{
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_i *= sol_D[icell].spec[i].erm_ddl;
}
if (dl_aq2)
{
for (it_sc = s_charge_p2.begin(); it_sc != s_charge_p2.end(); it_sc++)
{
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_j *= sol_D[jcell].spec[j].erm_ddl;
}
}
b_i = A1 * sol_D[icell].spec[i].Dwt;
b_j = A2 * sol_D[jcell].spec[j].Dwt;
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
//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;
dum1 = sol_D[jcell].spec[j].lg - sol_D[icell].spec[i].lg;
dum = 1;
if (ddlm)
dum += dum1 / ddlm;
if (dum > 0.2 && dum < 2.25)
{
ct[icell].v_m[k].grad *= dum;
if (implicit)
ct[icell].v_m[k].D = dum;
}
//if (ct[icell].v_m[k].b_ij > max_b)
// max_b = ct[icell].v_m[k].b_ij;
k++;
}
if (i < i_max)
i++;
if (j < j_max)
j++;
}
}
if (!dV_dcell && !ct[icell].Dz2c)
k = 0;
ct[icell].J_ij_count_spec = i_max = k;
ct[icell].J_ij_il_count_spec = k_il;
if (icell == count_cells)
ct[jcell].J_ij_count_spec = ct[icell].J_ij_count_spec;
if (implicit && stagnant < 2)
return(max_b);
/*
* fill in J_ij...
*/
if (dV_dcell)
{
current_cells[icell].ele = current_cells[icell].dif = 0;
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;
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 = 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);
}
// dV_dcell, 2nd pass, current_x has been calculated, and
// voltage was adapted to give equal current in the cells.
dV_dcell2:
//ct[icell].J_ij_sum = 0;
Sum_zM = 0.0;
for (i = 0; i < ct[icell].J_ij_count_spec; i++)
{
if (ct[icell].v_m[i].z)
Sum_zM += ct[icell].v_m[i].b_ij * ct[icell].v_m[i].z * ct[icell].v_m[i].grad;
}
std::map<std::string, class J_ij_save> J_map;
for (i = 0; i < ct[icell].J_ij_count_spec; i++)
{
class J_ij_save J_save;
ct[icell].J_ij[i].tot1 = -ct[icell].v_m[i].grad;
J_save.flux_c = ct[icell].J_ij[i].tot1;
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;
}
J_save.flux_t = ct[icell].J_ij[i].tot1;
if (stagnant)
{
ct[icell].J_ij[i].tot1 *= ct[icell].v_m[i].b_ij * 2 * mixf;
J_save.flux_c *= ct[icell].v_m[i].b_ij * 2 * mixf;
J_save.flux_t *= ct[icell].v_m[i].b_ij * 2 * mixf;
}
else
{
ct[icell].J_ij[i].tot1 *= ct[icell].v_m[i].b_ij * DDt;
J_save.flux_c *= ct[icell].v_m[i].b_ij * DDt;
J_save.flux_t *= ct[icell].v_m[i].b_ij * DDt;
}
J_map[ct[icell].J_ij[i].name] = J_save;
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;
}
cell_J_ij[icell] = J_map;
// assure that icell and jcell have dl water when checking negative conc's in MCD
ct[icell].dl_s = dl_aq1;
ct[jcell].dl_s = dl_aq2;
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...
//current_cells[icell].ele = current_cells[icell].dif = 0;
//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;
// current_cells[icell].dif -= ct[icell].v_m[i].b_ij * ct[icell].v_m[i].z *
// ct[icell].v_m[i].grad;
//}
//dum1 = (current_x - current_cells[icell].dif) / current_cells[icell].ele;
//if (isnan(dum1))
// dum1 = 1;
//dum *= dum1;
for (i = 0; i < ct[icell].J_ij_count_spec; i++)
{
if (!ct[icell].v_m[i].z)
continue;
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;
}
/*
* calculate interlayer mass transfer...
*/
if (il_calcs && ct[icell].Dz2c_il != 0 && ct[icell].J_ij_il_count_spec > 0)
{
cxxExchange *ex_ptr1 = Utilities::Rxn_find(Rxn_exchange_map, icell);
cxxExchange *ex_ptr2 = Utilities::Rxn_find(Rxn_exchange_map, jcell);
Sum_zM = 0.0;
i_max = k_il = ct[icell].J_ij_il_count_spec;
for (i = 0; i < i_max; i++)
Sum_zM += ct[icell].v_m_il[i].Dz * ct[icell].v_m_il[i].grad;
for (i = 0; i < i_max; i++)
{
ct[icell].J_ij_il[i].tot1 = -ct[icell].v_m_il[i].D * ct[icell].v_m_il[i].grad +
Sum_zM * ct[icell].v_m_il[i].Dzc / ct[icell].Dz2c_il;
if (stagnant)
ct[icell].J_ij_il[i].tot1 *= ct[icell].mixf_il;
else
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; // appt not used now
}
/* express the transfer in elemental moles... */
tot1_h = tot1_o = tot2_h = tot2_o = 0.0;
m_s = (class M_S *) free_check_null(m_s);
m_s = (class M_S *) PHRQ_malloc((size_t) count_moles_added *
sizeof(class M_S));
if (m_s == NULL)
malloc_error();
for (i1 = 0; i1 < count_moles_added; 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, icell, stagnant);
/* do the mass transfer... */
if (icell > 0 || stagnant)
{
size_t k;
for (k = 0; k < ex_ptr1->Get_exchange_comps().size(); k++)
{
cxxNameDouble nd(ex_ptr1->Get_exchange_comps()[k].Get_totals());
cxxNameDouble::iterator it = nd.begin();
i_max = 0;
for (; it != nd.end(); it++)
{
if (strcmp("X", it->first.c_str()) == 0)
i_max = 1;
}
if (i_max)
break;
}
if (k < ex_ptr1->Get_exchange_comps().size())
{
cxxExchComp &comp_ref = ex_ptr1->Get_exchange_comps()[k];
cxxNameDouble nd(comp_ref.Get_totals());
cxxNameDouble::iterator it = nd.begin();
/* transfer O and H... */
for (; it != nd.end(); it++)
{
class element *elt_ptr = element_store(it->first.c_str());
LDBLE coef = it->second;
if (strcmp("H", elt_ptr->name) == 0)
{
if (coef < rc1 * tot1_h)
{
tot1_h -= coef;
comp_ref.Get_totals().insert("H", 0);
}
else
{
comp_ref.Get_totals().insert("H", coef - rc1 * tot1_h);
tot1_h *= (1 - rc1);
}
}
else if (strcmp("O", elt_ptr->name) == 0)
{
if (coef < rc1 * tot1_o)
{
tot1_o -= coef;
comp_ref.Get_totals().insert("O", 0);
}
else
{
comp_ref.Get_totals().insert("O", coef - rc1 * tot1_o);
tot1_o *= (1 - rc1);
}
}
}
/* transfer other elements... */
j_max = 0; /* if j_max turns true, reallocate the exchange structure */
for (j = 0; j < count_m_s; j++)
{
// Make sure list includes each element
comp_ref.Get_totals().add(m_s[j].name, 0);
cxxNameDouble nd(comp_ref.Get_totals());
cxxNameDouble::iterator it = nd.begin();
for (; it != nd.end(); it++)
{
class element *elt_ptr = element_store(it->first.c_str());
LDBLE coef = it->second;
if (strcmp(m_s[j].name, elt_ptr->name) != 0)
continue;
/* rc1 part goes to exchange species... */
if (coef < rc1 * m_s[j].tot1)
{
m_s[j].tot1 -= coef;
comp_ref.Get_totals().insert(m_s[j].name, 0);
}
else
{
comp_ref.Get_totals().insert(m_s[j].name, coef - rc1 * m_s[j].tot1);
m_s[j].tot1 *= (1 - rc1);
}
}
}
}
}
if (icell < count_cells || stagnant)
{
size_t k;
for (k = 0; k < ex_ptr2->Get_exchange_comps().size(); k++)
{
cxxExchComp &comp_ref = ex_ptr2->Get_exchange_comps()[k];
cxxNameDouble nd(comp_ref.Get_totals());
cxxNameDouble::iterator it = nd.begin();
i_max = 0;
for (; it != nd.end(); it++)
{
if (strcmp("X", it->first.c_str()) == 0)
i_max = 1;
}
if (i_max)
break;
}
if (k < ex_ptr2->Get_exchange_comps().size())
{
cxxExchComp &comp_ref = ex_ptr2->Get_exchange_comps()[k];
cxxNameDouble nd(comp_ref.Get_totals());
cxxNameDouble::iterator it = nd.begin();
/* transfer O and H... */
for (; it != nd.end(); it++)
{
class element *elt_ptr = element_store(it->first.c_str());
LDBLE coef = it->second;
if (strcmp("H", elt_ptr->name) == 0)
{
if (coef < -rc2 * tot2_h)
{
tot2_h += coef;
comp_ref.Get_totals().insert("H", 0);
}
else
{
comp_ref.Get_totals().insert("H", coef + rc2 * tot2_h);
tot2_h *= (1 - rc2);
}
}
else if (strcmp("O", elt_ptr->name) == 0)
{
if (coef < -rc2 * tot2_o)
{
tot2_o += coef;
comp_ref.Get_totals().insert("O", 0);
}
else
{
comp_ref.Get_totals().insert("O", coef + rc2 * tot2_o);
tot2_o *= (1 - rc2);
}
}
}
/* transfer other elements... */
for (j = 0; j < count_m_s; j++)
{
// Make sure list includes each element
comp_ref.Get_totals().add(m_s[j].name, 0);
cxxNameDouble nd(comp_ref.Get_totals());
cxxNameDouble::iterator it = nd.begin();
for (; it != nd.end(); it++)
{
class element *elt_ptr = element_store(it->first.c_str());
LDBLE coef = it->second;
if (strcmp(m_s[j].name, elt_ptr->name) != 0)
continue;
/* rc2 part goes to exchange species... */
if (coef < -rc2 * m_s[j].tot2)
{
m_s[j].tot2 += coef;
comp_ref.Get_totals().insert(m_s[j].name, 0);
}
else
{
comp_ref.Get_totals().insert(m_s[j].name, coef + rc2 * m_s[j].tot2);
m_s[j].tot2 *= (1 - rc2);
}
}
}
}
}
}
/* do not transport charge imbalance */
//ct[icell].J_ij_sum = 0;
//V_M = (struct V_M *) free_check_null(V_M);
if (il_calcs)
ct[icell].v_m_il = (struct V_M *) free_check_null(ct[icell].v_m_il);
return (il_calcs);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
disp_surf(LDBLE DDt)
/* ---------------------------------------------------------------------- */
/*
* Disperse/diffuse surfaces:
* surface[n1] = SS(mixf * surface[n2]) + (1 - SS(mixf) * surface[i1])
* where SS is Sum of, f2 is a mixing factor.
* Implementation:
* Mobile surface comps and charges are mixed face by face and 1 by 1 in surface[n1]:
Step (from cell 1 to count_cells + 1):
* 1. Surface n1 is made a copy of cell[i1]'s surface, if it exists, or else
* b. a copy of the first encountered mobile surface[n2] from cell i2 = i1 - 1.
* 2 a. Column surfaces are mixed by dispersion/diffusion
* b. Column surfaces are mixed by MCD (if multi_Dflag is true)
* 3. Surfaces n1 and n2 are stored in a temporary surface, with n_user = i1 or i2.
* 4. When done for all cells, new surfaces are copied into the cells.
* NOTE... The surfaces' diffuse_layer, edl and phases/kinetics relations must be identical,
* but only mobile surface_comp's (Dw > 0) and their surface_charge's are transported.
*/
{
int i, i1, i2, k, k1/*, n1, n2*/;
int charge_done, surf1, surf2;
LDBLE f1, f2, mixf, mixf_store, mcd_mixf;
LDBLE lav, A_ij, por, Dp1, Dp2;
cxxMix * mix_ptr;
cxxSurface *surface_ptr1, *surface_ptr2;
LDBLE viscos_f;
/*
* temperature and viscosity correction for MCD coefficient, D_T = D_298 * viscos_298 / viscos
*/
viscos_f = viscos_0;
viscos_f = viscos_0_25 / viscos_f;
//n1 = 0;
//n2 = n1 + 1;
cxxSurface surface_n1, surface_n2;
cxxSurface *surface_n2_ptr;
std::map<int, cxxSurface> Rxn_temp_surface_map;
for (i1 = 1; i1 <= count_cells + 1; i1++)
{
if (i1 <= count_cells && cell_data[i1].por < multi_Dpor_lim)
continue;
if (i1 == 1 && bcon_first != 1)
continue;
if (i1 == count_cells + 1 && bcon_last != 1)
continue;
i2 = i1 - 1;
if (i2 > 0 && cell_data[i2].por < multi_Dpor_lim)
continue;
/*
* step 1. define surface n1 from cell i1, if it exists...
*/
surface_n1.Set_n_user(-99);
surface_n2.Set_n_user(-99);
surface_ptr1 = Utilities::Rxn_find(Rxn_surface_map, i1);
if (surface_ptr1 != NULL)
{
surface_n1 = *surface_ptr1;
}
surface_ptr2 = Utilities::Rxn_find(Rxn_surface_map, i2);
surf1 = surf2 = 0;
if (surface_ptr2 != NULL)
{
if (surface_ptr2->Get_transport())
surf2 = 1;
}
if (surface_ptr1 != NULL)
{
if (surface_ptr1->Get_transport())
surf1 = 1;
}
mixf = mixf_store = 0;
if (surf1 || surf2)
{
/*
* step 2a. Dispersive mixing of mobile surfaces from column cells i2 < i1...
*/
if (i1 <= count_cells)
mix_ptr = Utilities::Rxn_find(Dispersion_mix_map, i1);
else
mix_ptr = Utilities::Rxn_find(Dispersion_mix_map, count_cells);
std::vector<int> num;
std::vector<LDBLE> frac;
mix_ptr->Vectorize(num, frac);
for (size_t i3 = 0; i3 < num.size(); i3++)
{
if (i1 <= count_cells)
{
if (num[i3] == i2)
{
mixf = mixf_store = frac[i3];
break;
}
}
else if (num[i3] == i1)
{
mixf = mixf_store = frac[i3];
break;
}
}
/* step 1b. If cell i1 has no surface, get the mobile comps from surface i2... */
if (surface_n1.Get_n_user() == -99)
{
surface_n1 = mobile_surface_copy(surface_ptr2, i1, false);
/* limit charges to 1... */
if (surface_n1.Get_surface_charges().size() > 1)
{
std::string charge_name = surface_n1.Get_surface_charges()[0].Get_name();
surface_n1 = sum_surface_comp(&surface_n1, 0,
&surface_n1, charge_name, 1, surface_n1.Get_surface_comps()[0].Get_Dw());
}
f1 = 0;
}
else
f1 = 1;
/* find the (possibly modified) surface in the previous cell i2... */
f2 = 1;
if (i2 > 0 || bcon_first == 1)
{
surface_n2_ptr = Utilities::Rxn_find(Rxn_temp_surface_map, i2);
if (surface_n2_ptr != NULL)
{
surface_n2 = *surface_n2_ptr;
}
/* if not found... */
else
{
/* copy it from surface_ptr2... */
if (surface_ptr2 != NULL)
{
surface_n2 = *surface_ptr2;
surface_n2.Set_n_user_both(i2);
}
else
{
/* or make it a mobile copy of the surface in cell i1... */
surface_n2 = mobile_surface_copy(surface_ptr1, i2, false);
/* limit charges to 1... */
if (surface_n2.Get_surface_charges().size() > 1)
{
std::string charge_name = surface_n2.Get_surface_charges()[0].Get_name();
surface_n2 = sum_surface_comp(&surface_n2, 0, &surface_n2, charge_name, 1,
surface_n2.Get_surface_comps()[0].Get_Dw());
}
f2 = 0;
}
}
}
/*
* For MCD, step 2b. Find physical dimensions and porosity of the cells...
*/
if (i1 == 1)
{
por = cell_data[1].por;
lav = cell_data[1].length / 2;
A_ij = Utilities::Rxn_find(Rxn_solution_map, 1)->Get_mass_water() / cell_data[1].length;
}
else if (i1 == count_cells + 1)
{
por = cell_data[count_cells].por;
lav = cell_data[count_cells].length / 2;
A_ij =
Utilities::Rxn_find(Rxn_solution_map, count_cells)->Get_mass_water() /
cell_data[count_cells].length;
}
else
{
por = cell_data[i2].por;
lav =
(cell_data[i1].length + cell_data[i2].length) / 2;
A_ij =
Utilities::Rxn_find(Rxn_solution_map, i1)->Get_mass_water() /
(cell_data[i1].length * cell_data[i1].por);
A_ij +=
Utilities::Rxn_find(Rxn_solution_map, i2)->Get_mass_water() /
(cell_data[i2].length * cell_data[i2].por);
A_ij /= 2;
A_ij *=
(cell_data[i1].por <
cell_data[i2].por ? cell_data[i1].por : cell_data[i2].por);
}
/* mix in comps with the same charge structure... */
if (surf2)
{
for (k = 0; k < (int) surface_ptr2->Get_surface_comps().size(); k++)
{
cxxSurfaceComp *comp_k_ptr = &(surface_ptr2->Get_surface_comps()[k]);
std::string charge_name = comp_k_ptr->Get_charge_name();
if (comp_k_ptr->Get_Dw() == 0)
continue;
charge_done = FALSE;
for (k1 = 0; k1 < k; k1++)
{
cxxSurfaceComp *comp_k1_ptr = &(surface_ptr2->Get_surface_comps()[k1]);
if (comp_k_ptr->Get_charge_name() ==
comp_k1_ptr->Get_charge_name())
{
charge_done = TRUE;
break;
}
}
if (charge_done)
continue;
/* Define the MCD diffusion coefficient... */
mcd_mixf = 0;
if (multi_Dflag)
{
Dp2 = comp_k_ptr->Get_Dw() * pow(por, multi_Dn);
Dp1 = 0;
if (surface_ptr1 != NULL && surface_ptr1->Get_transport())
{
for (k1 = 0; k1 < (int) surface_ptr1->Get_surface_comps().size(); k1++)
{
cxxSurfaceComp *comp_k1_ptr = &(surface_ptr1->Get_surface_comps()[k1]);
if (strcmp(comp_k1_ptr->Get_formula().c_str(),
comp_k_ptr->Get_formula().c_str()) != 0)
continue;
Dp1 =
comp_k1_ptr->Get_Dw() *
pow(cell_data[i1].por, multi_Dn);
break;
}
}
if (Dp1 > 0)
Dp2 = (Dp2 + Dp1) / 2;
/* and the mixing factor... */
mcd_mixf = Dp2 * viscos_f * A_ij * DDt / lav;
}
mixf = mixf_store + mcd_mixf;
if (mixf < 1e-8)
mixf = 0;
if (mixf > 0.99999999)
mixf = 0.99999999;
if (i1 <= count_cells)
{
surface_n1 = sum_surface_comp(&surface_n1, f1, surface_ptr2, charge_name, mixf,
surface_ptr2->Get_surface_comps()[k].Get_Dw());
f1 = 1;
}
if (i2 > 0)
{
surface_n2 = sum_surface_comp(&surface_n2, f2, surface_ptr2, charge_name, -mixf,
surface_ptr2->Get_surface_comps()[k].Get_Dw());
f2 = 1;
}
surface_n1.Set_n_user_both(i1);
}
}
if (surf1)
{
for (k = 0; k < (int) surface_ptr1->Get_surface_comps().size(); k++)
{
cxxSurfaceComp * comp_k_ptr = &(surface_ptr1->Get_surface_comps()[k]);
std::string charge_name = comp_k_ptr->Get_charge_name();
if (comp_k_ptr->Get_Dw() == 0)
continue;
charge_done = FALSE;
for (k1 = 0; k1 < k; k1++)
{
cxxSurfaceComp * comp_k1_ptr = &(surface_ptr1->Get_surface_comps()[k1]);
if (comp_k_ptr->Get_charge_name() ==
comp_k1_ptr->Get_charge_name())
{
charge_done = TRUE;
break;
}
}
if (charge_done)
continue;
/* Define the MCD diffusion coefficient... */
mcd_mixf = 0;
if (multi_Dflag)
{
if (i1 <= count_cells)
{
Dp1 =
comp_k_ptr->Get_Dw() *
pow(cell_data[i1].por, multi_Dn);
}
else
{
Dp1 = comp_k_ptr->Get_Dw() * pow(por, multi_Dn);
}
Dp2 = 0;
if (surface_ptr2 != NULL && surface_ptr2->Get_transport())
{
for (k1 = 0; k1 < (int) surface_ptr2->Get_surface_comps().size(); k1++)
{
cxxSurfaceComp * comp_k1_ptr = &(surface_ptr2->Get_surface_comps()[k1]);
if (strcmp(comp_k1_ptr->Get_formula().c_str(),
comp_k_ptr->Get_formula().c_str()) != 0)
continue;
Dp2 = comp_k1_ptr->Get_Dw() * pow(por, multi_Dn);
break;
}
}
if (Dp2 > 0)
Dp1 = (Dp1 + Dp2) / 2;
/* and the mixing factor... */
mcd_mixf = Dp1 * viscos_f * A_ij * DDt / lav;
}
mixf = mixf_store + mcd_mixf;
if (mixf < 1e-8)
mixf = 0;
if (mixf > 0.99999999)
mixf = 0.99999999;
if (i2 > 0)
{
surface_n2 = sum_surface_comp
(&surface_n2, f2, surface_ptr1, charge_name, mixf,
surface_ptr1->Get_surface_comps()[k].Get_Dw());
f2 = 1;
}
if (i1 <= count_cells)
{
surface_n1 = sum_surface_comp
(&surface_n1, f1, surface_ptr1, charge_name, -mixf,
surface_ptr1->Get_surface_comps()[k].Get_Dw());
f1 = 1;
}
surface_n2.Set_n_user_both(i2);
}
}
}
/*
* Step 3. copy surface[n1] and [n2] in a new temporary surface...
*/
if (surface_n1.Get_n_user() == -99)
continue;
Rxn_temp_surface_map[i1] = surface_n1;
{
cxxSurface t;
surface_n1 = t;
}
surface_n1.Set_n_user_both(-99);
if (surface_n2.Get_n_user() == i2)
{
surface_n2.Set_n_user_both(i2);
Rxn_temp_surface_map[i2] = surface_n2;
{
cxxSurface t;
surface_n2 = t;
}
surface_n2.Set_n_user_both(-99);
}
}
/*
* Step 4. Dispersion/diffusion is done. New surfaces can be copied in the cell's surface...
*/
//n2 = 0;
std::map<int, cxxSurface>::iterator jit = Rxn_temp_surface_map.begin();
for (; jit != Rxn_temp_surface_map.end(); jit++)
{
i = jit->first;
assert(i == jit->second.Get_n_user());
if ((i == 0 && bcon_first == 1) || (i == count_cells + 1 && bcon_last == 1))
{
continue;
}
if (i >= 0 && i <= 1 + count_cells * (1 + stag_data.count_stag))
{
surface_ptr1 = Utilities::Rxn_find(Rxn_surface_map, i);
if (surface_ptr1 != NULL)
{
Rxn_surface_map[i] = jit->second;
}
else
{
//Add to map
Rxn_surface_map[i] = jit->second;
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
cxxSurface Phreeqc::
sum_surface_comp(cxxSurface *source1, LDBLE f1, cxxSurface *source2,
std::string charge_name, LDBLE f2, LDBLE new_Dw)
/* ---------------------------------------------------------------------- */
{
/*
* Takes fraction f1 of the 1st surface, adds fraction f2 of the 2nd surface's comps[k] and its charge.
* The result goes in target
*/
int new_n_user;
cxxSurface *surface_ptr1, *surface_ptr2;
std::string token;
/*
* Find surfaces
*/
surface_ptr1 = source1;
if (surface_ptr1 == NULL)
{
error_string = sformatf("Null pointer for surface 1 in sum_surface.");
error_msg(error_string, STOP);
input_error++;
return (ERROR);
}
surface_ptr2 = source2;
/*
* Store data for structure surface
*/
new_n_user = surface_ptr1->Get_n_user();
cxxSurface temp_surface(*surface_ptr1);
temp_surface.Set_n_user_both(new_n_user);
temp_surface.Set_description("Copy");
temp_surface.Set_solution_equilibria(false);
temp_surface.Set_n_solution(-99);
/*
* Multiply component compositions by f1
*/
temp_surface.multiply(f1);
/*
* Add in surface_ptr2
*/
// Only components with same charge as component k
cxxSurface addee(*surface_ptr2);
addee.Get_surface_comps().clear();
addee.Get_surface_charges().clear();
for (std::vector<cxxSurfaceComp>::iterator it = surface_ptr2->Get_surface_comps().begin();
it != surface_ptr2->Get_surface_comps().end(); it++)
{
if (it->Get_charge_name() == charge_name)
{
addee.Get_surface_comps().push_back(*it);
}
}
for (std::vector<cxxSurfaceCharge>::iterator it = surface_ptr2->Get_surface_charges().begin();
it != surface_ptr2->Get_surface_charges().end(); it++)
{
if (it->Get_name() == charge_name)
{
addee.Get_surface_charges().push_back(*it);
}
}
if (f2 == 0)
f2 = 1e-30;
temp_surface.add(addee, f2);
temp_surface.Set_transport(false);
for (size_t i = 0; i < temp_surface.Get_surface_comps().size(); i++)
{
if (temp_surface.Get_surface_comps()[i].Get_charge_name() == charge_name)
{
temp_surface.Get_surface_comps()[i].Set_Dw(new_Dw);
}
if (temp_surface.Get_surface_comps()[i].Get_Dw() > 0)
{
temp_surface.Set_transport(true);
}
}
/*
* Finish up
*/
temp_surface.Sort_comps();
return (temp_surface);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
check_surfaces(cxxSurface *surface_ptr1, cxxSurface *surface_ptr2)
/* ---------------------------------------------------------------------- */
{
/* checks if surfaces can be mixed...
*/
int n_user1, n_user2, return_code;
return_code = OK;
n_user1 = surface_ptr1->Get_n_user();
n_user2 = surface_ptr2->Get_n_user();
if (surface_ptr1->Get_dl_type() != surface_ptr2->Get_dl_type())
{
error_string = sformatf(
"Surfaces %d and %d differ in definition of diffuse layer. Cannot mix.",
n_user1, n_user2);
error_msg(error_string, STOP);
return_code = ERROR;
input_error++;
}
if (surface_ptr1->Get_type() != surface_ptr2->Get_type())
{
error_string = sformatf(
"Surfaces %d and %d differ in use of electrical double layer. Cannot mix.",
n_user1, n_user2);
error_msg(error_string, STOP);
return_code = ERROR;
input_error++;
}
if (surface_ptr1->Get_only_counter_ions() != surface_ptr2->Get_only_counter_ions())
{
error_string = sformatf(
"Surfaces %d and %d differ in use of only counter ions in the diffuse layer. Cannot mix.",
n_user1, n_user2);
error_msg(error_string, STOP);
return_code = ERROR;
input_error++;
}
if (surface_ptr1->Get_related_phases() != surface_ptr2->Get_related_phases())
{
error_string = sformatf(
"Surfaces %d and %d differ in use of related phases (sites proportional to moles of an equilibrium phase). Cannot mix.",
n_user1, n_user2);
error_msg(error_string, STOP);
return_code = ERROR;
input_error++;
}
if (surface_ptr1->Get_related_rate() != surface_ptr2->Get_related_rate())
{
error_string = sformatf(
"Surfaces %d and %d differ in use of related rate (sites proportional to moles of a kinetic reactant). Cannot mix.",
n_user1, n_user2);
error_msg(error_string, STOP);
return_code = ERROR;
input_error++;
}
return (return_code);
}
/* ---------------------------------------------------------------------- */
cxxSurface Phreeqc::
mobile_surface_copy(cxxSurface *surface_old_ptr,
int n_user_new, bool move_old)
/* ---------------------------------------------------------------------- */
{
/*
* Copies mobile comps from surface_old_ptr to surf_ptr1,
* comps and charges with Dw > 0 are moved if move_old == TRUE, else copied.
* NOTE... when all comps are moved, the old surface is deleted and surfaces are sorted again,
* which will modify pointers and surface numbers.
* User number of new surface structure is n_user_new, structure is freed when n_user_new is already defined
*/
cxxSurface temp_surface(*surface_old_ptr);
/*
* Store moving surface's properties in temp_surface
*/
temp_surface.Set_n_user_both(n_user_new);
std::ostringstream desc;
desc << "Surface defined in simulation " << simulation << ".";
temp_surface.Set_description(desc.str().c_str());
temp_surface.Set_solution_equilibria(false);
temp_surface.Set_transport(true);
size_t count_comps = surface_old_ptr->Get_surface_comps().size();
int i1, i2;
i1 = i2 = 0;
temp_surface.Get_surface_comps().clear();
temp_surface.Get_surface_charges().clear();
/* see if comps must be moved, Dw > 0 */
for (size_t i = 0; i < count_comps; i++)
{
cxxSurfaceComp &comp_ref = surface_old_ptr->Get_surface_comps()[i];
if (comp_ref.Get_Dw() > 0)
{
i1++;
// copy comp
temp_surface.Get_surface_comps().push_back(comp_ref);
// copy charge, if needed
cxxSurfaceCharge *charge_ptr = temp_surface.Find_charge(comp_ref.Get_charge_name());
if (charge_ptr == NULL)
{
i2++;
cxxSurfaceCharge *old_charge_ptr = surface_old_ptr->Find_charge(comp_ref.Get_charge_name());
temp_surface.Get_surface_charges().push_back(*old_charge_ptr);
}
}
}
if (i1 > 0)
{
/* OK, store moved parts from old surface, but first:
get immobile surface comps from new surface... */
cxxSurface *surf_ptr = Utilities::Rxn_find(Rxn_surface_map, n_user_new);
if (surf_ptr != NULL)
{
for (size_t k = 0; k < surf_ptr->Get_surface_comps().size(); k++)
{
cxxSurfaceComp &comp_ref = surf_ptr->Get_surface_comps()[k];
if (comp_ref.Get_Dw() > 0)
continue;
bool charge_done(false);
for (size_t k1 = 0; k1 < k; k1++)
{
if (surf_ptr->Get_surface_comps()[k1].Get_charge_name() ==
comp_ref.Get_charge_name())
{
charge_done = true;
break;
}
}
if (charge_done)
continue;
temp_surface = sum_surface_comp(&temp_surface, 1, surf_ptr, comp_ref.Get_charge_name(), 1, 0);
}
}
// Return value is temp_surface
temp_surface.Set_n_user_both(n_user_new);
}
/* delete moved parts from old surface */
if (move_old && i1 > 0)
{
cxxSurface replace_old(temp_surface);
int n_user_old = surface_old_ptr->Get_n_user();
if ((size_t) i1 != count_comps)
{
/* redefine old surface with only unmovable comps (Dw = 0) */
/* other temp_surface members were set above */
replace_old.Set_n_user_both(n_user_old);
replace_old.Set_transport(false);
replace_old.Get_surface_comps().clear();
replace_old.Get_surface_charges().clear();
i1 = i2 = 0;
for (size_t i = 0; i < count_comps; i++)
{
if (surface_old_ptr->Get_surface_comps()[i].Get_Dw() == 0)
{
i1++;
// copy surface comp
cxxSurfaceComp & comp_ref = surface_old_ptr->Get_surface_comps()[i];
replace_old.Get_surface_comps().push_back(comp_ref);
cxxSurfaceCharge *charge_ptr = replace_old.Find_charge(comp_ref.Get_charge_name());
// copy surface charge if necessary
if (charge_ptr == NULL)
{
i2++;
cxxSurfaceCharge *old_charge_ptr = surface_old_ptr->Find_charge(comp_ref.Get_charge_name());
replace_old.Get_surface_charges().push_back(*old_charge_ptr);
}
}
}
if (replace_old.Get_surface_comps().size() == 0)
{
Rxn_surface_map.erase(surface_old_ptr->Get_n_user());
}
else
{
replace_old.Sort_comps();
Rxn_surface_map[surface_old_ptr->Get_n_user()] = replace_old;
}
}
}
temp_surface.Sort_comps();
return (temp_surface);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
diff_stag_surf(int mobile_cell)
/* ---------------------------------------------------------------------- */
/*
* Diffuse stagnant and mobile surfaces, following the steps of disp_surf.
* First the mobile/stagnant surfaces are mixed, then the stagnant surfaces
* when not already done.
* If mixing factors among the cells are defined expicitly, it is assumed that
* mixing with a lower numbered cell was done when that cell was processed:
* for any cell in MCD, need only include the mixing factors for higher numbered cells.
*/
{
int i, i1, i2, k, k1, /*n1, n2,*/ ns;
int charge_done, surf1, surf2;
LDBLE f1, f2, mixf, mixf_store;
LDBLE Dp1, Dp2;
cxxMix *mix_ptr;
cxxSurface *surface_ptr1, *surface_ptr2;
LDBLE viscos_f;
/*
* temperature and viscosity correction for MCD coefficient, D_T = D_298 * Tk * viscos_298 / (298 * viscos)
*/
viscos_f = viscos_0;
viscos_f = viscos_0_25 / viscos_f;
cxxSurface surface_n1, surface_n2;
cxxSurface *surface_n1_ptr = &surface_n1;
cxxSurface *surface_n2_ptr;
std::map<int, cxxSurface> Rxn_temp_surface_map;
for (ns = 0; ns < stag_data.count_stag; ns++)
{
i1 = mobile_cell + 1 + ns * count_cells;
if (ns == 0)
i1--;
if (cell_data[i1].por < multi_Dpor_lim)
continue;
surface_n1.Set_n_user_both(-99);
surface_n2.Set_n_user_both(-99);
surface_ptr1 = Utilities::Rxn_find(Rxn_surface_map, i1);
/*
* step 2a. mix surfaces...
*/
mix_ptr = Utilities::Rxn_find(Rxn_mix_map, i1);
if (mix_ptr == NULL)
continue;
std::vector<int> num;
std::vector<LDBLE> frac;
mix_ptr->Vectorize(num, frac);
for (size_t i3 = 0; i3 < num.size(); i3++)
{
if ((i2 = num[i3]) <= i1)
continue;
if (cell_data[i2].por < multi_Dpor_lim)
continue;
surface_ptr2 = Utilities::Rxn_find(Rxn_surface_map, i2);
surf1 = surf2 = 0;
if (surface_ptr2 != NULL)
{
if (surface_ptr2->Get_transport())
surf2 = 1;
}
if (surface_ptr1 != NULL)
{
if (surface_ptr1->Get_transport())
surf1 = 1;
}
if (!surf1 && !surf2)
continue;
mixf = mixf_store = frac[i3];;
/* find the (possibly modified) surface in cell i1... */
f1 = 1;
surface_n1_ptr = Utilities::Rxn_find(Rxn_temp_surface_map, i1);
if (surface_n1_ptr != NULL)
{
surface_n1 = *surface_n1_ptr;
//n1 = 0;
}
/* if not found... */
else
{
/* copy it from surface_ptr1... */
if (surface_ptr1 != NULL)
{
surface_n1 = *surface_ptr1;
}
else
{
/* or make it a mobile copy of the surface in cell i2... */
surface_n1 = mobile_surface_copy(surface_ptr2, i1, false);
/* limit comps to 1... */
if (surface_n1.Get_surface_charges().size() > 1)
{
std::string charge_name = surface_n1.Get_surface_charges()[0].Get_name();
surface_n1 = sum_surface_comp(&surface_n1, 0, &surface_n1, charge_name, 1,
surface_n1.Get_surface_comps()[0].Get_Dw());
}
f1 = 0;
}
}
/* find the (possibly modified) surface in cell i2... */
f2 = 1;
surface_n2_ptr = Utilities::Rxn_find(Rxn_temp_surface_map, i2);
if (surface_n2_ptr != NULL)
{
surface_n2 = *surface_n2_ptr;
}
/* if not found... */
else
{
//n2 = 1;
/* copy it from surface_ptr2... */
if (surface_ptr2 != NULL)
{
surface_n2 = *surface_ptr2;
}
else
{
/* or make it a mobile copy of the surface in cell i1... */
surface_n2 = mobile_surface_copy(surface_ptr1, i2, false);
/* limit comps to 1... */
if (surface_n2.Get_surface_charges().size() > 1)
{
std::string charge_name = surface_n2.Get_surface_charges()[0].Get_name();
surface_n2 = sum_surface_comp(&surface_n2, 0, &surface_n2, charge_name, 1,
surface_n2.Get_surface_comps()[0].Get_Dw());
}
f2 = 0;
}
}
/* For MCD, step 2b. Adapt mixf to default values... */
if (multi_Dflag)
{
mixf_store *=
(cell_data[i1].por <= cell_data[i2].por ? cell_data[i1].por : cell_data[i2].por);
mixf_store /= (default_Dw * pow(multi_Dpor, multi_Dn) * multi_Dpor);
}
/* mix in comps with the same charge structure... */
if (surf2)
{
for (k = 0; k < (int) surface_ptr2->Get_surface_comps().size(); k++)
{
cxxSurfaceComp *comp_k_ptr = &(surface_ptr2->Get_surface_comps()[k]);
std::string charge_name = comp_k_ptr->Get_charge_name();
if (comp_k_ptr->Get_Dw() == 0)
continue;
charge_done = FALSE;
for (k1 = 0; k1 < k; k1++)
{
cxxSurfaceComp *comp_k1_ptr = &(surface_ptr2->Get_surface_comps()[k1]);
if (comp_k_ptr->Get_charge_name() == comp_k1_ptr->Get_charge_name())
{
charge_done = TRUE;
break;
}
}
if (charge_done)
continue;
/* find diffusion coefficients of surfaces... */
if (multi_Dflag)
{
Dp2 = comp_k_ptr->Get_Dw() *
pow(cell_data[i2].por, multi_Dn) * viscos_f;
Dp1 = 0;
if (surf1)
{
for (k1 = 0; k1 < (int) surface_ptr1->Get_surface_comps().size(); k1++)
{
cxxSurfaceComp *comp_k1_ptr = &(surface_ptr1->Get_surface_comps()[k1]);
if (strcmp(comp_k1_ptr->Get_formula().c_str(),
comp_k_ptr->Get_formula().c_str()) != 0)
continue;
Dp1 = comp_k1_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f;
break;
}
}
if (Dp1 > 0)
Dp2 = (Dp2 + Dp1) / 2;
/* and adapt the mixing factor... */
mixf = mixf_store * Dp2;
mixf /= Utilities::Rxn_find(Rxn_solution_map, i2)->Get_mass_water();
}
if (mixf < 1e-8)
mixf = 0;
if (mixf > 0.99999999)
mixf = 0.99999999;
surface_n1 = sum_surface_comp(&surface_n1, f1, surface_ptr2, charge_name, mixf,
surface_ptr2->Get_surface_comps()[k].Get_Dw());
f1 = 1;
surface_n2 = sum_surface_comp(&surface_n2, f2, surface_ptr2, charge_name, -mixf,
surface_ptr2->Get_surface_comps()[k].Get_Dw());
}
}
if (surf1)
{
for (k = 0; k < (int) surface_ptr1->Get_surface_comps().size(); k++)
{
cxxSurfaceComp *comp_k_ptr = &(surface_ptr1->Get_surface_comps()[k]);
std::string charge_name = comp_k_ptr->Get_charge_name();
if (comp_k_ptr->Get_Dw() == 0)
continue;
charge_done = FALSE;
for (k1 = 0; k1 < k; k1++)
{
cxxSurfaceComp *comp_k1_ptr = &(surface_ptr1->Get_surface_comps()[k1]);
if (comp_k_ptr->Get_charge_name() == comp_k1_ptr->Get_charge_name())
{
charge_done = TRUE;
break;
}
}
if (charge_done)
continue;
/* find diffusion coefficients of surfaces... */
if (multi_Dflag)
{
Dp1 = comp_k_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f;
Dp2 = 0;
if (surf2)
{
for (k1 = 0; k1 < (int) surface_ptr2->Get_surface_comps().size(); k1++)
{
cxxSurfaceComp *comp_k1_ptr = &(surface_ptr2->Get_surface_comps()[k1]);
if (strcmp(comp_k1_ptr->Get_formula().c_str(),
comp_k_ptr->Get_formula().c_str()) != 0)
continue;
Dp2 = comp_k1_ptr->Get_Dw() * pow(cell_data[i2].por, multi_Dn) * viscos_f;
break;
}
}
if (Dp2 > 0)
Dp1 = (Dp1 + Dp2) / 2;
/* and adapt the mixing factor... */
mixf = mixf_store * Dp1;
mixf /= Utilities::Rxn_find(Rxn_solution_map, i1)->Get_mass_water();
}
if (mixf < 1e-8)
mixf = 0;
if (mixf > 0.99999999)
mixf = 0.99999999;
surface_n2 = sum_surface_comp(&surface_n2, f2, surface_ptr1, charge_name, mixf,
surface_ptr1->Get_surface_comps()[k].Get_Dw());
f2 = 1;
surface_n1 = sum_surface_comp(&surface_n1, f1, surface_ptr1, charge_name, -mixf,
surface_ptr1->Get_surface_comps()[k].Get_Dw());
}
}
/*
* Step 3. copy surface[n1] and [n2] in a new temporary surface...
*/
if (surface_n1.Get_n_user() == -99)
continue;
surface_n1.Set_n_user_both(i1);
Rxn_temp_surface_map[i1] = surface_n1;
assert(surface_n2.Get_n_user() != -99);
assert(surface_n2.Get_n_user() == i2);
surface_n2.Set_n_user_both(i2);
Rxn_temp_surface_map[i2] = surface_n2;
}
}
/*
* Step 4. Diffusion is done. New surfaces can be copied in the cells...
*/
//n2 = 0;
std::map<int, cxxSurface>::iterator jit = Rxn_temp_surface_map.begin();
for (; jit != Rxn_temp_surface_map.end(); jit++)
{
i = jit->first;
assert(i == jit->second.Get_n_user());
if ((i == 0 && bcon_first == 1) || (i == count_cells + 1 && bcon_last == 1))
{
continue;
}
if (i >= 0 && i <= 1 + count_cells * (1 + stag_data.count_stag))
{
surface_ptr1 = Utilities::Rxn_find(Rxn_surface_map, i);
//if (surface_ptr1 != NULL)
//{
// Rxn_surface_map[i] = jit->second;
//}
//else
{
//Add to map
Rxn_surface_map[i] = jit->second;
}
}
}
return (OK);
}
/* ---------------------------------------------------------------------- */
int Phreeqc::
reformat_surf(const char *comp_name, LDBLE fraction, const char *new_comp_name,
LDBLE new_Dw, int l_cell)
/* ---------------------------------------------------------------------- */
{
cxxSurface *surface_ptr;
if ((surface_ptr = Utilities::Rxn_find(Rxn_surface_map, l_cell)) == NULL)
return (OK);
if (surface_ptr->Find_charge(comp_name) == NULL)
return (OK);
// Assume comp_name is charge name
std::string old_charge_name = comp_name;
std::string new_charge_name = new_comp_name;
if (fraction > 0.99999999)
fraction = 0.99999999;
cxxSurface temp_surface(*surface_ptr);
cxxSurface change;
for (size_t i = 0; i < temp_surface.Get_surface_comps().size(); i++)
{
cxxSurfaceComp *comp_ptr = &(temp_surface.Get_surface_comps()[i]);
std::string charge_name = comp_ptr->Get_charge_name();
if (charge_name == old_charge_name)
{
cxxSurfaceComp comp(*comp_ptr);
comp.multiply(fraction);
std::string std_comp_name = comp_ptr->Get_formula();
Utilities::replace(comp_name, new_comp_name, std_comp_name);
comp.Set_formula(std_comp_name.c_str());
comp.Set_charge_name(new_comp_name);
cxxNameDouble nd;
cxxNameDouble::iterator it;
for (it = comp.Get_totals().begin(); it != comp.Get_totals().end(); it++)
{
std::string tot_name(it->first);
Utilities::replace(comp_name, new_comp_name, tot_name);
nd[tot_name] = it->second;
}
comp.Set_totals(nd);
change.Get_surface_comps().push_back(comp);
comp_ptr->multiply(1 - fraction);
}
}
for (size_t i = 0; i < temp_surface.Get_surface_charges().size(); i++)
{
cxxSurfaceCharge *charge_ptr = &(temp_surface.Get_surface_charges()[i]);
std::string charge_name = charge_ptr->Get_name();
if (charge_name == old_charge_name)
{
cxxSurfaceCharge charge(*charge_ptr);
charge.multiply(fraction);
std::string std_charge_name = charge_ptr->Get_name();
Utilities::replace(comp_name, new_comp_name, std_charge_name);
charge.Set_name(std_charge_name.c_str());
change.Get_surface_charges().push_back(charge);
charge_ptr->multiply(1 - fraction);
}
}
temp_surface.add(change, 1.0);
// change Dw
for (size_t i = 0; i < temp_surface.Get_surface_comps().size(); i++)
{
cxxSurfaceComp *comp_ptr = &(temp_surface.Get_surface_comps()[i]);
std::string charge_name = comp_ptr->Get_charge_name();
if (charge_name == new_charge_name)
{
comp_ptr->Set_Dw(new_Dw);
}
}
temp_surface.Set_transport(false);
for (size_t i = 0; i < temp_surface.Get_surface_comps().size(); i++)
{
if (temp_surface.Get_surface_comps()[i].Get_Dw() > 0)
{
temp_surface.Set_transport(true);
break;
}
}
temp_surface.Sort_comps();
Rxn_surface_map[l_cell] = temp_surface;
return OK;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
viscosity(void)
/* ---------------------------------------------------------------------- */
{
/* from Atkins, 1994. Physical Chemistry, 5th ed. */
//viscos =
// pow((LDBLE) 10.,
// -(1.37023 * (tc_x - 20) +
// 0.000836 * (tc_x - 20) * (tc_x - 20)) / (109 + tc_x));
/* Huber et al., 2009, J. Phys. Chem. Ref. Data, Vol. 38, 101-125 */
LDBLE H[4] = { 1.67752, 2.20462, 0.6366564, -0.241605 };
LDBLE Tb = tk_x / 647.096, denom = H[0], mu0;
int i, j, i1;
for (i = 1; i < 4; i++)
denom += H[i] / pow(Tb, i);
mu0 = 100.0 * sqrt(Tb) / denom;
LDBLE H2[42] =
{ 5.20094e-1, 2.22531e-1, -2.81378e-1, 1.61913e-1, -3.25372e-2, 0, 0,
8.50895e-2, 9.99115e-1, -9.06851e-1, 2.57399e-1, 0, 0, 0,
-1.08374, 1.88797, -7.72479e-1, 0, 0, 0, 0,
-2.89555e-1, 1.26613, -4.89837e-1, 0, 6.98452e-2, 0, -4.35673e-3,
0, 0, -2.57040e-1, 0, 0, 8.72102e-3, 0,
0, 1.20573e-1, 0, 0, 0, 0, -5.93264e-4 };
LDBLE Rb = rho_0 / 0.322, Tb_1, Rb_1, S1 = 0.0, S2, mu1;
Tb_1 = 1.0 / Tb - 1.0; Rb_1 = Rb - 1.0;
for (i = 0; i < 6; i++)
{
S2 = 0.0;
for (j = 0; j < 7; j++)
{
i1 = 7 * i;
if (!H2[i1 + j])
continue;
if (j)
S2 += H2[i1 + j] * pow(Rb_1, j);
else
S2 += H2[i1];
}
if (i)
S1 += S2 * pow(Tb_1, i);
else
S1 += S2;
}
mu1 = exp(Rb * S1);
viscos_0 = viscos = mu0 * mu1 / 1e3;
viscos_0_25 = 0.8900239182946;
//#define OLD_VISCOSITY
#ifdef OLD_VISCOSITY
/* from Atkins, 1994. Physical Chemistry, 5th ed. */
viscos =
pow((LDBLE) 10.,
-(1.37023 * (tc_x - 20) +
0.000836 * (tc_x - 20) * (tc_x - 20)) / (109 + tc_x));
viscos_0_25 = 0.88862;
#endif
//return viscos;
if (!print_viscosity)
return viscos;
/* (modified) Jones-Dole eqn for viscosity:
viscos / viscos_0 =
1 + A * eq_tot^0.5 +
f_an * (Sum(B_i * m_i) +
Sum(D_i * m_i * ((1 + f_I) * mu_x^d3_i + (m_i * f_z)^d3_i) / (2 + f_I)))
A calculated from Falkenhagen-Dole
B_i = b0 + b1*exp(b2 * tc), b0..2 in Jones_Dole[0..2], read in SOLUTION_SPECIES
D_i = d1 * exp(d2 * tc), d1, 2 in Jones_Dole[3, 4]
d3_i in Jones_Dole[5]
Jones_Dole[6] contains the anion factor, 1 for Cl-, variable for other anions
f_z = (z * z + |z|) / 2, the contribution of the ion to mu_x, if z = 0: f_z = mu_x / m_i
f_I = variable, depends on d3_i > 1, or d3_i < 1.
tc is limited to 200 C.
A from Falkenhagen-Dole for a salt:
A = 4.3787e-14 * TK**1.5 / (eps_r**0.5)* (z1 + z2)**-0.5 / (D1 * D2) * psi
psi = (D1*z2 + D2*z1)/4 - z1*z2 * (D1-D2)**2 / ((D1*z1 + D2*z2)**0.5 + ((D1 + D2) * (z1 + z2))**0.5)**2
D1, z1 for the cation, D2, |z2| for the anion of the salt.
We use the harmonic mean of the Dw's, and the arithmetic mean of the z's,
both weighted by the equivalent concentration.
*/
LDBLE D1, D2, z1, z2, m_plus, m_min, eq_plus, eq_min, eq_dw_plus, eq_dw_min, t1, t2, t3, fan = 1;
LDBLE A, psi, Bc = 0, Dc = 0, Dw = 0.0, l_z, f_z, lm, V_an, m_an, V_Cl, tc;
m_plus = m_min = eq_plus = eq_min = eq_dw_plus = eq_dw_min = V_an = m_an = V_Cl = 0;
tc = (tc_x > 200) ? 200 : tc_x;
for (i = 0; i < (int)this->s_x.size(); i++)
{
if (s_x[i]->type != AQ && s_x[i]->type > HPLUS)
continue;
if ((lm = s_x[i]->lm) < -9)
continue;
if (s_x[i]->Jones_Dole[0] || s_x[i]->Jones_Dole[1] || s_x[i]->Jones_Dole[3])
{
s_x[i]->dw_t_visc = 0;
t1 = s_x[i]->moles / mass_water_aq_x;
l_z = fabs(s_x[i]->z);
if (l_z)
f_z = (l_z * l_z + l_z) / 2;
else
f_z = mu_x / t1;
//if data at tc's other than 25 are scarce, put the values found for 25 C in [7] and [8], optimize [1], [2], and [4]...
if (s_x[i]->Jones_Dole[7] || s_x[i]->Jones_Dole[8])
{
s_x[i]->Jones_Dole[0] = s_x[i]->Jones_Dole[7] -
s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * 25.0);
s_x[i]->Jones_Dole[3] =
s_x[i]->Jones_Dole[8] / exp(-s_x[i]->Jones_Dole[4] * 25.0);
}
// find B * m and D * m * mu^d3
s_x[i]->dw_t_visc = (s_x[i]->Jones_Dole[0] +
s_x[i]->Jones_Dole[1] * exp(-s_x[i]->Jones_Dole[2] * tc)) * t1;
Bc += s_x[i]->dw_t_visc;
// define f_I from the exponent of the D * m^d3 term...
if (s_x[i]->Jones_Dole[5] >= 1)
t2 = mu_x / 3 / s_x[i]->Jones_Dole[5];
else if (s_x[i]->Jones_Dole[5] > 0.4)
t2 = -0.8 / s_x[i]->Jones_Dole[5];
else
t2 = -1;
t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) *
t1 * (pow(mu_x, s_x[i]->Jones_Dole[5])*(1 + t2) + pow(t1 * f_z, s_x[i]->Jones_Dole[5])) / (2 + t2);
s_x[i]->dw_t_visc += t3;
Dc += t3;
//output_msg(sformatf("\t%s\t%e\t%e\t%e\n", s_x[i]->name, t1, Bc, Dc ));
}
// parms for A...
if ((l_z = s_x[i]->z) == 0)
continue;
Dw = s_x[i]->dw;
if (Dw)
{
Dw *= (0.89 / viscos_0 * tk_x / 298.15);
if (s_x[i]->dw_t)
Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15);
}
if (l_z < 0)
{
if (!strcmp(s_x[i]->name, "Cl-"))
// volumina for f_an...
{
V_Cl = s_x[i]->logk[vm_tc];
V_an += V_Cl * s_x[i]->moles;
m_an += s_x[i]->moles;
}
else// if (s_x[i]->Jones_Dole[6])
{
V_an += s_x[i]->logk[vm_tc] * s_x[i]->Jones_Dole[6] * s_x[i]->moles;
m_an += s_x[i]->moles;
}
if (Dw)
{
// anions for A...
m_min += s_x[i]->moles;
t1 = s_x[i]->moles * l_z;
eq_min -= t1;
eq_dw_min -= t1 / Dw;
}
}
else if (Dw)
{
// cations for A...
m_plus += s_x[i]->moles;
t1 = s_x[i]->moles * l_z;
eq_plus += t1;
eq_dw_plus += t1 / Dw;
}
}
if (m_plus && m_min && eq_dw_plus && eq_dw_min)
{
z1 = eq_plus / m_plus; z2 = eq_min / m_min;
D1 = eq_plus / eq_dw_plus; D2 = eq_min / eq_dw_min;
t1 = (D1 - D2) / (sqrt(D1 * z1 + D2 * z2) + sqrt((D1 + D2) * (z1 + z2)));
psi = (D1 * z2 + D2 * z1) / 4.0 - z1 * z2 * t1 * t1;
// Here A is A * viscos_0, avoids multiplication later on...
A = 4.3787e-14 * pow(tk_x, 1.5) / (sqrt(eps_r * (z1 + z2) / ((z1 > z2) ? z1 : z2)) * (D1 * D2)) * psi;
}
else
A = 0;
viscos = viscos_0 + A * sqrt((eq_plus + eq_min) / 2 / mass_water_aq_x);
if (m_an)
V_an /= m_an;
if (!V_Cl)
V_Cl = calc_vm_Cl();
if (V_an)
fan = 2 - V_an / V_Cl;
//else
// fan = 1;
viscos += viscos_0 * fan * (Bc + Dc);
if (viscos < 0)
{
viscos = viscos_0;
warning_msg("viscosity < 0, reset to viscosity of pure water");
}
for (i = 0; i < (int)this->s_x.size(); i++)
{
s_x[i]->dw_t_visc /= (Bc + Dc);
}
return viscos;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
calc_vm_Cl(void)
/* ---------------------------------------------------------------------- */
{
/*
* Calculate molar volume of Cl- with a Redlich type eqn:
Vm = Vm0(tc) + (Av / 2) * z^2 * I^0.5 + coef(tc) * I^(b4).
* Vm0(tc) is calc'd using supcrt parms, or from millero[0] + millero[1] * tc + millero[2] * tc^2
* for Av * z^2 * I^0.5, see Redlich and Meyer, Chem. Rev. 64, 221.
Av is in (cm3/mol)(mol/kg)^-0.5, = DH_Av.
If b_Av != 0, the extended DH formula is used: I^0.5 /(1 + b_Av * DH_B * I^0.5).
DH_Av and DH_B are from calc_dielectrics(tc, pa).
* coef(tc) = logk[vmi1] + logk[vmi2] / (TK - 228) + logk[vmi3] * (TK - 228).
* b4 = logk[vmi4], or
* coef(tc) = millero[3] + millero[4] * tc + millero[5] * tc^2
*/
LDBLE V_Cl = 0;
LDBLE pb_s = 2600. + patm_x * 1.01325, TK_s = tc_x + 45.15, sqrt_mu = sqrt(mu_x);
class species *s_ptr;
s_ptr = s_search("Cl-");
if (!s_ptr)
return V_Cl;
if (s_ptr->logk[vma1])
{
/* supcrt volume at I = 0... */
V_Cl = s_ptr->logk[vma1] + s_ptr->logk[vma2] / pb_s +
(s_ptr->logk[vma3] + s_ptr->logk[vma4] / pb_s) / TK_s -
s_ptr->logk[wref] * QBrn;
/* the ionic strength term * I^0.5... */
if (s_ptr->logk[b_Av] < 1e-5)
V_Cl += s_ptr->z * s_ptr->z * 0.5 * DH_Av * sqrt_mu;
else
{
/* limit the Debye-Hueckel slope by b... */
/* pitzer... */
//s_ptr->rxn_x.logk[vm_tc] += s_ptr->z * s_ptr->z * 0.5 * DH_Av *
// log(1 + s_ptr->logk[b_Av] * sqrt(mu_x)) / s_ptr->logk[b_Av];
/* extended DH... */
V_Cl += s_ptr->z * s_ptr->z * 0.5 * DH_Av *
sqrt_mu / (1 + s_ptr->logk[b_Av] * DH_B * sqrt_mu);
}
/* plus the volume terms * I... */
if (s_ptr->logk[vmi1] != 0.0 || s_ptr->logk[vmi2] != 0.0 || s_ptr->logk[vmi3] != 0.0)
{
LDBLE bi = s_ptr->logk[vmi1] + s_ptr->logk[vmi2] / TK_s + s_ptr->logk[vmi3] * TK_s;
if (s_ptr->logk[vmi4] == 1.0)
V_Cl += bi * mu_x;
else
V_Cl += bi * pow(mu_x, s_ptr->logk[vmi4]);
}
}
else if (s_ptr->millero[0])
{
/* Millero volume at I = 0... */
V_Cl = s_ptr->millero[0] + tc_x * (s_ptr->millero[1] + tc_x * s_ptr->millero[2]);
if (s_ptr->z)
{
/* the ionic strength terms... */
V_Cl += s_ptr->z * s_ptr->z * 0.5 * DH_Av * sqrt_mu +
(s_ptr->millero[3] + tc_x * (s_ptr->millero[4] + tc_x * s_ptr->millero[5])) * mu_x;
}
}
return V_Cl;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
flux_mcd(const char* species_name, int option)
/* ---------------------------------------------------------------------- */
{
class species* s_ptr;
double f = 0.0;
if (state == TRANSPORT && multi_Dflag)
{
s_ptr = s_search(species_name);
if (s_ptr != NULL && s_ptr->in != FALSE && s_ptr->type < EMINUS)
{
int n = cell_no;
std::map<int, std::map<std::string, class J_ij_save> >::iterator
j_map_iter = cell_J_ij.find(n);
if (j_map_iter != cell_J_ij.end())
{
std::map<std::string, class J_ij_save>::iterator
j_save_iter = j_map_iter->second.find(species_name);
if (j_save_iter != j_map_iter->second.end())
{
if (option == 1)
{
f = j_save_iter->second.flux_t;
}
else if (option == 2)
{
f = j_save_iter->second.flux_c;
}
}
}
}
}
return (f);
}
Reference in New Issue View Git Blame Copy Permalink