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

Squashed 'phreeqcpp/' changes from 7284fed..89d028d

Browse Source 
89d028d Tony's fix for Valgrind error in Debye1.
9716b89 fix for new Valgrind problem. Added a newer CEMDATA database for the database collection.
eee3969 Merge branch 'master' into viscosity
b575463 Merge branch 'master' into viscosity
b74423e Issue resolution to be tested: Valgrind-Conditional jump or move depends on uninitialised value (SC) #48
67a69ae Tony changes 20240414, with correction to CH4 Vm. Changes to src. seaw_SC expanded.
34b880c Merge branch 'master' into viscosity
1819d3a Merge branch 'master' into viscosity

git-subtree-dir: phreeqcpp
git-subtree-split: 89d028d8321339d90539869f46da1143f23d4025
This commit is contained in:
Darth Vader 2024-04-18 18:16:55 +00:00
parent 1ebe8191c2
commit c378c74d7f
6 changed files with 158 additions and 152 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
PBasic.cpp
View File

@ -3380,20 +3380,22 @@ factor(struct LOC_exec * LINK)
case toksetdiff_c:
{
double d;
double d, d_v_d = 0;
require(toklp, LINK);
const char* str = stringfactor(STR1, LINK);
require(tokcomma, LINK);
// double arugument
d = realexpr(LINK);
if (LINK->t != NULL && LINK->t->kind == tokcomma)
{
LINK->t = LINK->t->next;
d_v_d = realexpr(LINK);
}
require(tokrp, LINK);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->setdiff_c(str, d);
//PhreeqcPtr->PHRQ_free((void *) str);
n.UU.val = (parse_all) ? 1 : PhreeqcPtr->setdiff_c(str, d, d_v_d);
}
break;

2
Phreeqc.h
View File

@ -110,7 +110,7 @@ public:
LDBLE aqueous_vm(const char* species_name);
LDBLE phase_vm(const char* phase_name);
LDBLE diff_c(const char* species_name);
LDBLE setdiff_c(const char* species_name, double d);
LDBLE setdiff_c(const char * species_name, double d, double d_v_d);
LDBLE flux_mcd(const char* species_name, int option);
LDBLE sa_declercq(double type, double sa, double d, double m, double m0, double gfw);
LDBLE calc_SC(void);

137
basicsubs.cpp
View File

@ -178,7 +178,6 @@ sa_declercq(double sa_type, double Sa, double d, double m, double m0, double gfw
error_msg(error_string, CONTINUE);
input_error++;
return (MISSING);
}
/* ---------------------------------------------------------------------- */
@ -187,97 +186,58 @@ diff_c(const char *species_name)
/* ---------------------------------------------------------------------- */
{
class species *s_ptr;
LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av;
sqrt_mu = sqrt(mu_x);
LDBLE Dw;
s_ptr = s_search(species_name);
if (s_ptr == NULL)
return(0);
if ((Dw = s_ptr->dw) == 0)
{
if (correct_Dw)
Dw = default_Dw;
else
return(0);
}
return(0);
if (correct_Dw)
{
if ((l_z = fabs(s_ptr->z)) == 0)
{
//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
}
else
{
if (print_viscosity)
{
a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6);
a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73);
av = (s_ptr->dw_a_visc ? pow((viscos_0 / viscos), s_ptr->dw_a_visc) : 1);
a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : s_ptr->dw_a_visc ? 1 : pow(mu_x, 0.75));
}
else
{
a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6);
a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73);
av = 1.0;
a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : pow(mu_x, 0.75));
}
ka = DH_B * a2 * sqrt_mu / (1 + a3);
ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka));
Dw *= ff;
}
calc_SC();
Dw = s_ptr->dw_corr;
}
if (tk_x != 298.15 && s_ptr->dw_t)
Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
else
{
if (tk_x != 298.15 && s_ptr->dw_t)
Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
s_ptr->dw_corr = Dw * viscos_0_25 / viscos_0;
return s_ptr->dw_corr;
Dw *= viscos_0_25 / viscos_0;
}
if (s_ptr->dw_a_v_dif)
Dw *= pow(viscos_0 / viscos, s_ptr->dw_a_v_dif);
return Dw;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
setdiff_c(const char *species_name, double d)
setdiff_c(const char *species_name, double d, double d_v_d)
/* ---------------------------------------------------------------------- */
{
class species *s_ptr;
LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av;
sqrt_mu = sqrt(mu_x);
LDBLE Dw;
s_ptr = s_search(species_name);
if (s_ptr == NULL)
return(0);
Dw = s_ptr->dw = d;
s_ptr->dw_a_v_dif = d_v_d;
if (correct_Dw)
{
if ((l_z = fabs(s_ptr->z)) == 0)
{
//l_z = 1; // only a 1st approximation for correct_Dw in electrical field
}
else
{
if (print_viscosity)
{
a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6);
a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73);
av = (s_ptr->dw_a_visc ? pow((viscos_0 / viscos), s_ptr->dw_a_visc) : 1);
a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : 1);
}
else
{
a = (s_ptr->dw_a ? s_ptr->dw_a : 1.6);
a2 = (s_ptr->dw_a2 ? s_ptr->dw_a2 : 4.73);
av = 1.0;
a3 = (s_ptr->dw_a3 ? pow(mu_x, s_ptr->dw_a3) : s_ptr->dw_a_visc ? 1 : pow(mu_x, 0.75));
}
ka = DH_B * a2 * sqrt_mu / (1 + a3);
ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka));
Dw *= ff;
}
calc_SC();
Dw = s_ptr->dw_corr;
}
if (tk_x != 298.15 && s_ptr->dw_t)
else
{
if (tk_x != 298.15 && s_ptr->dw_t)
Dw *= exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
s_ptr->dw_corr = Dw * viscos_0_25 / viscos_0;
return s_ptr->dw_corr;
Dw *= viscos_0_25 / viscos_0;
}
if (d_v_d)
Dw *= pow(viscos_0 / viscos, d_v_d);
return Dw;
}
/* ---------------------------------------------------------------------- */
LDBLE Phreeqc::
@ -287,6 +247,7 @@ calc_SC(void)
class species *s_ptr;
int i;
LDBLE ka, l_z, Dw, ff, sqrt_mu, a, a2, a3, av, v_Cl = 1;
SC = 0;
sqrt_mu = sqrt(mu_x);
bool Falk = false;
s_ptr = s_search("H+");
@ -294,7 +255,6 @@ calc_SC(void)
return(0);
else if (s_ptr->dw_a3 > 24) Falk = true;
SC = 0;
//LDBLE ta1, ta2, ta3, ta4;
//for (i = 0; i < (int)this->s_x.size(); i++)
//{
@ -320,7 +280,10 @@ calc_SC(void)
if (correct_Dw)
Dw = default_Dw;
else
{
s_x[i]->dw_corr = 0;
continue;
}
}
if (s_x[i]->lm < min_dif_LM)
continue;
@ -343,27 +306,6 @@ calc_SC(void)
}
else
{
//if (s_x[i]->dw_a2)
// ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.259));
// //ka = DH_B * s_x[i]->dw_a2 * sqrt_mu / (1 + pow(mu_x, 0.75));
//else
//{
// ka = DH_B * 4.73 * sqrt_mu / (1 + pow(mu_x, 0.259));
// //ka = DH_B * ta1 * sqrt_mu / (1 + pow(mu_x, ta2));
// //ka = DH_B * ta1 * sqrt_mu / (1 + mu_x / ta2);
//}
//if (s_x[i]->dw_a)
//{
// ff = exp(-s_x[i]->dw_a * DH_A * l_z * sqrt_mu / (1 + ka));
// if (print_viscosity && s_x[i]->dw_a_visc)
// ff *= pow((viscos_0 / viscos), s_x[i]->dw_a_visc);
//}
//else
//{
// ff = exp(-1.6 * DH_A * l_z * sqrt_mu / (1 + ka));
// //ff = exp(-ta3 * DH_A * l_z * sqrt_mu / (1 + ka));
//}
//Dw *= ff;
s_ptr = s_x[i];
if (print_viscosity)
{
@ -384,6 +326,8 @@ calc_SC(void)
}
Dw *= ff;
if (correct_Dw)
s_x[i]->dw_corr = Dw;
s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw;
SC += s_x[i]->dw_t_SC;
}
@ -416,8 +360,13 @@ calc_SC(void)
continue;
if ((Dw = s_x[i]->dw) == 0)
{
if (correct_Dw) Dw = default_Dw; // or charge based...Dw = l_z > 0 ? 1.6e-9 / l_z : 2e-9 / -l_z;
else continue;
if (correct_Dw)
Dw = default_Dw; // or charge based...Dw = l_z > 0 ? 1.6e-9 / l_z : 2e-9 / -l_z;
else
{
s_x[i]->dw_corr = 0;
continue;
}
}
if (tk_x != 298.15 && s_x[i]->dw_t)
Dw *= exp(s_x[i]->dw_t / tk_x - s_x[i]->dw_t / 298.15);
@ -496,6 +445,8 @@ calc_SC(void)
ff = av * exp(-a * DH_A * l_z * sqrt_mu / (1 + ka));
Dw *= ff;
if (correct_Dw)
s_x[i]->dw_corr = Dw;
s_x[i]->dw_t_SC = s_x[i]->moles / mass_water_aq_x * l_z * l_z * Dw;
SC += s_x[i]->dw_t_SC * 1e3 * Dw_SC;
}
@ -540,7 +491,9 @@ calc_SC(void)
//t1 = (Dw - B2 * l_z * sqrt_mu / (1 + ka)) *
// (1 - B1 * sqrt_mu / ((1 + ka) *(1 + ka * sqrt_q + ka * ka / 6))); // S.cm2/eq / (kgw/L)
if (av)
t1 *= pow(viscos_0 / viscos, av);
t1 *= pow(viscos_0 / viscos, av);
if (correct_Dw)
s_x[i]->dw_corr *= t1 / Dw;
// fractional contribution in mu, and correct for charge imbalance
a2 = 2 / (eq_plus + eq_min);

24
global_structures.h
View File

@ -14,7 +14,7 @@
# define NAN nan("1")
# endif
#endif
#define MISSING -9999.999
#define MISSING -9999.999
#include "NA.h" /* NA = not available */
#define F_C_MOL 96493.5 /* C/mol or joule/volt-eq */
@ -712,18 +712,18 @@ public:
secondary = NULL;
gfw = 0; // gram formula wt of species
z = 0; // charge of species
// tracer diffusion coefficient in water at 25oC, m2/s
dw = 0;
// correct Dw for temperature: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15)
dw_t = 0;
dw = 0; // tracer diffusion coefficient in water at 25oC, m2/s
dw_t = 0; // correct Dw for temperature: Dw(TK) = Dw(298.15) * exp(dw_t / TK - dw_t / 298.15)
// parms for calc'ng SC = SC0 * exp(-dw_a * z * mu^0.5 / (1 + DH_B * dw_a2 * mu^0.5) / (1 + mu^dw_a3))
// with DHO: ka = DH_B * dw_a * (1 + DD(V_apparent)^dw_a2 * sqrt_mu, dw_a3 is a switch, see calc_SC in PBasic
dw_a = 0;
dw_a2 = 0;
dw_a3 = 0;
dw_a_visc = 0; // viscosity correction of SC
dw_a_visc = 0; // exponent in viscosity correction of SC
dw_a_v_dif = 0; // exponent in viscosity correction of D, the diffusion coefficient of the species
dw_t_SC = 0; // contribution to SC, for calc'ng transport number with BASIC
dw_t_visc = 0; // contribution to viscosity
dw_corr = 0; // dw corrected for TK and mu
dw_corr = 0; // dw corrected for mu and TK
erm_ddl = 0; // enrichment factor in DDL
equiv = 0; // equivalents in exchange species
alk = 0; // alkalinity of species, used for cec in exchange
@ -784,6 +784,7 @@ public:
LDBLE dw_a2;
LDBLE dw_a3;
LDBLE dw_a_visc;
LDBLE dw_a_v_dif;
LDBLE dw_t_SC;
LDBLE dw_t_visc;
LDBLE dw_corr;
@ -1489,6 +1490,8 @@ public:
Dwt = 0;
// temperature factor for Dw
dw_t = 0;
// viscosity factor for Dw
dw_a_v_dif = 0;
// enrichment factor in ddl
erm_ddl = 0;
}
@ -1502,6 +1505,7 @@ public:
LDBLE z;
LDBLE Dwt;
LDBLE dw_t;
LDBLE dw_a_v_dif;
LDBLE erm_ddl;
};
@ -1517,7 +1521,9 @@ public:
count_exch_spec = 0;
// total moles of X-, max X- in transport step in sol_D[1], tk
exch_total = 0, x_max = 0, tk_x = 0;
// (tk_x * viscos_0_25) / (298 * viscos)
// (tk_x * viscos_0_25) / (298 * viscos_0)
viscos_f0 = 0;
// (viscos_0) / (298 * viscos)
viscos_f = 0;
spec = NULL;
spec_size = 0;
@ -1525,7 +1531,7 @@ public:
int count_spec;
int count_exch_spec;
LDBLE exch_total, x_max, tk_x;
LDBLE viscos_f;
LDBLE viscos_f0, viscos_f;
class spec* spec;
int spec_size;
};

6
read.cpp
View File

@ -5499,9 +5499,9 @@ read_species(void)
input_error++;
break;
}
s_ptr->dw_t = 0; s_ptr->dw_a = 0; s_ptr->dw_a2 = 0; s_ptr->dw_a3 = 0; s_ptr->dw_a_visc = 0;
i = sscanf(next_char, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT,
&s_ptr->dw, &s_ptr->dw_t, &s_ptr->dw_a, &s_ptr->dw_a2, &s_ptr->dw_a_visc, &s_ptr->dw_a3);
s_ptr->dw_t = 0; s_ptr->dw_a = 0; s_ptr->dw_a2 = 0; s_ptr->dw_a3 = 0; s_ptr->dw_a_visc = 0; s_ptr->dw_a_v_dif = 0;
i = sscanf(next_char, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT,
&s_ptr->dw, &s_ptr->dw_t, &s_ptr->dw_a, &s_ptr->dw_a2, &s_ptr->dw_a_visc, &s_ptr->dw_a3, &s_ptr->dw_a_v_dif);
if (i < 1)
{
input_error++;

129
transport.cpp
View File

@ -115,6 +115,7 @@ transport(void)
sol_D[i].count_exch_spec = 0;
sol_D[i].exch_total = 0;
sol_D[i].x_max = 0;
sol_D[i].viscos_f0 = 1.0;
sol_D[i].viscos_f = 1.0;
sol_D[i].tk_x = 298.15;
sol_D[i].spec = NULL;
@ -1640,14 +1641,14 @@ init_heat_mix(int l_nmix)
{
if (implicit)
{
LDBLE viscos_f;
LDBLE viscos_f0;
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);
viscos_f0 = sol_D[i - 1].viscos_f0 * exp(heat_diffc / sol_D[i - 1].tk_x - heat_diffc / 298.15);
viscos_f0 += sol_D[i].viscos_f0 * exp(heat_diffc / sol_D[i].tk_x - heat_diffc / 298.15);
heat_mix_array[i - 1] *= (viscos_f0 / 2);
}
}
else
@ -1821,7 +1822,7 @@ fill_spec(int l_cell_no, int ref_cell)
class master *master_ptr;
LDBLE dum, dum2;
LDBLE lm;
LDBLE por, por_il, viscos_f, viscos_il_f, viscos;
LDBLE por, por_il, viscos_f0, viscos_f, viscos_il_f0, viscos_il_f, viscos;
bool x_max_done = false;
std::set <std::string> loc_spec_names;
@ -1860,13 +1861,14 @@ fill_spec(int l_cell_no, int ref_cell)
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].dw_a_v_dif = 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;
viscos_f0 = viscos_il_f0 = viscos_f = viscos_il_f = 1.0;
if (l_cell_no == 0)
{
por = cell_data[1].por;
@ -1883,10 +1885,10 @@ fill_spec(int l_cell_no, int ref_cell)
por_il = cell_data[l_cell_no].por_il;
}
if (por < multi_Dpor_lim)
por = viscos_f = 0.0;
por = viscos_f0 = viscos_f = 0.0;
if (por_il < interlayer_Dpor_lim)
por_il = viscos_il_f = 0.0;
por_il = viscos_il_f0 = 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
@ -1898,10 +1900,14 @@ fill_spec(int l_cell_no, int ref_cell)
/*
* 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;
dum = viscos_0_25 / viscos_0;
dum2 = viscos_0 / viscos;
viscos_f0 *= dum;
viscos_il_f0 *= dum;
sol_D[l_cell_no].viscos_f0 = dum;
viscos_f *= dum2;
viscos_il_f *= dum2;
sol_D[l_cell_no].viscos_f = dum2;
count_spec = count_exch_spec = 0;
/*
@ -2002,6 +2008,10 @@ fill_spec(int l_cell_no, int ref_cell)
else
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr2->dw * viscos_il_f;
}
if (s_ptr2->dw_a_v_dif)
sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = s_ptr2->dw_a_v_dif;
else
sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = 0.0;
//if (implicit) // && name_ret.second && (l_cell_no > 1 || (l_cell_no == 1 && bcon_first != 2)))
//{
@ -2089,22 +2099,29 @@ fill_spec(int l_cell_no, int ref_cell)
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;
sol_D[l_cell_no].spec[count_spec].Dwt = default_Dw;
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;
exp(s_ptr->dw_t / tk_x - s_ptr->dw_t / 298.15);
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;
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw;
}
if (s_ptr->dw_a_v_dif)
{
sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = s_ptr->dw_a_v_dif;
sol_D[l_cell_no].spec[count_spec].Dwt *= pow(viscos_0 / viscos, s_ptr->dw_a_v_dif);
}
else
sol_D[l_cell_no].spec[count_spec].dw_a_v_dif = 0.0;
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;
sol_D[l_cell_no].spec[count_spec].Dwt = s_ptr->dw_corr;
}
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);
@ -4272,7 +4289,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
{
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1;
dum = ct[icell].visc1;
if (sol_D[icell].spec[i].dw_a_v_dif)
dum = pow(dum, sol_D[icell].spec[i].dw_a_v_dif);
g_i *= sol_D[icell].spec[i].erm_ddl / dum;
//g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1;
}
if (dl_aq2)
{
@ -4295,7 +4316,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
}
}
}
g_j *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc2;
dum = ct[icell].visc2;
if (sol_D[icell].spec[i].dw_a_v_dif)
dum = pow(dum, sol_D[icell].spec[i].dw_a_v_dif);
g_j *= sol_D[icell].spec[i].erm_ddl / dum;
//g_j *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc2;
}
}
@ -4305,11 +4330,13 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
b_j *= sol_D[icell].spec[i].Dwt;
else
{
dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f;
dum2 = sol_D[icell].spec[i].Dwt / sol_D[icell].viscos_f0;
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;
dum2 *= sol_D[jcell].viscos_f0;
b_j *= dum2;
}
if (sol_D[icell].spec[i].dw_a_v_dif)
b_j *= pow(sol_D[jcell].viscos_f / sol_D[icell].viscos_f, sol_D[icell].spec[i].dw_a_v_dif);
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
k++;
@ -4388,7 +4415,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
}
}
}
g_i *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc1;
dum = ct[icell].visc1;
if (sol_D[jcell].spec[j].dw_a_v_dif)
dum = pow(dum, sol_D[icell].spec[j].dw_a_v_dif);
g_i *= sol_D[jcell].spec[j].erm_ddl / dum;
//g_i *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc1;
}
if (dl_aq2)
{
@ -4396,7 +4427,12 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
{
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2;
dum = ct[icell].visc2;
if (sol_D[jcell].spec[j].dw_a_v_dif)
dum = pow(dum, sol_D[jcell].spec[j].dw_a_v_dif);
g_j *= sol_D[jcell].spec[j].erm_ddl / dum;
//g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2;
}
}
b_i = A1;
@ -4405,11 +4441,14 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
b_i *= sol_D[jcell].spec[j].Dwt;
else
{
dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f;
dum2 = sol_D[jcell].spec[j].Dwt / sol_D[jcell].viscos_f0;
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;
dum2 *= sol_D[icell].viscos_f0;
b_i *= dum2;
}
if (sol_D[icell].spec[i].dw_a_v_dif)
b_i *= pow(sol_D[icell].viscos_f / sol_D[jcell].viscos_f, sol_D[icell].spec[i].dw_a_v_dif);
calc_b_ij(icell, jcell, k, b_i, b_j, g_i, g_j, f_free_i, f_free_j, stagnant);
k++;
@ -4482,7 +4521,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
{
g_i += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1;
dum = ct[icell].visc1;
if (sol_D[icell].spec[i].dw_a_v_dif)
dum = pow(dum, sol_D[icell].spec[i].dw_a_v_dif);
g_i *= sol_D[icell].spec[i].erm_ddl / dum;
//g_i *= sol_D[icell].spec[i].erm_ddl / ct[icell].visc1;
}
if (dl_aq2)
{
@ -4490,7 +4533,11 @@ find_J(int icell, int jcell, LDBLE mixf, LDBLE DDt, int stagnant)
{
g_j += it_sc->Get_z_gMCD_map()[ct[icell].v_m[k].z];
}
g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2;
dum = ct[icell].visc2;
if (sol_D[jcell].spec[j].dw_a_v_dif)
dum = pow(dum, sol_D[jcell].spec[j].dw_a_v_dif);
g_j *= sol_D[jcell].spec[j].erm_ddl / dum;
//g_j *= sol_D[jcell].spec[j].erm_ddl / ct[icell].visc2;
}
}
b_i = A1 * sol_D[icell].spec[i].Dwt;
@ -4874,12 +4921,11 @@ Step (from cell 1 to count_cells + 1):
LDBLE lav, A_ij, por, Dp1, Dp2;
cxxMix * mix_ptr;
cxxSurface *surface_ptr1, *surface_ptr2;
LDBLE viscos_f;
LDBLE viscos_f0;
/*
* 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;
viscos_f0 = viscos_0_25 / viscos_0;
//n1 = 0;
//n2 = n1 + 1;
@ -5089,7 +5135,7 @@ Step (from cell 1 to count_cells + 1):
Dp2 = (Dp2 + Dp1) / 2;
/* and the mixing factor... */
mcd_mixf = Dp2 * viscos_f * A_ij * DDt / lav;
mcd_mixf = Dp2 * viscos_f0 * A_ij * DDt / lav;
}
mixf = mixf_store + mcd_mixf;
@ -5165,7 +5211,7 @@ Step (from cell 1 to count_cells + 1):
Dp1 = (Dp1 + Dp2) / 2;
/* and the mixing factor... */
mcd_mixf = Dp1 * viscos_f * A_ij * DDt / lav;
mcd_mixf = Dp1 * viscos_f0 * A_ij * DDt / lav;
}
mixf = mixf_store + mcd_mixf;
@ -5545,12 +5591,11 @@ diff_stag_surf(int mobile_cell)
LDBLE Dp1, Dp2;
cxxMix *mix_ptr;
cxxSurface *surface_ptr1, *surface_ptr2;
LDBLE viscos_f;
LDBLE viscos_f0;
/*
* 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;
viscos_f0 = viscos_0_25 / viscos_0;
cxxSurface surface_n1, surface_n2;
cxxSurface *surface_n1_ptr = &surface_n1;
@ -5698,7 +5743,7 @@ diff_stag_surf(int mobile_cell)
if (multi_Dflag)
{
Dp2 = comp_k_ptr->Get_Dw() *
pow(cell_data[i2].por, multi_Dn) * viscos_f;
pow(cell_data[i2].por, multi_Dn) * viscos_f0;
Dp1 = 0;
if (surf1)
{
@ -5708,7 +5753,7 @@ diff_stag_surf(int mobile_cell)
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;
Dp1 = comp_k1_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f0;
break;
}
}
@ -5757,7 +5802,7 @@ diff_stag_surf(int mobile_cell)
/* find diffusion coefficients of surfaces... */
if (multi_Dflag)
{
Dp1 = comp_k_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f;
Dp1 = comp_k_ptr->Get_Dw() * pow(cell_data[i1].por, multi_Dn) * viscos_f0;
Dp2 = 0;
if (surf2)
@ -5768,7 +5813,7 @@ diff_stag_surf(int mobile_cell)
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;
Dp2 = comp_k1_ptr->Get_Dw() * pow(cell_data[i2].por, multi_Dn) * viscos_f0;
break;
}
}
@ -6121,7 +6166,7 @@ viscosity(cxxSurface *surf_ptr)
t2 = -1;
t3 = (s_x[i]->Jones_Dole[3] * exp(-s_x[i]->Jones_Dole[4] * tc)) *
t1 * (pow(l_mu_x, s_x[i]->Jones_Dole[5])*(1 + t2) + pow(t1 * f_z, s_x[i]->Jones_Dole[5])) / (2 + t2);
if (t3 < -1e-5) // add this check
if (t3 < -1e-5)
t3 = 0;
Dc += t3;
if (!surf_ptr) s_x[i]->dw_t_visc = dw_t_visc + t3;
@ -6186,10 +6231,10 @@ viscosity(cxxSurface *surf_ptr)
V_an /= m_an;
if (!V_Cl)
V_Cl = calc_vm_Cl();
if (V_an)
if (V_an && V_Cl)
fan = 2 - V_an / V_Cl;
//else
// fan = 1;
else
fan = 1;
if (Dc < 0)
Dc = 0; // provisional...
viscos += viscos_0 * fan * (Bc + Dc);