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

still produces different residuals

Browse Source
This commit is contained in:
David Parkhurst 2021-07-07 23:55:03 -06:00
parent 9022ded877
commit 71cf2a97b7
1 changed files with 130 additions and 56 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

186
pitzer.cpp
View File

@ -1947,24 +1947,87 @@ int Phreeqc::
jacobian_pz(void) jacobian_pz(void)
/* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */
{ // calculate the derivatives numerically { // calculate the derivatives numerically
std::vector<double> base; 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;
LDBLE d, d1, d2; LDBLE d, d1, d2;
int i, j; int i, j;
calculating_deriv = 1;
Restart: Restart:
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++)
{
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);
}
// Debugging
//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);
for (i = 0; i < count_unknowns; i++) base_x.resize(count_unknowns, 0.0);
for (size_t i1 = 0; i1 < count_unknowns; i1++)
{ {
base[i] = residual[i]; base[i1] = residual[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:
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 < base_gas_moles.size(); g++)
{
base_gas_moles[g] = phase_ptrs[g]->moles_x;
}
break;
case MU:
base_x[i1] = mu_x;
break;
case PP:
continue;
break;
case SS_MOLES:
break;
}
} }
d = 0.0001; d = 0.0001;
d1 = d * LOG_10; d1 = d * LOG_10;
@ -2075,55 +2138,50 @@ 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;
} }
switch (x[i]->type) for (size_t i1 = 0; i1 < count_unknowns; i1++)
{ {
case MB: switch (x[i1]->type)
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)
{ {
my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] = case MB:
exp(s_h2->lm * LOG_10) * 2; case ALK:
} case CB:
break; case SOLUTION_PHASE_BOUNDARY:
case PITZER_GAMMA: case EXCH:
x[i]->s->lg -= d; case SURFACE:
break; case SURFACE_CB:
case MH2O: case SURFACE_CB1:
//mass_water_aq_x /= (1 + d); case SURFACE_CB2:
//x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; case AH2O:
//break; x[i1]->master[0]->s->la = base_x[i1];
//DL_pitz break;
mass_water_aq_x -= d1; case PITZER_GAMMA:
if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL) x[i1]->s->lg = base_x[i1];
mass_water_bulk_x -= d1; break;
x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; case MH2O:
break; mass_water_aq_x = base_x[i1];
case MU: mass_water_bulk_x = base_mass_water_bulk_x;
mu_x -= d2; x[i1]->master[0]->s->moles = base_moles_h2o;
//k_temp(tc_x, patm_x); break;
gammas_pz(false); case MH:
break; s_eminus->la = base_x[i1];
case GAS_MOLES: break;
if (gas_in == FALSE) 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];
}
break;
case MU:
mu_x = base_x[i1];
gammas_pz(false);
break;
case PP:
continue; continue;
x[i]->moles -= d2; break;
break; case SS_MOLES:
case SS_MOLES: break;
delta[i] = -d2; }
reset();
break;
} }
} }
molalities(TRUE); molalities(TRUE);
@ -2131,6 +2189,22 @@ Restart:
pitzer(); pitzer();
mb_sums(); mb_sums();
residuals(); 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;
//output_flush();
}
residual[i] = base[i];
my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i];
}
}
base.clear(); base.clear();
calculating_deriv = 0; calculating_deriv = 0;
return OK; return OK;