From df0d68b9d5cfb3df06c0a44118d0b8c39ecda91f Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Fri, 16 Jul 2021 22:27:03 -0600 Subject: [PATCH] Runs all the test cases. Numerical derivatives work, but still some changes in residuals before and after jacobian calculations. --- model.cpp | 49 +++++++------ pitzer.cpp | 11 ++- sit.cpp | 197 ++++++++++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 235 insertions(+), 22 deletions(-) diff --git a/model.cpp b/model.cpp index 2c691a75..88d33609 100644 --- a/model.cpp +++ b/model.cpp @@ -5304,6 +5304,7 @@ ss_f(LDBLE xb, LDBLE l_a0, LDBLE l_a1, LDBLE l_kc, LDBLE l_kb, LDBLE xcaq, f = xcaq * (xb / r + xc) + xbaq * (xb + r * xc) - 1; return (f); } +//#define ORIGINAL #ifdef ORIGINAL /* ---------------------------------------------------------------------- */ int Phreeqc:: @@ -5545,7 +5546,7 @@ numerical_jacobian(void) calculating_deriv = FALSE; return OK; } -#endif +#else /* ---------------------------------------------------------------------- */ int Phreeqc:: numerical_jacobian(void) @@ -5569,7 +5570,7 @@ numerical_jacobian(void) return(OK); calculating_deriv = TRUE; - jacobian_sums(); + //jacobian_sums(); if (use.Get_gas_phase_ptr() != NULL) { //cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); @@ -5583,7 +5584,7 @@ numerical_jacobian(void) base_phases[i] = *phase_ptr; } } - //gammas(mu_x); + gammas(mu_x); molalities(TRUE); mb_sums(); //mb_gases(); @@ -5652,7 +5653,7 @@ numerical_jacobian(void) case MU: d2 = d * mu_x; mu_x += d2; - gammas(mu_x); + //gammas(mu_x); break; case PP: for (j = 0; j < count_unknowns; j++) @@ -5695,7 +5696,7 @@ numerical_jacobian(void) x[i]->moles += d2; break; } - //gammas(mu_x); + gammas(mu_x); molalities(TRUE); mb_sums(); //mb_gases(); @@ -5746,7 +5747,7 @@ numerical_jacobian(void) break; case MU: mu_x -= d2; - gammas(mu_x); + //gammas(mu_x); break; case PP: delta[i] = -d2; @@ -5768,25 +5769,33 @@ numerical_jacobian(void) *phase_ptrs[g] = base_phases[g]; } } - gammas(mu_x); - molalities(TRUE); - mb_sums(); - //mb_gases(); - //mb_ss(); - residuals(); - //for (size_t i2 = 0; i2 < count_unknowns; i2++) - //{ - // if (fabs(2.0 * (residual[i2] - base[i2]) / (residual[i2] + base[i2])) > 1e-6) - // { - // //std::cerr << "iter: " << iterations << ". " << std::endl; - // std::cerr << i2 << ": " << x[i2]->description << " " << residual[i2] << " " << base[i2] << std::endl; - // } - //} + //gammas(mu_x); + //molalities(TRUE); + //mb_sums(); + ////mb_gases(); + ////mb_ss(); + //residuals(); } + gammas(mu_x); + molalities(TRUE); + mb_sums(); + //mb_gases(); + //mb_ss(); + residuals(); + //for (i = 0; i < count_unknowns; i++) + //{ + // //Debugging + // if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 && + // fabs(residual[i]) + fabs(base[i]) > 1e-8) + // { + // std::cerr << i << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl; + // } + //} base.clear(); calculating_deriv = FALSE; return OK; } +#endif /* ---------------------------------------------------------------------- */ void Phreeqc:: set_inert_moles(void) diff --git a/pitzer.cpp b/pitzer.cpp index 81f5f098..21d667aa 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1954,6 +1954,7 @@ jacobian_pz(void) LDBLE d, d1, d2; int i, j; Restart: + calculating_deriv = 1; molalities(TRUE); if (full_pitzer == TRUE) { @@ -1961,7 +1962,6 @@ Restart: } mb_sums(); residuals(); - calculating_deriv = 1; if (use.Get_gas_phase_ptr() != NULL) { cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); @@ -2155,6 +2155,15 @@ Restart: mb_sums(); residuals(); } + //for (i = 0; i < count_unknowns; i++) + //{ + // //Debugging + // if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 && + // fabs(residual[i]) + fabs(base[i]) > 1e-8) + // { + // std::cerr << i << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl; + // } + //} base.clear(); calculating_deriv = 0; return OK; diff --git a/sit.cpp b/sit.cpp index f410b773..a1b140b6 100644 --- a/sit.cpp +++ b/sit.cpp @@ -798,7 +798,8 @@ sit_revise_guesses(void) /*gammas(mu_x); */ return (OK); } - +//#define ORIGINAL +#ifdef ORIGINAL /* ---------------------------------------------------------------------- */ int Phreeqc:: jacobian_sit(void) @@ -952,7 +953,201 @@ Restart: residuals(); return OK; } +#else +/* ---------------------------------------------------------------------- */ +int Phreeqc:: +jacobian_sit(void) +/* ---------------------------------------------------------------------- */ +{ + std::vector base; + LDBLE d, d1, d2; + int i, j; + cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); + std::vector phase_ptrs; + std::vector base_phases; + double base_mass_water_bulk_x = 0, base_moles_h2o = 0; + cxxGasPhase base_gas_phase; +Restart: + if (use.Get_gas_phase_ptr() != NULL) + { + cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); + 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]); + class phase* phase_ptr = phase_bsearch(gas_comp_ptr->Get_phase_name().c_str(), &j, FALSE); + phase_ptrs.push_back(phase_ptr); + base_phases[i] = *phase_ptr; + } + } + size_t pz_max_unknowns = max_unknowns; + //k_temp(tc_x, patm_x); + calculating_deriv = 1; + molalities(TRUE); + if (full_pitzer == TRUE) + { + + sit(); + } + mb_sums(); + residuals(); + base = residual; // std::vectors + d = 0.0001; + d1 = d * LOG_10; + d2 = 0; + for (i = 0; i < count_unknowns; i++) + { + switch (x[i]->type) + { + case MB: + case ALK: + case CB: + case SOLUTION_PHASE_BOUNDARY: + case EXCH: + case SURFACE: + case SURFACE_CB: + case SURFACE_CB1: + case SURFACE_CB2: + x[i]->master[0]->s->la += d; + d2 = d1; + break; + case AH2O: + x[i]->master[0]->s->la += d; + d2 = d1; + break; + case PITZER_GAMMA: + if (!full_pitzer) + continue; + x[i]->s->lg += d; + d2 = d; + break; + case MH2O: + mass_water_aq_x *= (1.0 + d); + x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; + d2 = log(1.0 + d); + break; + case MH: + s_eminus->la += d; + d2 = d1; + break; + /* + if (pitzer_pe == TRUE) + { + s_eminus->la += d; + d2 = d1; + break; + } + else + { + continue; + } + */ + case GAS_MOLES: + if (gas_in == FALSE) + continue; + d2 = (x[i]->moles > 1 ? 1 : 20); + d2 *= d * x[i]->moles; + if (d2 < 1e-14) + d2 = 1e-14; + x[i]->moles += d2; + break; + case MU: + //continue; + d2 = d * mu_x; + mu_x += d2; + //k_temp(tc_x, patm_x); + gammas_sit(); + break; + case PP: + case SS_MOLES: + continue; + break; + } + molalities(TRUE); + if (max_unknowns > pz_max_unknowns) + { + gammas_sit(); + jacobian_sums(); + goto Restart; + } + if (full_pitzer == TRUE) + sit(); + mb_sums(); + residuals(); + for (j = 0; j < count_unknowns; j++) + { + my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] = + -(residual[j] - base[j]) / d2; + } + switch (x[i]->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) + { + my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] = + exp(s_h2->lm * LOG_10) * 2; + } + break; + case PITZER_GAMMA: + x[i]->s->lg -= d; + break; + case MH2O: + mass_water_aq_x /= (1 + d); + x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; + break; + case MU: + mu_x -= d2; + //k_temp(tc_x, patm_x); + gammas_sit(); + break; + case GAS_MOLES: + if (gas_in == FALSE) + continue; + x[i]->moles -= d2; + break; + } + if (use.Get_gas_phase_ptr() != NULL) + { + *use.Get_gas_phase_ptr() = base_gas_phase; + for (size_t g = 0; g < base_phases.size(); g++) + { + *phase_ptrs[g] = base_phases[g]; + } + } + } + molalities(TRUE); + if (full_pitzer == TRUE) + sit(); + mb_sums(); + residuals(); + //for (i = 0; i < count_unknowns; i++) + //{ + // //Debugging + // if (fabs(2.0 * (residual[i] - base[i]) / (residual[i] + base[i])) > 1e-2 && + // fabs(residual[i]) + fabs(base[i]) > 1e-8) + // { + // std::cerr << i << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl; + // } + //} + calculating_deriv = 0; + return OK; +} +#endif /* ---------------------------------------------------------------------- */ int Phreeqc:: model_sit(void)