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

revised jacobian_pz with new logic. Works with fixed_pressure examples H2S, H2S_pz, H2S_pz_appt, H2S_NaCl_Na2SO4.

Browse Source
This commit is contained in:
David Parkhurst 2021-07-10 16:37:37 -06:00
parent 71cf2a97b7
commit 8be1ba8327
2 changed files with 80 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
model.cpp
View File

@ -3227,12 +3227,12 @@ reset(void)
}
else if (x[i]->type == GAS_MOLES)
{
up = 1000. * x[i]->moles;
up = 10. * x[i]->moles;
if (up <= 0.0)
up = 1e-1;
if (up >= 1.0)
up = 1.;
down = x[i]->moles;
down = 0.3*x[i]->moles;
}
else if (x[i]->type == SS_MOLES)
{

96
pitzer.cpp
View File

@ -1949,26 +1949,29 @@ jacobian_pz(void)
{ // calculate the derivatives numerically
std::vector<double> base, base_x, base_gas_moles;
std::vector<class phase*> phase_ptrs;
double base_mass_water_bulk_x=0, base_moles_h2o=0;
std::vector<class phase> base_phases;
double base_mass_water_bulk_x = 0, base_moles_h2o = 0;
cxxGasPhase base_gas_phase;
LDBLE d, d1, d2;
int i, j;
Restart:
molalities(TRUE);
if (full_pitzer == TRUE)
pitzer();
mb_sums();
residuals();
//molalities(TRUE);
//if (full_pitzer == TRUE)
// pitzer();
//mb_sums();
//residuals();
if (use.Get_gas_phase_ptr() != NULL)
{
cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr();
for (size_t i = 0; i < gas_phase_ptr->Get_gas_comps().size(); i++)
for (size_t i = 0; i < use.Get_gas_phase_ptr()->Get_gas_comps().size(); 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);
phase_ptrs.push_back(phase_ptr);
}
base_gas_moles.resize(gas_phase_ptr->Get_gas_comps().size(), 0.0);
//base_gas_moles.resize(gas_phase_ptr->Get_gas_comps().size(), 0.0);
base_phases.resize(gas_phase_ptr->Get_gas_comps().size());
}
// Debugging
//output_msg("\n\nBefore jacobian_pz===================================================\n");
@ -1986,7 +1989,7 @@ Restart:
base_x.resize(count_unknowns, 0.0);
for (size_t i1 = 0; i1 < count_unknowns; i1++)
{
base[i1] = residual[i1];
//base[i1] = residual[i1];
switch (x[i1]->type)
{
case MB:
@ -2014,9 +2017,10 @@ Restart:
break;
case GAS_MOLES:
base_x[i1] = x[i1]->moles;
for (size_t g = 0; g < base_gas_moles.size(); g++)
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_gas_moles[g] = phase_ptrs[g]->moles_x;
base_phases[g] = *phase_ptrs[g];
}
break;
case MU:
@ -2029,7 +2033,62 @@ Restart:
break;
}
}
d = 0.0001;
//molalities(TRUE);
//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;
d2 = 0;
for (i = 0; i < count_unknowns; i++)
@ -2169,7 +2228,8 @@ Restart:
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]->moles_x = base_gas_moles[g];
*phase_ptrs[g] = base_phases[g];
}
break;
case MU:
@ -2184,11 +2244,11 @@ Restart:
}
}
}
molalities(TRUE);
if (full_pitzer == TRUE)
pitzer();
mb_sums();
residuals();
//molalities(TRUE);
//if (full_pitzer == TRUE)
// pitzer();
//mb_sums();
//residuals();
// Debugging
//output_msg("\n\nAfter jacobian_pz--------------------------------------------------------\n");
//print_all();