From 9022ded8773669859d20745bf64776cd0c0741f5 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Tue, 6 Jul 2021 21:52:06 -0600 Subject: [PATCH 01/11] Tony H2S. Amm.dat, phreeqc.dat, pitzer.dat, utf8, updated test cases --- Solution.cxx | 18 ++++++++++++++---- class_main.cpp | 14 ++++++++------ gases.cpp | 16 ++++++++-------- global_structures.h | 1 + inverse.cpp | 4 ++-- kinetics.cpp | 4 ++-- model.cpp | 7 +++---- pitzer.cpp | 7 ++++--- prep.cpp | 34 +++++++++++++--------------------- print.cpp | 4 ++-- read.cpp | 14 +++++++------- sit.cpp | 2 +- step.cpp | 12 ++++++++++-- transport.cpp | 23 +++++++++++++++++++---- 14 files changed, 94 insertions(+), 66 deletions(-) diff --git a/Solution.cxx b/Solution.cxx index d29a788f..85055223 100644 --- a/Solution.cxx +++ b/Solution.cxx @@ -114,6 +114,17 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions, this->n_user = this->n_user_end = l_n_user; this->new_def = false; this->ah2o = 0; + // potV is an external variable, imposed in a given solution, not mixed. + std::map < int, cxxSolution >::const_iterator sol = solutions.find(mix.Get_n_user()); + const cxxSolution *cxxsoln_ptr1; + if (sol != solutions.end()) + { + cxxsoln_ptr1 = &(sol->second); + if (cxxsoln_ptr1->new_def) + this->potV = 0.0; + else + this->potV = cxxsoln_ptr1->potV; + } // // Mix solutions // @@ -121,8 +132,7 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions, std::map < int, LDBLE >::const_iterator it; for (it = mixcomps.begin(); it != mixcomps.end(); it++) { - std::map < int, cxxSolution >::const_iterator sol = - solutions.find(it->first); + sol = solutions.find(it->first); if (sol == solutions.end()) { std::ostringstream msg; @@ -131,7 +141,7 @@ cxxSolution::cxxSolution(std::map < int, cxxSolution > &solutions, } else { - const cxxSolution *cxxsoln_ptr1 = &(sol->second); + cxxsoln_ptr1 = &(sol->second); this->add(*cxxsoln_ptr1, it->second); } } @@ -1388,7 +1398,7 @@ cxxSolution::add(const cxxSolution & addee, LDBLE extensive) this->cb += addee.cb * extensive; this->density = f1 * this->density + f2 * addee.density; this->patm = f1 * this->patm + f2 * addee.patm; - this->potV = f1 * this->potV + f2 * addee.potV; + // this->potV = f1 * this->potV + f2 * addee.potV; // appt this->mass_water += addee.mass_water * extensive; this->soln_vol += addee.soln_vol * extensive; this->total_alkalinity += addee.total_alkalinity * extensive; diff --git a/class_main.cpp b/class_main.cpp index 9e6a0d0c..ee43e385 100644 --- a/class_main.cpp +++ b/class_main.cpp @@ -147,8 +147,8 @@ main_method(int argc, char *argv[]) } Phreeqc MyCopy = *this; this->clean_up(); - //this->init(); - //this->initialize(); + // this->init(); + // this->initialize(); /* * Read input data for simulation */ @@ -167,7 +167,7 @@ main_method(int argc, char *argv[]) * Display successful status */ pr.headings = TRUE; - //errors = do_status(); + // errors = do_status(); if (errors != 0) { return errors; @@ -292,7 +292,7 @@ write_banner(void) /* version */ #ifdef NPP - len = sprintf(buffer, "* PHREEQC-%s *", "3.6.5"); + len = sprintf(buffer, "* PHREEQC-%s *", "3.7.1"); #else len = sprintf(buffer, "* PHREEQC-%s *", "@VERSION@"); #endif @@ -316,7 +316,7 @@ write_banner(void) /* date */ #ifdef NPP - len = sprintf(buffer, "%s", "February 24, 2021"); + len = sprintf(buffer, "%s", "July 5, 2021"); #else len = sprintf(buffer, "%s", "@VER_DATE@"); #endif @@ -329,6 +329,7 @@ write_banner(void) return 0; } + /* ---------------------------------------------------------------------- */ int Phreeqc:: process_file_names(int argc, char *argv[], std::istream **db_cookie, @@ -499,7 +500,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie, output_msg(sformatf(" Input file: %s\n", in_file.c_str())); output_msg(sformatf(" Output file: %s\n", out_file.c_str())); #ifdef NPP - output_msg(sformatf("Using PHREEQC: version 3.6.5, compiled February 24, 2021\n")); + output_msg(sformatf("Using PHREEQC: version 3.7.1, compiled July 5, 2021\n")); #endif output_msg(sformatf("Database file: %s\n\n", token.c_str())); #ifdef NPP @@ -508,6 +509,7 @@ process_file_names(int argc, char *argv[], std::istream **db_cookie, /* * local cleanup */ + line = (char *) free_check_null(line); line_save = (char *) free_check_null(line_save); diff --git a/gases.cpp b/gases.cpp index e2c5495c..dec63178 100644 --- a/gases.cpp +++ b/gases.cpp @@ -436,22 +436,22 @@ calc_PR(void) { if (!strcmp(phase_ptr1->name, "CO2(g)")) a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 - else if (!strcmp(phase_ptr1->name, "H2S(g)")) + else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) a_aa *= 0.81; - else if (!strcmp(phase_ptr1->name, "CH4(g)")) + else if (!strcmp(phase_ptr1->name, "CH4(g)") || !strcmp(phase_ptr1->name, "Mtg(g)")) a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "N2(g)")) + else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) a_aa *= 0.51; } if (!strcmp(phase_ptr1->name, "H2O(g)")) { if (!strcmp(phase_ptr->name, "CO2(g)")) a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "H2S(g)")) + else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "CH4(g)")) + else if (!strcmp(phase_ptr->name, "CH4(g)") || !strcmp(phase_ptr->name, "Mtg(g)")) a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "N2(g)")) + else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) a_aa *= 0.51; } a_aa_sum += phase_ptr->fraction_x * phase_ptr1->fraction_x * a_aa; @@ -598,10 +598,10 @@ calc_PR(void) { phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) * log((rz + 2.41421356 * B) / (rz - 0.41421356 * B)); - phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi)); + //phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi)); } else - phi = -3.0; // fugacity coefficient > 0.05 + phi = -3.0; // fugacity coefficient = 0.05 phase_ptr->pr_phi = exp(phi); phase_ptr->pr_si_f = phi / LOG_10; // pr_si_f updated // **** diff --git a/global_structures.h b/global_structures.h index e66b374c..b217d965 100644 --- a/global_structures.h +++ b/global_structures.h @@ -41,6 +41,7 @@ #define DISP 2 #define STAG 3 #define NOMIX 4 +#define MIX_BS 5 // mix boundary solutions in electromigration #define CONVERGED 2 #define MASS_BALANCE 3 diff --git a/inverse.cpp b/inverse.cpp index 7ce99fbb..81a33781 100644 --- a/inverse.cpp +++ b/inverse.cpp @@ -272,7 +272,7 @@ setup_inverse(class inverse *inv_ptr) /* * Define inv_zero and inv_zero array, delta */ - for (i = 0; i < max; i++) + for (i = 0; i < (int) max; i++) inv_zero[i] = 0.0; memcpy((void *) &(delta[0]), (void *) &(inv_zero[0]), @@ -3595,7 +3595,7 @@ check_isotopes(class inverse *inv_ptr) /* use solution-defined uncertainties second */ } #ifdef NPP - else if (inv_ptr->i_u[i].count_uncertainties > 0 + else if (inv_ptr->i_u[i].uncertainties.size() > 0 && !isnan(inv_ptr->i_u[i].uncertainties[inv_ptr->i_u[i].uncertainties.size() - 1])) #else else if (inv_ptr->i_u[i].uncertainties.size() > 0 diff --git a/kinetics.cpp b/kinetics.cpp index dab9ab62..27ce4bdc 100644 --- a/kinetics.cpp +++ b/kinetics.cpp @@ -1595,7 +1595,7 @@ set_transport(int i, int use_mix, int use_kinetics, int nsaver) /* * i --user number for soln, reaction, etc. * use_mix --integer flag - state == TRANSPORT: DISP, STAG, NOMIX + state == TRANSPORT: DISP, STAG, NOMIX, MIX_BS state == REACTION: TRUE, FALSE * use_kinetics --true or false flag to calculate kinetic reactions * nsaver --user number to store solution @@ -1615,7 +1615,7 @@ set_transport(int i, int use_mix, int use_kinetics, int nsaver) use.Set_n_mix_user(i); use.Set_n_mix_user_orig(i); } - else if (use_mix == STAG && multi_Dflag != TRUE) + else if ((use_mix == STAG && multi_Dflag != TRUE) || use_mix == MIX_BS) { use.Set_mix_ptr(Utilities::Rxn_find(Rxn_mix_map, i)); if (use.Get_mix_ptr() != NULL) diff --git a/model.cpp b/model.cpp index cebc94e8..56104b73 100644 --- a/model.cpp +++ b/model.cpp @@ -2649,8 +2649,8 @@ calc_gas_pressures(void) lp += rxn_ptr->s->la * rxn_ptr->coef; } phase_ptr->p_soln_x = exp(LOG_10 * (lp - phase_ptr->pr_si_f)); - if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x > 90) - phase_ptr->p_soln_x = 90; + //if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x > 90) + // phase_ptr->p_soln_x = 90; if (gas_phase_ptr->Get_type() == cxxGasPhase::GP_PRESSURE) { @@ -5425,8 +5425,7 @@ numerical_jacobian(void) case GAS_MOLES: if (gas_in == FALSE) continue; - - d2 = d * x[i]->moles; + d2 = d * 20 * x[i]->moles; if (d2 < 1e-14) d2 = 1e-14; x[i]->moles += d2; diff --git a/pitzer.cpp b/pitzer.cpp index 59882063..2f4e8e3f 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -2026,9 +2026,10 @@ Restart: case GAS_MOLES: if (gas_in == FALSE) continue; - d2 = d * x[i]->moles; - if (d2 < 1e-14) - d2 = 1e-14; + d2 = d * 30 * x[i]->moles; + d2 = (d2 < ineq_tol ? ineq_tol : d2); + //if (d2 < 1e-14) + // d2 = 1e-14; x[i]->moles += d2; break; case MU: diff --git a/prep.cpp b/prep.cpp index 5791517b..82b3125c 100644 --- a/prep.cpp +++ b/prep.cpp @@ -3841,15 +3841,11 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) { if (!strcmp(phase_ptr1->name, "CO2(g)")) a_aa *= 0.81; // Soreide and Whitson, 1992, FPE 77, 217 - else if (!strcmp(phase_ptr1->name, "H2S(g)")) + else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) a_aa *= 0.81; - else if (!strcmp(phase_ptr1->name, "CH4(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; - else if (!strcmp(phase_ptr1->name, "Mtg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr1->name, "N2(g)")) + else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) a_aa *= 0.51; else if (!strcmp(phase_ptr1->name, "Ethane(g)")) a_aa *= 0.51; @@ -3860,15 +3856,11 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) { if (!strcmp(phase_ptr->name, "CO2(g)")) a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "H2S(g)")) + else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) a_aa *= 0.81; - else if (!strcmp(phase_ptr->name, "CH4(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; - else if (!strcmp(phase_ptr->name, "Mtg(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "Methane(g)")) - a_aa *= 0.51; - else if (!strcmp(phase_ptr->name, "N2(g)")) + else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) a_aa *= 0.51; else if (!strcmp(phase_ptr->name, "Ethane(g)")) a_aa *= 0.51; @@ -4002,17 +3994,17 @@ calc_PR(std::vector phase_ptrs, LDBLE P, LDBLE TK, LDBLE V_m) { phi = B_r * (rz - 1) - log(rz - B) + A / (2.828427 * B) * (B_r - 2.0 * phase_ptr->pr_aa_sum2 / a_aa_sum) * log((rz + 2.41421356 * B) / (rz - 0.41421356 * B)); - phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi)); + //phi = (phi > 4.44 ? 4.44 : (phi < -3 ? -3 : phi)); //if (phi > 4.44) // phi = 4.44; } else - phi = -3.0; // fugacity coefficient > 0.05 - if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3) - { - // avoid such phi... - phi = -3; - } + phi = -3.0; // fugacity coefficient = 0.05 + //if (/*!strcmp(phase_ptr->name, "H2O(g)") && */phi < -3) + //{ + //// avoid such phi... + // phi = -3; + //} phase_ptr->pr_phi = exp(phi); phase_ptr->pr_si_f = phi / LOG_10; // for initial equilibrations, adapt log_k of the gas phase... diff --git a/print.cpp b/print.cpp index cb01ea9f..d4ad7b79 100644 --- a/print.cpp +++ b/print.cpp @@ -708,8 +708,8 @@ print_gas_phase(void) (double) initial_moles, (double) moles, (double) delta_moles)); - if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x == 90) - output_msg(" WARNING: The pressure of H2O(g) is above the program limit: use the polynomial for log_k.\n"); + //if (!strcmp(phase_ptr->name, "H2O(g)") && phase_ptr->p_soln_x == 90) + // output_msg(" WARNING: The pressure of H2O(g) is fixed to the program limit.\n"); } output_msg("\n"); diff --git a/read.cpp b/read.cpp index fb9e78e1..4ade4fe4 100644 --- a/read.cpp +++ b/read.cpp @@ -2619,7 +2619,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v) /* * Read supcrt parms and Ionic strength terms */ - for (j = 0; j < 9; j++) + for (j = 0; j < 11; j++) { delta_v[j] = 0.0; } @@ -2627,7 +2627,7 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v) /* Vmax, dmax... delta_v[10] = 999.0; delta_v[11] = 1.0; */ - j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT /*SCANFORMAT SCANFORMAT */, + j = sscanf(cptr, SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT SCANFORMAT , /* a1..a4 */ &(delta_v[0]), &(delta_v[1]), &(delta_v[2]), &(delta_v[3]), /* wref */ @@ -2635,9 +2635,9 @@ read_aq_species_vm_parms(const char* cptr, LDBLE * delta_v) /* b_Av */ &(delta_v[5]), /* c1..c4 */ - &(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9])); - /* vmax, dmax - &(delta_v[10]), &(delta_v[11])); */ + &(delta_v[6]), &(delta_v[7]), &(delta_v[8]), &(delta_v[9]), + //vmP, exP + &(delta_v[10]), &(delta_v[11])); if (j < 1) { input_error++; @@ -2665,8 +2665,8 @@ read_vm_only(const char* cptr, LDBLE * delta_v, DELTA_V_UNIT * units) char token[MAX_LENGTH]; /* * Read analytical expression - */ - for (j = 0; j < 8; j++) + */ + for (j = 0; j < 9; j++) { delta_v[j] = 0.0; } diff --git a/sit.cpp b/sit.cpp index 27825c26..236793d3 100644 --- a/sit.cpp +++ b/sit.cpp @@ -870,7 +870,7 @@ Restart: case GAS_MOLES: if (gas_in == FALSE) continue; - d2 = d * x[i]->moles; + d2 = d * 20 * x[i]->moles; if (d2 < 1e-14) d2 = 1e-14; x[i]->moles += d2; diff --git a/step.cpp b/step.cpp index b7204ecd..a5503f02 100644 --- a/step.cpp +++ b/step.cpp @@ -57,11 +57,19 @@ step(LDBLE step_fraction) */ if (use.Get_mix_ptr() != NULL) { - add_mix(use.Get_mix_ptr()); + add_mix(use.Get_mix_ptr()); + int n = use.Get_n_mix_user_orig(); + if (n == 0 || n == count_cells + 1) + { + cxxSolution *solution_ptr = Utilities::Rxn_find(Rxn_solution_map, n); + if (solution_ptr != NULL && !solution_ptr->Get_new_def()) + potV_x = solution_ptr->Get_potV(); + } } else if (use.Get_solution_ptr() != NULL) { add_solution(use.Get_solution_ptr(), 1.0, 1.0); + potV_x = use.Get_solution_ptr()->Get_potV(); cell_no = use.Get_n_solution_user(); } else @@ -367,7 +375,7 @@ add_solution(cxxSolution *solution_ptr, LDBLE extensive, LDBLE intensive) tc_x += solution_ptr->Get_tc() * intensive; ph_x += solution_ptr->Get_ph() * intensive; patm_x += solution_ptr->Get_patm() * intensive; - potV_x += solution_ptr->Get_potV() * intensive; + //potV_x += solution_ptr->Get_potV() * intensive; solution_pe_x += solution_ptr->Get_pe() * intensive; mu_x += solution_ptr->Get_mu() * intensive; ah2o_x += solution_ptr->Get_ah2o() * intensive; diff --git a/transport.cpp b/transport.cpp index ca2317c0..00fedd8a 100644 --- a/transport.cpp +++ b/transport.cpp @@ -1,4 +1,4 @@ -#include "Utils.h" +#include "Utils.h" #include "Phreeqc.h" #include "phqalloc.h" #include "Exchange.h" @@ -605,7 +605,10 @@ transport(void) if (i == 0 || i == count_cells + 1) { - run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i + if (dV_dcell) + run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i + else + run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i } else { @@ -868,7 +871,12 @@ transport(void) status(0, token); if (i == 0 || i == count_cells + 1) - run_reactions(i, kin_time, NOMIX, step_fraction); + { + if (dV_dcell) + run_reactions(i, kin_time, MIX_BS, step_fraction); // nsaver = i + else + run_reactions(i, kin_time, NOMIX, step_fraction); // nsaver = i + } else { run_reactions(i, kin_time, DISP, step_fraction); @@ -1070,6 +1078,12 @@ print_punch(int i, boolean active) if (!active) run_reactions(i, 0, NOMIX, 0); cell_no = i; + if (dV_dcell) + { + use.Set_n_solution_user(i); + use.Get_solution_ptr()->Set_potV(cell_data[i].potV); + potV_x = cell_data[i].potV; + } use.Set_kinetics_ptr(Utilities::Rxn_find(Rxn_kinetics_map, i)); if (use.Get_kinetics_ptr() != NULL) { @@ -1413,7 +1427,7 @@ mix_stag(int i, LDBLE kin_time, int l_punch, LDBLE step_fraction) for (std::map < int, LDBLE >::const_iterator it = use.Get_mix_ptr()->Get_mixComps().begin(); it != use.Get_mix_ptr()->Get_mixComps().end(); it++) { - if (it->first > i && it->first < all_cells) + if (it->first > i && it->first < all_cells && it->first != count_cells + 1) { k = it->first; ptr_imm = Utilities::Rxn_find(Rxn_solution_map, k); @@ -6053,6 +6067,7 @@ calc_vm_Cl(void) V_Cl = s_ptr->logk[vma1] + s_ptr->logk[vma2] / pb_s + (s_ptr->logk[vma3] + s_ptr->logk[vma4] / pb_s) / TK_s - s_ptr->logk[wref] * QBrn; + /* the ionic strength term * I^0.5... */ if (s_ptr->logk[b_Av] < 1e-5) V_Cl += s_ptr->z * s_ptr->z * 0.5 * DH_Av * sqrt_mu; From 71cf2a97b782edefc376fbd62cde3442cbc33b4b Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Wed, 7 Jul 2021 23:55:03 -0600 Subject: [PATCH 02/11] still produces different residuals --- pitzer.cpp | 186 +++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 130 insertions(+), 56 deletions(-) diff --git a/pitzer.cpp b/pitzer.cpp index 2f4e8e3f..43ab52bc 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1947,24 +1947,87 @@ int Phreeqc:: jacobian_pz(void) /* ---------------------------------------------------------------------- */ { // calculate the derivatives numerically - std::vector base; + std::vector base, base_x, base_gas_moles; + std::vector phase_ptrs; + double base_mass_water_bulk_x=0, base_moles_h2o=0; LDBLE d, d1, d2; int i, j; - calculating_deriv = 1; 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; //k_temp(tc_x, patm_x); - if (full_pitzer == TRUE) - { - molalities(TRUE); - pitzer(); - residuals(); - } + //if (full_pitzer == TRUE) + //{ + // molalities(TRUE); + // pitzer(); + // residuals(); + //} 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; d1 = d * LOG_10; @@ -2075,55 +2138,50 @@ Restart: if (x[i]->type == MH2O) // DL_pitz 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: - 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) + switch (x[i1]->type) { - 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; - //DL_pitz - mass_water_aq_x -= d1; - if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL) - mass_water_bulk_x -= d1; - 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_pz(false); - break; - case GAS_MOLES: - if (gas_in == FALSE) + 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]; + } + break; + case MU: + mu_x = base_x[i1]; + gammas_pz(false); + break; + case PP: continue; - x[i]->moles -= d2; - break; - case SS_MOLES: - delta[i] = -d2; - reset(); - break; + break; + case SS_MOLES: + break; + } } } molalities(TRUE); @@ -2131,6 +2189,22 @@ Restart: 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; + //output_flush(); + } + residual[i] = base[i]; + my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i]; + } + } base.clear(); calculating_deriv = 0; return OK; From 8be1ba8327a23d9c755e25df1aad0de008b70777 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Sat, 10 Jul 2021 16:37:37 -0600 Subject: [PATCH 03/11] revised jacobian_pz with new logic. Works with fixed_pressure examples H2S, H2S_pz, H2S_pz_appt, H2S_NaCl_Na2SO4. --- model.cpp | 4 +-- pitzer.cpp | 96 ++++++++++++++++++++++++++++++++++++++++++++---------- 2 files changed, 80 insertions(+), 20 deletions(-) diff --git a/model.cpp b/model.cpp index 56104b73..d31ff093 100644 --- a/model.cpp +++ b/model.cpp @@ -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) { diff --git a/pitzer.cpp b/pitzer.cpp index 43ab52bc..cdf91324 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1949,26 +1949,29 @@ jacobian_pz(void) { // calculate the derivatives numerically std::vector base, base_x, base_gas_moles; std::vector phase_ptrs; - double base_mass_water_bulk_x=0, base_moles_h2o=0; + std::vector 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(); From 13ec2fcd3ea5a10b09010a3b8ceda6708c25ddf8 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Tue, 13 Jul 2021 14:24:00 -0600 Subject: [PATCH 04/11] best I could do for H2S while maintaining old tests. Used INCREMENTAL reactions --- gases.cpp | 16 +++- model.cpp | 12 +-- pitzer.cpp | 257 +++++++++++++++-------------------------------------- sit.cpp | 3 +- 4 files changed, 92 insertions(+), 196 deletions(-) diff --git a/gases.cpp b/gases.cpp index dec63178..d64d2c97 100644 --- a/gases.cpp +++ b/gases.cpp @@ -101,7 +101,7 @@ build_fixed_volume_gas(void) } /* All elements in gas */ - for (j = 0; j < count_elts; j++) + for (j = 0; j < (int) count_elts; j++) { unknown_ptr = NULL; 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", phase_ptr->name)); } - for (j = 0; j < count_elts; j++) + for (j = 0; j < (int) count_elts; j++) { unknown_ptr = NULL; 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 else if (!strcmp(phase_ptr1->name, "H2S(g)") || !strcmp(phase_ptr1->name, "H2Sg(g)")) 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; else if (!strcmp(phase_ptr1->name, "N2(g)") || !strcmp(phase_ptr1->name, "Ntg(g)")) 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)")) { @@ -449,10 +453,14 @@ calc_PR(void) a_aa *= 0.81; else if (!strcmp(phase_ptr->name, "H2S(g)") || !strcmp(phase_ptr->name, "H2Sg(g)")) 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; else if (!strcmp(phase_ptr->name, "N2(g)") || !strcmp(phase_ptr->name, "Ntg(g)")) 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_sum2 += phase_ptr1->fraction_x * a_aa; diff --git a/model.cpp b/model.cpp index d31ff093..adeea063 100644 --- a/model.cpp +++ b/model.cpp @@ -2613,13 +2613,12 @@ calc_gas_pressures(void) V_m = (2. * gas_phase_ptr->Get_v_m() + V_m) / 3; else 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 */ numerical_fixed_volume = true; //switch_numerical = true; - warning_msg - ("Numerical method failed, switching to numerical derivatives."); + if (!pitzer_model) warning_msg("Numerical method failed, switching to numerical derivatives."); prep(); //switch_numerical = false; } @@ -3227,12 +3226,12 @@ reset(void) } else if (x[i]->type == GAS_MOLES) { - up = 10. * x[i]->moles; + up = 1000. * x[i]->moles; if (up <= 0.0) up = 1e-1; if (up >= 1.0) up = 1.; - down = 0.3*x[i]->moles; + down = x[i]->moles; } else if (x[i]->type == SS_MOLES) { @@ -5425,7 +5424,8 @@ numerical_jacobian(void) case GAS_MOLES: if (gas_in == FALSE) continue; - d2 = d * 20 * x[i]->moles; + d2 = (x[i]->moles > 1 ? 1 : 20); + d2 *= d * x[i]->moles; if (d2 < 1e-14) d2 = 1e-14; x[i]->moles += d2; diff --git a/pitzer.cpp b/pitzer.cpp index cdf91324..ce7576ba 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1955,140 +1955,35 @@ jacobian_pz(void) LDBLE d, d1, d2; int i, j; -Restart: - //molalities(TRUE); - //if (full_pitzer == TRUE) - // pitzer(); - //mb_sums(); - //residuals(); + calculating_deriv = 1; if (use.Get_gas_phase_ptr() != NULL) { 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]); 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; } - //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"); - //print_all(); - calculating_deriv = 1; +Restart: size_t pz_max_unknowns = max_unknowns; //k_temp(tc_x, patm_x); - //if (full_pitzer == TRUE) - //{ - // molalities(TRUE); - // pitzer(); - // residuals(); - //} + if (full_pitzer == TRUE) + { + molalities(TRUE); + pitzer(); + residuals(); + } base.resize(count_unknowns); - base_x.resize(count_unknowns, 0.0); - for (size_t i1 = 0; i1 < count_unknowns; i1++) + for (i = 0; i < count_unknowns; 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 < 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; - } + base[i] = residual[i]; } - //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; + d = 0.0001; d1 = d * LOG_10; d2 = 0; for (i = 0; i < count_unknowns; i++) @@ -2148,10 +2043,9 @@ Restart: case GAS_MOLES: if (gas_in == FALSE) continue; - d2 = d * 30 * x[i]->moles; - d2 = (d2 < ineq_tol ? ineq_tol : d2); - //if (d2 < 1e-14) - // d2 = 1e-14; + d2 = d * x[i]->moles; + if (d2 < 1e-14) + d2 = 1e-14; x[i]->moles += d2; break; case MU: @@ -2197,74 +2091,67 @@ Restart: if (x[i]->type == MH2O) // DL_pitz 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: - 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: + 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; + //DL_pitz + mass_water_aq_x -= d1; + if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL) + mass_water_bulk_x -= d1; + 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_pz(false); + break; + case GAS_MOLES: + if (gas_in == FALSE) continue; - break; - case SS_MOLES: - break; - } - } - } - //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) + x[i]->moles -= d2; + *use.Get_gas_phase_ptr() = base_gas_phase; + for (size_t g = 0; g < base_phases.size(); g++) { - //std::cerr << i << " " << residual[i] << " " << base[i] << std::endl; - //output_flush(); + *phase_ptrs[g] = base_phases[g]; } - residual[i] = base[i]; - my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i]; + break; + case SS_MOLES: + delta[i] = -d2; + reset(); + break; } } + molalities(TRUE); + if (full_pitzer == TRUE) + pitzer(); + mb_sums(); + residuals(); base.clear(); calculating_deriv = 0; return OK; diff --git a/sit.cpp b/sit.cpp index 236793d3..f410b773 100644 --- a/sit.cpp +++ b/sit.cpp @@ -870,7 +870,8 @@ Restart: case GAS_MOLES: if (gas_in == FALSE) continue; - d2 = d * 20 * x[i]->moles; + d2 = (x[i]->moles > 1 ? 1 : 20); + d2 *= d * x[i]->moles; if (d2 < 1e-14) d2 = 1e-14; x[i]->moles += d2; From 20281a0e7aaf3903dd0a0c3976940747e9b20b0e Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Wed, 14 Jul 2021 13:53:28 -0600 Subject: [PATCH 05/11] always reset gases --- pitzer.cpp | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/pitzer.cpp b/pitzer.cpp index ce7576ba..abfa14b7 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -2135,17 +2135,20 @@ Restart: if (gas_in == FALSE) continue; x[i]->moles -= d2; - *use.Get_gas_phase_ptr() = base_gas_phase; - for (size_t g = 0; g < base_phases.size(); g++) - { - *phase_ptrs[g] = base_phases[g]; - } break; case SS_MOLES: delta[i] = -d2; reset(); 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) From aef51fa3e6f5edbf2448ffda3aaa73bcb703993a Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Thu, 15 Jul 2021 11:38:17 -0600 Subject: [PATCH 06/11] Finally have derivatives right, I think --- pitzer.cpp | 50 ++++++++++++++++++++++++++++++++++---------------- 1 file changed, 34 insertions(+), 16 deletions(-) diff --git a/pitzer.cpp b/pitzer.cpp index abfa14b7..ea5c2bd2 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1954,7 +1954,14 @@ jacobian_pz(void) cxxGasPhase base_gas_phase; LDBLE d, d1, d2; int i, j; - +Restart: + molalities(TRUE); + if (full_pitzer == TRUE) + { + pitzer(); + } + mb_sums(); + residuals(); calculating_deriv = 1; if (use.Get_gas_phase_ptr() != NULL) { @@ -1969,15 +1976,14 @@ jacobian_pz(void) base_phases[i] = *phase_ptr; } } -Restart: size_t pz_max_unknowns = max_unknowns; //k_temp(tc_x, patm_x); - if (full_pitzer == TRUE) - { - molalities(TRUE); - pitzer(); - residuals(); - } + //if (full_pitzer == TRUE) + //{ + // molalities(TRUE); + // pitzer(); + // residuals(); + //} base.resize(count_unknowns); for (i = 0; i < count_unknowns; i++) { @@ -2043,9 +2049,11 @@ Restart: case GAS_MOLES: if (gas_in == FALSE) continue; - d2 = d * x[i]->moles; - if (d2 < 1e-14) - d2 = 1e-14; + d2 = (x[i]->moles > 1 ? 1 : 30); + d2 *= d * x[i]->moles; + d2 = (d2 < ineq_tol ? ineq_tol : d2); + //if (d2 < 1e-14) + // d2 = 1e-14; x[i]->moles += d2; break; case MU: @@ -2149,12 +2157,22 @@ Restart: *phase_ptrs[g] = base_phases[g]; } } + molalities(TRUE); + if (full_pitzer == TRUE) + pitzer(); + mb_sums(); + residuals(); } - molalities(TRUE); - if (full_pitzer == TRUE) - pitzer(); - mb_sums(); - residuals(); + //molalities(TRUE); + //if (full_pitzer == TRUE) + // pitzer(); + //mb_sums(); + //residuals(); + //for (i = 0; i < count_unknowns; i++) + //{ + // residual[i] = base[i]; + // my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i]; + //} base.clear(); calculating_deriv = 0; return OK; From 0dde2b010f074212772cc35f5aa6fca7ae706634 Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Thu, 15 Jul 2021 11:44:28 -0600 Subject: [PATCH 07/11] removed comments --- pitzer.cpp | 17 ----------------- 1 file changed, 17 deletions(-) diff --git a/pitzer.cpp b/pitzer.cpp index ea5c2bd2..905adeb8 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1977,13 +1977,6 @@ Restart: } } size_t pz_max_unknowns = max_unknowns; - //k_temp(tc_x, patm_x); - //if (full_pitzer == TRUE) - //{ - // molalities(TRUE); - // pitzer(); - // residuals(); - //} base.resize(count_unknowns); for (i = 0; i < count_unknowns; i++) { @@ -2163,16 +2156,6 @@ Restart: mb_sums(); residuals(); } - //molalities(TRUE); - //if (full_pitzer == TRUE) - // pitzer(); - //mb_sums(); - //residuals(); - //for (i = 0; i < count_unknowns; i++) - //{ - // residual[i] = base[i]; - // my_array[((size_t)i + 1) * (count_unknowns + 1) - 1] = residual[i]; - //} base.clear(); calculating_deriv = 0; return OK; From 6bd936e0d7ba3497eed22420c7f7a00220d615da Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Thu, 15 Jul 2021 20:37:58 -0600 Subject: [PATCH 08/11] Fixed numerical derivative (non-pitzer) --- model.cpp | 243 ++++++++++++++++++++++++++++++++++++++++++++++++++++- pitzer.cpp | 3 +- 2 files changed, 243 insertions(+), 3 deletions(-) diff --git a/model.cpp b/model.cpp index adeea063..2c691a75 100644 --- a/model.cpp +++ b/model.cpp @@ -5304,7 +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); } - +#ifdef ORIGINAL /* ---------------------------------------------------------------------- */ int Phreeqc:: numerical_jacobian(void) @@ -5545,7 +5545,248 @@ numerical_jacobian(void) calculating_deriv = FALSE; return OK; } +#endif +/* ---------------------------------------------------------------------- */ +int Phreeqc:: +numerical_jacobian(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; + if (! + (numerical_deriv || + (use.Get_surface_ptr() != NULL && use.Get_surface_ptr()->Get_type() == cxxSurface::CD_MUSIC) || + (gas_phase_ptr != NULL && gas_phase_ptr->Get_type() == cxxGasPhase::GP_VOLUME && + (gas_phase_ptr->Get_pr_in() || force_numerical_fixed_volume) && numerical_fixed_volume) + )) + return(OK); + + calculating_deriv = TRUE; + jacobian_sums(); + 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; + } + } + //gammas(mu_x); + molalities(TRUE); + mb_sums(); + //mb_gases(); + //mb_ss(); + residuals(); + /* + * Clear array, note residuals are in array[i, count_unknowns+1] + */ + + //for (i = 0; i < count_unknowns; i++) + //{ + // my_array[i] = 0.0; + //} + //for (i = 1; i < count_unknowns; i++) + //{ + // memcpy((void*)&(my_array[(size_t)i * (count_unknowns + 1)]), + // (void*)&(my_array[0]), count_unknowns * sizeof(LDBLE)); + //} + + base.resize(count_unknowns); + base = residual; + 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 = d;// *LOG_10; + break; + case MH: + s_eminus->la += d; + d2 = d;// *LOG_10; + break; + case AH2O: + x[i]->master[0]->s->la += d; + d2 = d;// *LOG_10; + break; + case PITZER_GAMMA: + x[i]->s->lg += d; + d2 = d; + break; + case MH2O: + //mass_water_aq_x *= (1 + d); + //x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; + //d2 = log(1.0 + d); + //break; + // DL_pitz + d1 = mass_water_aq_x * d; + mass_water_aq_x += d1; + if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL) + mass_water_bulk_x += d1; + x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; + d2 = d1; + break; + case MU: + d2 = d * mu_x; + mu_x += d2; + gammas(mu_x); + break; + case PP: + for (j = 0; j < count_unknowns; j++) + { + delta[j] = 0.0; + } + d2 = -1e-8; + delta[i] = d2; + reset(); + d2 = delta[i]; + break; + case SS_MOLES: + if (x[i]->ss_in == FALSE) + continue; + for (j = 0; j < count_unknowns; j++) + { + delta[j] = 0.0; + } + /*d2 = -1e-8; */ + d2 = d * 10 * x[i]->moles; + //d2 = -.1 * x[i]->moles; + /* + if (d2 > -1e-10) d2 = -1e-10; + calculating_deriv = FALSE; + */ + delta[i] = d2; + /*fprintf (stderr, "delta before reset %e\n", delta[i]); */ + reset(); + d2 = delta[i]; + /*fprintf (stderr, "delta after reset %e\n", delta[i]); */ + break; + case GAS_MOLES: + if (gas_in == FALSE) + continue; + d2 = (x[i]->moles > 1 ? 1 : 30); + d2 *= d * x[i]->moles; + d2 = (d2 < ineq_tol ? ineq_tol : d2); + //if (d2 < 1e-14) + // d2 = 1e-14; + x[i]->moles += d2; + break; + } + //gammas(mu_x); + molalities(TRUE); + mb_sums(); + //mb_gases(); + //mb_ss(); + residuals(); + + for (j = 0; j < count_unknowns; j++) + { + my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] = -(residual[j] - base[j]) / d2; + if (x[i]->type == MH2O) // DL_pitz + my_array[(size_t)j * (count_unknowns + 1) + (size_t)i] *= mass_water_aq_x; + } + 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 -= d2; + break; + case MH: + s_eminus->la -= d2; + if (my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] == 0) + { + /*output_msg(sformatf( "Zero diagonal for MH\n")); */ + my_array[(size_t)i * (count_unknowns + 1) + (size_t)i] = + under(s_h2->lm) * 2; + } + break; + case PITZER_GAMMA: + x[i]->s->lg -= d2; + break; + case MH2O: + //mass_water_aq_x /= (1 + d); + //x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; + //break; + //DL_pitz + mass_water_aq_x -= d2; + if (use.Get_surface_in() && dl_type_x == cxxSurface::DONNAN_DL) + mass_water_bulk_x -= d2; + x[i]->master[0]->s->moles = mass_water_aq_x / gfw_water; + break; + case MU: + mu_x -= d2; + gammas(mu_x); + break; + case PP: + delta[i] = -d2; + reset(); + break; + case SS_MOLES: + delta[i] = -d2; + reset(); + break; + case GAS_MOLES: + 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]; + } + } + 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; + // } + //} + } + base.clear(); + calculating_deriv = FALSE; + return OK; +} /* ---------------------------------------------------------------------- */ void Phreeqc:: set_inert_moles(void) diff --git a/pitzer.cpp b/pitzer.cpp index 905adeb8..81f5f098 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1947,10 +1947,9 @@ int Phreeqc:: jacobian_pz(void) /* ---------------------------------------------------------------------- */ { // calculate the derivatives numerically - std::vector base, base_x, base_gas_moles; + std::vector base; std::vector phase_ptrs; std::vector base_phases; - double base_mass_water_bulk_x = 0, base_moles_h2o = 0; cxxGasPhase base_gas_phase; LDBLE d, d1, d2; int i, j; From df0d68b9d5cfb3df06c0a44118d0b8c39ecda91f Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Fri, 16 Jul 2021 22:27:03 -0600 Subject: [PATCH 09/11] 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) From 56975a7e0c6ee92cec5342a9e036fab72914f75f Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Sat, 17 Jul 2021 08:49:17 -0600 Subject: [PATCH 10/11] Saved surface for numerical derivatives --- model.cpp | 9 +++++++++ pitzer.cpp | 27 ++++++++++++++++++--------- sit.cpp | 10 +++++++++- 3 files changed, 36 insertions(+), 10 deletions(-) diff --git a/model.cpp b/model.cpp index 88d33609..39955f26 100644 --- a/model.cpp +++ b/model.cpp @@ -5560,6 +5560,7 @@ numerical_jacobian(void) std::vector base_phases; double base_mass_water_bulk_x = 0, base_moles_h2o = 0; cxxGasPhase base_gas_phase; + cxxSurface base_surface; if (! (numerical_deriv || @@ -5571,6 +5572,10 @@ numerical_jacobian(void) calculating_deriv = TRUE; //jacobian_sums(); + if (use.Get_surface_ptr() != NULL) + { + base_surface = *use.Get_surface_ptr(); + } if (use.Get_gas_phase_ptr() != NULL) { //cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); @@ -5761,6 +5766,10 @@ numerical_jacobian(void) x[i]->moles -= d2; break; } + if (use.Get_surface_ptr() != NULL) + { + base_surface = *use.Get_surface_ptr(); + } if (use.Get_gas_phase_ptr() != NULL) { *use.Get_gas_phase_ptr() = base_gas_phase; diff --git a/pitzer.cpp b/pitzer.cpp index 21d667aa..b84aa850 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1951,6 +1951,7 @@ jacobian_pz(void) std::vector phase_ptrs; std::vector base_phases; cxxGasPhase base_gas_phase; + cxxSurface base_surface; LDBLE d, d1, d2; int i, j; Restart: @@ -1962,6 +1963,10 @@ Restart: } mb_sums(); residuals(); + if (use.Get_surface_ptr() != NULL) + { + base_surface = *use.Get_surface_ptr(); + } if (use.Get_gas_phase_ptr() != NULL) { cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); @@ -2141,6 +2146,10 @@ Restart: reset(); break; } + if (use.Get_surface_ptr() != NULL) + { + base_surface = *use.Get_surface_ptr(); + } if (use.Get_gas_phase_ptr() != NULL) { *use.Get_gas_phase_ptr() = base_gas_phase; @@ -2155,15 +2164,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; - // } - //} + 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-6) + { + 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 a1b140b6..ec6d328f 100644 --- a/sit.cpp +++ b/sit.cpp @@ -967,7 +967,12 @@ jacobian_sit(void) std::vector base_phases; double base_mass_water_bulk_x = 0, base_moles_h2o = 0; cxxGasPhase base_gas_phase; + cxxSurface base_surface; Restart: + if (use.Get_surface_ptr() != NULL) + { + base_surface = *use.Get_surface_ptr(); + } if (use.Get_gas_phase_ptr() != NULL) { cxxGasPhase* gas_phase_ptr = use.Get_gas_phase_ptr(); @@ -1120,7 +1125,10 @@ Restart: x[i]->moles -= d2; break; } - + if (use.Get_surface_ptr() != NULL) + { + base_surface = *use.Get_surface_ptr(); + } if (use.Get_gas_phase_ptr() != NULL) { *use.Get_gas_phase_ptr() = base_gas_phase; From 07a864d1a6ed582ddfae5e1c3bf5b0620c3dd75c Mon Sep 17 00:00:00 2001 From: David Parkhurst Date: Sat, 17 Jul 2021 12:21:02 -0600 Subject: [PATCH 11/11] all jacobians are consistent. Looks pretty good. --- model.cpp | 6 +++--- pitzer.cpp | 52 +++++++++++++++++++++++++++++----------------------- sit.cpp | 4 ++-- 3 files changed, 34 insertions(+), 28 deletions(-) diff --git a/model.cpp b/model.cpp index 39955f26..1d31226a 100644 --- a/model.cpp +++ b/model.cpp @@ -5570,7 +5570,6 @@ numerical_jacobian(void) )) return(OK); - calculating_deriv = TRUE; //jacobian_sums(); if (use.Get_surface_ptr() != NULL) { @@ -5589,6 +5588,7 @@ numerical_jacobian(void) base_phases[i] = *phase_ptr; } } + calculating_deriv = TRUE; gammas(mu_x); molalities(TRUE); mb_sums(); @@ -5768,7 +5768,7 @@ numerical_jacobian(void) } if (use.Get_surface_ptr() != NULL) { - base_surface = *use.Get_surface_ptr(); + *use.Get_surface_ptr() = base_surface; } if (use.Get_gas_phase_ptr() != NULL) { @@ -5797,7 +5797,7 @@ numerical_jacobian(void) // 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; + // std::cerr << iterations << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl; // } //} base.clear(); diff --git a/pitzer.cpp b/pitzer.cpp index b84aa850..68a57a39 100644 --- a/pitzer.cpp +++ b/pitzer.cpp @@ -1955,14 +1955,6 @@ jacobian_pz(void) LDBLE d, d1, d2; int i, j; Restart: - calculating_deriv = 1; - molalities(TRUE); - if (full_pitzer == TRUE) - { - pitzer(); - } - mb_sums(); - residuals(); if (use.Get_surface_ptr() != NULL) { base_surface = *use.Get_surface_ptr(); @@ -1980,6 +1972,15 @@ Restart: base_phases[i] = *phase_ptr; } } + calculating_deriv = 1; + molalities(TRUE); + if (full_pitzer == TRUE) + { + pitzer(); + } + mb_sums(); + residuals(); + size_t pz_max_unknowns = max_unknowns; base.resize(count_unknowns); for (i = 0; i < count_unknowns; i++) @@ -2148,7 +2149,7 @@ Restart: } if (use.Get_surface_ptr() != NULL) { - base_surface = *use.Get_surface_ptr(); + *use.Get_surface_ptr() = base_surface; } if (use.Get_gas_phase_ptr() != NULL) { @@ -2158,21 +2159,26 @@ Restart: *phase_ptrs[g] = base_phases[g]; } } - molalities(TRUE); - if (full_pitzer == TRUE) - pitzer(); - 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-6) - { - std::cerr << i << ": " << x[i]->description << " " << residual[i] << " " << base[i] << std::endl; - } + //molalities(TRUE); + //if (full_pitzer == TRUE) + // pitzer(); + //mb_sums(); + //residuals(); } + molalities(TRUE); + if (full_pitzer == TRUE) + pitzer(); + 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-6) + // { + // std::cerr << iterations << ": " << 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 ec6d328f..fb372c3b 100644 --- a/sit.cpp +++ b/sit.cpp @@ -986,9 +986,9 @@ Restart: base_phases[i] = *phase_ptr; } } + calculating_deriv = 1; size_t pz_max_unknowns = max_unknowns; //k_temp(tc_x, patm_x); - calculating_deriv = 1; molalities(TRUE); if (full_pitzer == TRUE) { @@ -1127,7 +1127,7 @@ Restart: } if (use.Get_surface_ptr() != NULL) { - base_surface = *use.Get_surface_ptr(); + *use.Get_surface_ptr() = base_surface; } if (use.Get_gas_phase_ptr() != NULL) {