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

best I could do for H2S while maintaining old tests. Used INCREMENTAL reactions

Browse Source
This commit is contained in:
David Parkhurst 2021-07-13 14:24:00 -06:00
parent 8be1ba8327
commit 13ec2fcd3e
4 changed files with 92 additions and 196 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
gases.cpp
View File

@ -101,7 +101,7 @@ build_fixed_volume_gas(void)
} }
/* All elements in gas */ /* All elements in gas */
for (j = 0; j < count_elts; j++) for (j = 0; j < (int) count_elts; j++)
{ {
unknown_ptr = NULL; unknown_ptr = NULL;
if (strcmp(elt_list[j].elt->name, "H") == 0) if (strcmp(elt_list[j].elt->name, "H") == 0)
@ -149,7 +149,7 @@ build_fixed_volume_gas(void)
output_msg(sformatf( "\n\tJacobian summations %s.\n\n", output_msg(sformatf( "\n\tJacobian summations %s.\n\n",
phase_ptr->name)); phase_ptr->name));
} }
for (j = 0; j < count_elts; j++) for (j = 0; j < (int) count_elts; j++)
{ {
unknown_ptr = NULL; unknown_ptr = NULL;
if (strcmp(elt_list[j].elt->name, "H") == 0) if (strcmp(elt_list[j].elt->name, "H") == 0)
@ -438,10 +438,14 @@ calc_PR(void)
a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217
else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)"))
a_aa *= 0.81; a_aa *= 0.81;
else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)")) else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)") || !strcmp(phase_ptr1->name, "Methane(g)"))
a_aa *= 0.51; a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)"))
a_aa *= 0.51; a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "Ethane(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr1->name, "Propane(g)"))
a_aa *= 0.45;
} }
if (!strcmp(phase_ptr1->name, "H2O(g)")) if (!strcmp(phase_ptr1->name, "H2O(g)"))
{ {
@ -449,10 +453,14 @@ calc_PR(void)
a_aa *= 0.81; a_aa *= 0.81;
else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)"))
a_aa *= 0.81; a_aa *= 0.81;
else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)")) else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)") || !strcmp(phase_ptr->name, "Methane(g)"))
a_aa *= 0.51; a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)"))
a_aa *= 0.51; a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "Ethane(g)"))
a_aa *= 0.51;
else if (!strcmp(phase_ptr->name, "Propane(g)"))
a_aa *= 0.45;
} }
a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa; a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa;
a_aa_sum2 += phase_ptr1->fraction_x * a_aa; a_aa_sum2 += phase_ptr1->fraction_x * a_aa;

12
model.cpp
View File

@ -2613,13 +2613,12 @@ calc_gas_pressures(void)
V_m = (2. * gas_phase_ptr->Get_v_m() + V_m) / 3; V_m = (2. * gas_phase_ptr->Get_v_m() + V_m) / 3;
else else
V_m = (1. * gas_phase_ptr->Get_v_m() + V_m) / 2; V_m = (1. * gas_phase_ptr->Get_v_m() + V_m) / 2;
if (iterations > 99 && numerical_fixed_volume == false) if ((pitzer_model || iterations > 99) && numerical_fixed_volume == false)
{ {
//V_m *= 1; /* debug */ //V_m *= 1; /* debug */
numerical_fixed_volume = true; numerical_fixed_volume = true;
//switch_numerical = true; //switch_numerical = true;
warning_msg if (!pitzer_model) warning_msg("Numerical method failed, switching to numerical derivatives.");
("Numerical method failed, switching to numerical derivatives.");
prep(); prep();
//switch_numerical = false; //switch_numerical = false;
} }
@ -3227,12 +3226,12 @@ reset(void)
} }
else if (x[i]->type == GAS_MOLES) else if (x[i]->type == GAS_MOLES)
{ {
up = 10. * x[i]->moles; up = 1000. * x[i]->moles;
if (up <= 0.0) if (up <= 0.0)
up = 1e-1; up = 1e-1;
if (up >= 1.0) if (up >= 1.0)
up = 1.; up = 1.;
down = 0.3*x[i]->moles; down = x[i]->moles;
} }
else if (x[i]->type == SS_MOLES) else if (x[i]->type == SS_MOLES)
{ {
@ -5425,7 +5424,8 @@ numerical_jacobian(void)
case GAS_MOLES: case GAS_MOLES:
if (gas_in == FALSE) if (gas_in == FALSE)
continue; continue;
d2 = d * 20 * x[i]->moles; d2 = (x[i]->moles > 1 ? 1 : 20);
d2 *= d * x[i]->moles;
if (d2 < 1e-14) if (d2 < 1e-14)
d2 = 1e-14; d2 = 1e-14;
x[i]->moles += d2; x[i]->moles += d2;

257
pitzer.cpp
View File

@ -1955,140 +1955,35 @@ jacobian_pz(void)
LDBLE d, d1, d2; LDBLE d, d1, d2;
int i, j; int i, j;
Restart: calculating_deriv = 1;
//molalities(TRUE);
//if (full_pitzer == TRUE)
// pitzer();
//mb_sums();
//residuals();
if (use.Get_gas_phase_ptr() != NULL) if (use.Get_gas_phase_ptr() != NULL)
{ {
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
for (size_t i = 0; i < use.Get_gas_phase_ptr()->Get_gas_comps().size(); i++) base_gas_phase = *gas_phase_ptr;
base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
{ {
const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]); const cxxGasComp* gas_comp_ptr = &(gas_phase_ptr->Get_gas_comps()[i]);
class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE); class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE);
phase_ptrs.push_back(phase_ptr); phase_ptrs.push_back(phase_ptr);
base_phases[i] = *phase_ptr;
} }
//base_gas_moles.resize(gas_phase_ptr->Get_gas_comps().size(), 0.0);
base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
} }
// Debugging Restart:
//output_msg("\n\nBefore jacobian_pz===================================================\n");
//print_all();
calculating_deriv = 1;
size_t pz_max_unknowns = max_unknowns; size_t pz_max_unknowns = max_unknowns;
//k_temp(tc_x, patm_x); //k_temp(tc_x, patm_x);
//if (full_pitzer == TRUE) if (full_pitzer == TRUE)
//{ {
// molalities(TRUE); molalities(TRUE);
// pitzer(); pitzer();
// residuals(); residuals();
//} }
base.resize(count_unknowns); base.resize(count_unknowns);
base_x.resize(count_unknowns, 0.0); for (i = 0; i < count_unknowns; i++)
for (size_t i1 = 0; i1 < count_unknowns; i1++)
{ {
//base[i1] = residual[i1]; base[i] = residual[i];
switch (x[i1]->type)
{
case MB:
case ALK:
case CB:
case SOLUTION_PHASE_BOUNDARY:
case EXCH:
case SURFACE:
case SURFACE_CB:
case SURFACE_CB1:
case SURFACE_CB2:
case AH2O:
base_x[i1] = x[i1]->master[0]->s->la;
break;
case PITZER_GAMMA:
base_x[i1] = x[i1]->s->lg;
break;
case MH2O:
base_x[i1] = mass_water_aq_x;
base_mass_water_bulk_x = mass_water_bulk_x;
base_moles_h2o = x[i1]->master[0]->s->moles;
break;
case MH:
base_x[i1] = s_eminus->la;
break;
case GAS_MOLES:
base_x[i1] = x[i1]->moles;
for (size_t g = 0; g < use.Get_gas_phase_ptr()->Get_gas_comps().size(); g++)
{
//base_gas_moles[g] = phase_ptrs[g]->moles_x;
base_phases[g] = *phase_ptrs[g];
}
break;
case MU:
base_x[i1] = mu_x;
break;
case PP:
continue;
break;
case SS_MOLES:
break;
}
} }
//molalities(TRUE); d = 0.0001;
//if (full_pitzer == TRUE)
// pitzer();
//mb_sums();
//residuals();
for (size_t i1 = 0; i1 < count_unknowns; i1++)
{
base[i1] = residual[i1];
}
for (size_t i1 = 0; i1 < count_unknowns; i1++)
{
switch (x[i1]->type)
{
case MB:
case ALK:
case CB:
case SOLUTION_PHASE_BOUNDARY:
case EXCH:
case SURFACE:
case SURFACE_CB:
case SURFACE_CB1:
case SURFACE_CB2:
case AH2O:
x[i1]->master[0]->s->la = base_x[i1];
break;
case PITZER_GAMMA:
x[i1]->s->lg = base_x[i1];
break;
case MH2O:
mass_water_aq_x = base_x[i1];
mass_water_bulk_x = base_mass_water_bulk_x;
x[i1]->master[0]->s->moles = base_moles_h2o;
break;
case MH:
s_eminus->la = base_x[i1];
break;
case GAS_MOLES:
x[i1]->moles = base_x[i1];
for (size_t g = 0; g < base_gas_moles.size(); g++)
{
//phase_ptrs[g]->moles_x = base_gas_moles[g];
*phase_ptrs[g] = base_phases[g];
}
break;
case MU:
mu_x = base_x[i1];
gammas_pz(false);
break;
case PP:
continue;
break;
case SS_MOLES:
break;
}
}
d = 0.001;
d1 = d * LOG_10; d1 = d * LOG_10;
d2 = 0; d2 = 0;
for (i = 0; i < count_unknowns; i++) for (i = 0; i < count_unknowns; i++)
@ -2148,10 +2043,9 @@ Restart:
case GAS_MOLES: case GAS_MOLES:
if (gas_in == FALSE) if (gas_in == FALSE)
continue; continue;
d2 = d * 30 * x[i]->moles; d2 = d * x[i]->moles;
d2 = (d2 < ineq_tol ? ineq_tol : d2); if (d2 < 1e-14)
//if (d2 < 1e-14) d2 = 1e-14;
// d2 = 1e-14;
x[i]->moles += d2; x[i]->moles += d2;
break; break;
case MU: case MU:
@ -2197,74 +2091,67 @@ Restart:
if (x[i]->type == MH2O) // DL_pitz if (x[i]->type == MH2O) // DL_pitz
my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] *= mass_water_aq_x; my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] *= mass_water_aq_x;
} }
for (size_t i1 = 0; i1 < count_unknowns; i1++) switch (x[i]->type)
{ {
switch (x[i1]->type) case MB:
case ALK:
case CB:
case SOLUTION_PHASE_BOUNDARY:
case EXCH:
case SURFACE:
case SURFACE_CB:
case SURFACE_CB1:
case SURFACE_CB2:
case AH2O:
x[i]->master[0]->s->la -= d;
break;
case MH:
s_eminus->la -= d;
if (my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] == 0)
{ {
case MB: my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] =
case ALK: exp(s_h2->lm * LOG_10) * 2;
case CB: }
case SOLUTION_PHASE_BOUNDARY: break;
case EXCH: case PITZER_GAMMA:
case SURFACE: x[i]->s->lg -= d;
case SURFACE_CB: break;
case SURFACE_CB1: case MH2O:
case SURFACE_CB2: //mass_water_aq_x /= (1 + d);
case AH2O: //x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
x[i1]->master[0]->s->la = base_x[i1]; //break;
break; //DL_pitz
case PITZER_GAMMA: mass_water_aq_x -= d1;
x[i1]->s->lg = base_x[i1]; if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL)
break; mass_water_bulk_x -= d1;
case MH2O: x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water;
mass_water_aq_x = base_x[i1]; break;
mass_water_bulk_x = base_mass_water_bulk_x; case MU:
x[i1]->master[0]->s->moles = base_moles_h2o; mu_x -= d2;
break; //k_temp(tc_x, patm_x);
case MH: gammas_pz(false);
s_eminus->la = base_x[i1]; break;
break; case GAS_MOLES:
case GAS_MOLES: if (gas_in == FALSE)
x[i1]->moles = base_x[i1];
for (size_t g = 0; g < base_gas_moles.size(); g++)
{
//phase_ptrs[g]->moles_x = base_gas_moles[g];
*phase_ptrs[g] = base_phases[g];
}
break;
case MU:
mu_x = base_x[i1];
gammas_pz(false);
break;
case PP:
continue; continue;
break; x[i]->moles -= d2;
case SS_MOLES: *use.Get_gas_phase_ptr() = base_gas_phase;
break; for (size_t g = 0; g < base_phases.size(); g++)
}
}
}
//molalities(TRUE);
//if (full_pitzer == TRUE)
// pitzer();
//mb_sums();
//residuals();
// Debugging
//output_msg("\n\nAfter jacobian_pz--------------------------------------------------------\n");
//print_all();
for (i = 0; i < count_unknowns; i++)
{
if (residual[i] != base[i])
{
if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 0.1)
{ {
//std::cerr << i << " " << residual[i] << " " << base[i] << std::endl; *phase_ptrs[g] = base_phases[g];
//output_flush();
} }
residual[i] = base[i]; break;
my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i]; case SS_MOLES:
delta[i] = -d2;
reset();
break;
} }
} }
molalities(TRUE);
if (full_pitzer == TRUE)
pitzer();
mb_sums();
residuals();
base.clear(); base.clear();
calculating_deriv = 0; calculating_deriv = 0;
return OK; return OK;

3
sit.cpp
View File

@ -870,7 +870,8 @@ Restart:
case GAS_MOLES: case GAS_MOLES:
if (gas_in == FALSE) if (gas_in == FALSE)
continue; continue;
d2 = d * 20 * x[i]->moles; d2 = (x[i]->moles > 1 ? 1 : 20);
d2 *= d * x[i]->moles;
if (d2 < 1e-14) if (d2 < 1e-14)
d2 = 1e-14; d2 = 1e-14;
x[i]->moles += d2; x[i]->moles += d2;