max/iphreeqc
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/iphreeqc.git synced 2025-12-16 08:38:23 +01:00
Code Issues Packages Projects Releases Wiki Activity

Tony's updates to allow fixed current.

Browse Source 
git-svn-id: svn://136.177.114.72/svn_GW/phreeqc3/branches/concrete@11133 1feff8c3-07ed-0310-ac33-dd36852eb9cd
This commit is contained in:
David L Parkhurst 2016-05-03 16:27:32 +00:00
parent 8604336e09
commit 03c9d56c09
4 changed files with 67 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Phreeqc.cpp
View File

@ -470,6 +470,7 @@ void Phreeqc::init(void)
multi_Dn = 0;
interlayer_tortf = 100.0;
cell_no = 0;
fix_current = 0.0;
/*----------------------------------------------------------------------
* Advection data
*---------------------------------------------------------------------- */
@ -1385,6 +1386,7 @@ Phreeqc::InternalCopy(const Phreeqc *pSrc)
multi_Dn = pSrc->multi_Dn;
interlayer_tortf = pSrc->interlayer_tortf;
cell_no = pSrc->cell_no;
fix_current = pSrc->fix_current;
/*----------------------------------------------------------------------
* Advection data
*---------------------------------------------------------------------- */

2
Phreeqc.h
View File

@ -1955,7 +1955,7 @@ protected:
int nmix, heat_nmix;
LDBLE heat_mix_f_imm, heat_mix_f_m;
int warn_MCD_X, warn_fixed_Surf;
LDBLE current_x, current_A; // current: coulomb * s, ampere
LDBLE current_x, current_A, fix_current; // current: coulomb / s, Ampere, fixed current (Ampere)
#ifdef PHREEQ98
int AutoLoadOutputFile, CreateToC;

25
readtr.cpp
View File

@ -91,9 +91,10 @@ read_transport(void)
"multi_d", /* 40 */
"interlayer_d", /* 41 */
"porosities", /* 42 */
"porosity" /* 43 */
"porosity", /* 43 */
"fix_current" /* 44 */
};
int count_opt_list = 44;
int count_opt_list = 45;
strcpy(file_name, "phreeqc.dmp");
/*
@ -665,6 +666,19 @@ read_transport(void)
}
opt_save = 42;
break;
case 44: /* fix_current */
if (copy_token(token, &next_char, &l) == DIGIT)
{
sscanf(token, SCANFORMAT, &fix_current);
fix_current = fabs(fix_current);
}
else
{
warning_msg("Expected the fixed value for the current (Ampere).");
fix_current = 0.0;
}
opt_save = OPTION_DEFAULT;
break;
}
if (return_value == EOF || return_value == KEYWORD)
break;
@ -938,13 +952,6 @@ read_transport(void)
cell_data[i].por_il = interlayer_Dpor;
}
#endif
//{
// for (int i = 0; i < all_cells; i++)
// {
// std::cerr << i << " " << cell_data[i].por << std::endl;
// }
//}
/*
* Calculate dump_modulus
*/

73
transport.cpp
View File

@ -105,8 +105,8 @@ transport(void)
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 = -1;
ct[i].J_ij_il_count_spec = -1;
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;
@ -186,6 +186,12 @@ transport(void)
}
}
if (fix_current && !dV_dcell)
{
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 (/*stag_data->count_stag || */!multi_Dflag || ishift)
@ -1815,7 +1821,7 @@ fill_spec(int l_cell_no)
}
if (correct_Dw)
{
if ((l_z = fabs(s_x[i]->z)) == 0)
if ((l_z = fabs(s_ptr->z)) == 0)
{ // first approximation for neutral species (HTO), but is viscosity dependent
l_z = 1;
l_g = -DH_A * (sqrt(mu_x) / (1 + sqrt(mu_x)));
@ -1912,7 +1918,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant)
current_x = sum_R = sum_Rd = 0.0;
if (dV_dcell && !stagnant)
find_current = loop_f_c = 1; // calculate J_ij once for all cells, find smallest j_x, next with this j_x.
find_current = loop_f_c = 1; // calculate J_ij once for all cells, loop to dV_dcell2.
else
find_current = loop_f_c = 0;
@ -1984,16 +1990,27 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant)
else
{
LDBLE dVc, j_0e;
// distribute dV_dcell according to relative resistance, calculate current_x
// 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 = (signbit(current_x) == 0 ? 1 : -1);
current_x = sign * fix_current;
j_0e = current_x - current_cells[0].dif;
}
dVc = j_0e * current_cells[0].R;
cell_data[1].potV = cell_data[0].potV + dVc;
current_x = j_0e + current_cells[0].dif;
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;
}
@ -2030,7 +2047,7 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant)
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 && i > 0 && !stagnant)
if (dV_dcell && (i > 0 || fix_current) && !stagnant)
{
use.Get_solution_ptr()->Set_potV(cell_data[icell].potV);
}
@ -2062,6 +2079,10 @@ multi_D(LDBLE DDt, int mobile_cell, int stagnant)
use.Get_solution_ptr()->Set_total_h(dummy + tot1_h);
dummy = use.Get_solution_ptr()->Get_total_o();
use.Get_solution_ptr()->Set_total_o(dummy + tot1_o);
if (i == count_cells && fix_current && !stagnant)
{
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);
@ -2329,29 +2350,31 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
ct[icell].v_m = (struct V_M *) free_check_null(ct[icell].v_m);
ct[icell].v_m_il = (struct V_M *) free_check_null(ct[icell].v_m_il);
//ct[icell].v_m = ct[icell].v_m_il = NULL;
if (stagnant)
if (!il_calcs)
{
if (!il_calcs && (cell_data[icell].por < multi_Dpor_lim
|| cell_data[jcell].por < multi_Dpor_lim))
return (OK);
}
else
{ /* regular column... */
if (icell == 0)
if (stagnant)
{
if (!il_calcs && cell_data[1].por < multi_Dpor_lim)
return (OK);
}
else if (icell == count_cells)
{
if (!il_calcs && cell_data[count_cells].por < multi_Dpor_lim)
if (cell_data[icell].por < multi_Dpor_lim
|| cell_data[jcell].por < multi_Dpor_lim)
return (OK);
}
else
{
if (!il_calcs && (cell_data[icell].por < multi_Dpor_lim
|| cell_data[jcell].por < multi_Dpor_lim))
{ /* 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;
}
return (OK);
}
}
}
@ -3283,7 +3306,7 @@ dV_dcell2:
ct[icell].J_ij_sum += ct[icell].v_m[i].z * ct[icell].J_ij[i].tot1;
}
}
// ensure that icell has dl water when checking negative conc's in MCD
// assure that icell has dl water when checking negative conc's in MCD
ct[icell].dl_s = dl_aq_i;
ct[jcell].dl_s = dl_aq_j;