/Software/PhreeqC/modelengine.cpp
C++ | 12907 lines | 9019 code | 2338 blank | 1550 comment | 3087 complexity | 3812a53d6a5af9c4679ab2ee6b2180c9 MD5 | raw file
Large files files are truncated, but you can click here to view the full file
- #include <math.h>
- #include <float.h>
- #include "modelengine.h"
- #include "conio.h"
-
- //===========================
- // Constructors & Destructors
- //===========================
- ModelEngine::ModelEngine(GlobalData *gd, ModelData *md):
- SSAssemblageModel(gd, md)
- {
- this->gd = gd;
- this->md = md;
-
- SetData(md);
-
- r_temp = NULL;
- r_temp = new Reaction;
-
- count_unknowns = 0;
-
- error = NULL;
- error_string[0] = '\0';
- }
-
- ModelEngine::~ModelEngine()
- {
- delete r_temp;
- }
-
- //===========================
- // Public functions
- //===========================
- bool ModelEngine::ConvertUnits(Solution *sol_p)
- {
- int i;
- LDBLE sum_solutes;
- Master *m;
- Conc *conc;
-
- sum_solutes = exp(-sol_p->ph * md->LOG_10);
-
- for (i = 0; i < sol_p->totals->Count(); i++)
- {
- conc = (*sol_p->totals)[i];
-
- conc->moles = 0.0;
- if (conc->name == "H(1)" || conc->name == "E")
- continue;
-
- if (conc->input_conc <= 0)
- continue;
-
- if (conc->gfw <= 0.0)
- {
- if (!conc->as.IsEmpty())
- {
- if (!ComputeGFW(conc->as, conc->gfw))
- return false;
-
- if (conc->name == "Alkalinity" && conc->as == "CaCO3")
- conc->gfw /= 2.0;
- }
- else
- {
- if ((m = gd->master_list.Search(&conc->name(0, " "), true)) != NULL)
- conc->gfw = m->gfw;
- else
- return false;
- }
- }
-
- conc->moles = conc->input_conc;
-
- if (Check_l_(sol_p->units))
- conc->moles *= (LDBLE)1.0 / (sol_p->density);
-
- conc->moles *= ConversionFactorForMoles(conc->units);
-
- if (Check_g_kgs_(conc->units) || Check_g_l_(conc->units))
- sum_solutes += conc->moles;
- else if (Check_mol_kgs_(conc->units) || Check_mol_l_(conc->units) || Check_eq_l_(conc->units))
- sum_solutes += (conc->moles * conc->gfw);
-
- if (Check_g_(conc->units) && conc->gfw != 0.0)
- conc->moles /= conc->gfw;
- }
-
- if (Check_kgs_(sol_p->units) || Check_l_(sol_p->units))
- {
- md->mass_water_aq_x = (LDBLE)1.0 - ((LDBLE)1e-3 * sum_solutes);
- for (i = 0; i < sol_p->totals->Count(); i++)
- (*sol_p->totals)[i]->moles /= md->mass_water_aq_x;
- }
-
- md->mass_water_aq_x = sol_p->mass_water;
- for (i = 0; i < sol_p->totals->Count(); i++)
- (*sol_p->totals)[i]->moles *= md->mass_water_aq_x;
-
- return true;
- }
-
-
- bool ModelEngine::Clear()
- {
- int i;
- Solution *sol_p = md->use.sol_p;
-
- // Clear species solution-dependent data
- for (i = 0; i < gd->species_list.Count(); i++)
- {
- gd->species_list[i]->in = false;
- //gd->species_list[i]->rewrite = false;
- }
-
- // Set pe structure
- sol_p->pe->CopyTo(md->pe_x);
- md->default_pe_x = sol_p->default_pe;
-
- // Clear master species solution-dependent data
- Master *m;
- PEData *ped;
- for (i = 0; i < gd->master_list.Count(); i++)
- {
- m = gd->master_list[i];
-
- m->in = false;
- m->rewrite = false;
- m->u = NULL;
-
- ped = (*md->pe_x)[sol_p->default_pe];
- m->pe_rxn = &(ped->rxn);
-
- // copy primary reaction to secondary reaction
- m->rxn_primary->CopyTo(m->rxn_secondary);
- }
-
- if (md->state == INITIAL_SOLUTION)
- {
- gd->s_h2o->secondary->in = true;
- gd->s_h2o->secondary->rewrite = false;
- gd->s_hplus->secondary->in = true;
- gd->s_hplus->secondary->rewrite = false;
- }
- else
- {
- gd->s_h2o->primary->in = true;
- gd->s_h2o->primary->rewrite = false;
- gd->s_hplus->primary->in = true;
- gd->s_hplus->primary->rewrite = false;
- }
- gd->s_eminus->primary->in = true;
- gd->s_eminus->primary->rewrite = false;
-
- // Set all unknown pointers to NULL
- md->mb_unknown = NULL;
- md->ah2o_unknown = NULL;
- md->mass_hydrogen_unknown = NULL;
- md->mass_oxygen_unknown = NULL;
- md->mu_unknown = NULL;
- md->alkalinity_unknown = NULL;
- md->carbon_unknown = NULL;
- md->ph_unknown = NULL;
- md->pe_unknown = NULL;
- md->charge_balance_unknown = NULL;
- md->solution_phase_boundary_unknown = NULL;
- md->pure_phase_unknown = NULL;
- md->exchange_unknown = NULL;
- md->surface_unknown = NULL;
- md->gas_unknown = NULL;
- md->s_s_unknown = NULL;
-
- // Free arrays used in model
- FreeModelAllocs(); //free_model_allocs ();
-
- return true;
- }
-
-
- //===========================
- // Private functions
- //===========================
- bool ModelEngine::GetMasterList(String &list, Master *m, ListOfPointers<Master> *copy)
- {
- String *token;
- Master *m_ptr, *m_1;
- int i;
-
- tl.SetString(list);
- copy->Clear();
-
- // Make list of master species pointers
- if (m == m->s->primary)
- {
- // First in list is primary species
- for (i = 0; i < gd->master_list.Count(); i++)
- {
- m_1 = gd->master_list[i];
- if (m_1 == m)
- break;
- }
- i++;
-
- if (i >= gd->master_list.Count())
- {
- // Element has only one valence
- copy->AddNew(m);
- }
- else
- {
- m_1 = gd->master_list[i];
-
- if (m != m_1->e->primary)
- {
- // Element has only one valence
- copy->AddNew(m);
- }
- else
- {
- // Element has multiple valences
- if (m->s->secondary == NULL)
- return false;
-
- copy->AddNew(m->s->secondary);
-
- while (i < gd->master_list.Count() && (m_ptr = gd->master_list[i])->e->primary == m)
- {
- if (m_ptr->s->primary == NULL)
- copy->AddNew(m_ptr);
-
- i++;
- }
- }
- }
- }
- else
- {
- // First in list is secondary species, Include all valences from input
- copy->AddNew(m);
-
- while (!tl.EOTL())
- {
- token = tl.NextToken();
- m_ptr = gd->master_list.Search(token, true);
-
- if (m_ptr != NULL)
- copy->AddNew(m_ptr);
- }
- }
-
- return true;
- }
-
- bool ModelEngine::RewriteMasterToSecondary(Master *m1, Master *m2)
- {
- //
- // Write equation for secondary master species in terms of another secondary master species
- // Store result in rxn_secondary of master_ptr.
- //
-
- LDBLE coef1, coef2;
- Master *m_p1, *m_p2;
-
- // Check that the two master species have the same primary master species
- m_p1 = m1->e->primary;
- m_p2 = m2->e->primary;
-
- if (m_p1 != m_p2 || m_p1 == NULL)
- return false;
-
- // Find coefficient of primary master in reaction
- if (!GetReactionCoef(m1->rxn_primary, m_p1->s->name, coef1) || !GetReactionCoef(m2->rxn_primary, m_p1->s->name, coef2))
- return false;
-
- if (Equal(coef1, (LDBLE)0.0, (LDBLE)TOLERANCE) || Equal(coef2, (LDBLE)0.0, (LDBLE)TOLERANCE))
- return false;
-
- // Rewrite equation to secondary master species
- r_temp->AddTRXN(m1->rxn_primary, 1.0, false);
- r_temp->AddTRXN(m2->rxn_primary, -coef1 / coef2, true);
-
- return false;
- }
-
- bool ModelEngine::GetReactionCoef(Reaction *r, String name, LDBLE &coef, int start)
- {
- //
- // Finds coefficient of token in reaction.
- // input: r, pointer to a reaction structure
- // name, string to find as reaction token
- // start, search start position in r (default = 1)
- //
- // Return: 0.0, if token not found
- // coefficient of token, if found.
- //
-
- coef = 0.0;
-
- ReactionToken *r_token;
- for (int i = start; i < r->token_list->Count(); i++)
- {
- r_token = r->token_list->Element(i);
-
- if (r_token->s->name == name)
- {
- coef = r_token->coef;
- break;
- }
- }
-
- return true;
- }
- bool ModelEngine::CheckIn()
- {
- //
- // Routine goes through trxn to determine if each master species is in this model.
- // Assumes equation is written in terms of primary and secondary species
- // Checks to see if in is TRUE or REWRITE for each species
- // Returns TRUE if in model FALSE if not
- //
-
- ReactionToken *t;
-
- for (int i = 1; i < r_temp->token_list->Count(); i++)
- {
- t = r_temp->token_list->Element(i);
-
- // Check primary master species in
- if (t->s->primary != NULL && t->s->primary->in)
- continue;
-
- // Check secondary master species
- if (t->s->secondary != NULL && (t->s->secondary->in || t->s->secondary->rewrite))
- continue;
-
- // Must be primary master species that is out
- return false;
- }
-
- return true;
- }
-
- bool ModelEngine::WriteMassActionEqnX()
- {
- //
- // Reduce mass-action equation to the master species that are in the model
- //
-
- LDBLE coef_e;
- int count;
- int i, count_rxn_orig;
-
-
- // Rewrite any secondary master species flagged REWRITE. Replace pe if necessary
- count = 0;
- bool repeat = true;
- while (repeat)
- {
- count++;
- if (count > MAX_ADD_EQUATIONS)
- return false;
-
- repeat = false;
-
- ReactionToken *r;
- count_rxn_orig = r_temp->token_list->Count();
- for (i = 1; i < count_rxn_orig; i++)
- {
- r = r_temp->token_list->Element(i);
- if (r->s->secondary == NULL)
- continue;
-
- if (r->s->secondary->rewrite)
- {
- repeat = true;
-
- if (!GetReactionCoef(r->s->secondary->rxn_secondary, "e-", coef_e))
- return false;
-
- r_temp->AddTRXN(r->s->secondary->rxn_secondary, r->coef, false);
-
- if (!Equal(coef_e, (LDBLE)0.0, (LDBLE)TOLERANCE))
- r_temp->AddTRXN(*r->s->secondary->pe_rxn, r->coef * coef_e, false);
- }
- }
-
- r_temp->TRXNCombine();
- }
-
- return true;
- }
-
- bool ModelEngine::AddPotentialFactor()
- {
- //
- // Add the potential factor to surface mass-action equations.
- // Factor is essentially the activity coefficient, representing
- // the work required to bring charged ions to the surface
- //
-
- int i;
- LDBLE sum_z;
- Master *master_ptr;
- Unknown *unknown_ptr;
-
- if (md->use.sur_p->type != DDL)
- return true;
-
- sum_z = 0.0;
- master_ptr = NULL;
-
- // Find sum of charge of aqueous species and surface master species
- ReactionToken *r;
- for (i = 1; i < r_temp->token_list->Count(); i++)
- {
- r = r_temp->token_list->Element(i);
-
- if (r->s->type == AQ || r->s == gd->s_hplus || r->s == gd->s_eminus)
- sum_z += r->s->z * r->coef;
-
- if (r->s->type == SURF)
- master_ptr = r->s->primary;
- }
-
- if (master_ptr == NULL)
- return false;
-
- // Find potential unknown for surface species
- unknown_ptr = FindSurfaceChargeUnknown(master_ptr->name, SURF_PSI);
-
- if (unknown_ptr == NULL)
- return false;
-
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- // Include psi in mass action equation
- if (master_ptr != NULL)
- {
- r = r_temp->token_list->AddNew();
-
- r->name = master_ptr->s->name;
- r->s = master_ptr->s;
- r->coef = (LDBLE)-2.0 * sum_z;
- }
- else
- return false;
-
- return true;
- }
-
- bool ModelEngine::AddCDMusicFactors(int n)
- {
- //
- // Add the potential factors for cd_music to surface mass-action equations.
- // Factors are essentially the activity coefficient, representing
- // the work required to bring charged ions to the three charge layers
- // of the cd_music model
- //
-
- int i;
- String name;
- Master *master_ptr;
- Unknown *unknown_ptr;
- ReactionToken *r;
-
- if (md->use.sur_p->type != CD_MUSIC)
- return true;
-
- master_ptr = NULL;
-
- // Find sum of charge of aqueous species and surface master species
- for (i = 1; i < r_temp->token_list->Count(); i++)
- {
- r = r_temp->token_list->Element(i);
-
- if (r->s->type == SURF)
- master_ptr = r->s->primary;
- }
-
- if (master_ptr == NULL)
- return false;
-
- name = master_ptr->e->name;
-
- // Find potential unknown for surface species
- unknown_ptr = FindSurfaceChargeUnknown (name, SURF_PSI);
-
- if (unknown_ptr == NULL)
- return false;
-
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- // Include psi in mass action equation
- r = r_temp->token_list->AddNew();
-
- r->name = master_ptr->s->name;
- r->s = master_ptr->s;
- r->coef = gd->species_list[n]->dz[0];
-
- // Plane 1
- unknown_ptr = FindSurfaceChargeUnknown(name, SURF_PSI1);
-
- if (unknown_ptr == NULL)
- return false;
-
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- // Include psi in mass action equation
- r = r_temp->token_list->AddNew();
-
- r->name = master_ptr->s->name;
- r->s = master_ptr->s;
- r->coef = gd->species_list[n]->dz[1];
-
- // Plane 2
- unknown_ptr = FindSurfaceChargeUnknown(name, SURF_PSI2);
-
- if (unknown_ptr == NULL)
- return false;
-
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- // Include psi in mass action equation
- r = r_temp->token_list->AddNew();
-
- r->name = master_ptr->s->name;
- r->s = master_ptr->s;
- r->coef = gd->species_list[n]->dz[2];
-
- return true;
- }
-
- Unknown *ModelEngine::FindSurfaceChargeUnknown(String &name, int plane)
- {
- int i;
- String token;
-
- token = name;
- token.Replace ("_", " ");
- token = token(0, " ");
-
- if (plane == SURF_PSI)
- token += "_CB";
- else if (plane == SURF_PSI1)
- token += "_CBb";
- else if (plane == SURF_PSI2)
- token += "_CBd";
-
- for (i = 0; i < count_unknowns; i++)
- if (token == (*unknown_list)[i]->name)
- return (*unknown_list)[i];
-
- return NULL;
- }
- bool ModelEngine::WriteMBEqnX()
- {
- //
- // Rewrite any secondary master species flagged REWRITE
- // Don't add in any pe reactions
- //
-
- int count;
- bool repeat;
- int i, count_rxn_orig;
- int j, k;
- char *ptr, token[DEFAULT_STRING_LENGTH];
- Master *master_ptr;
- ReactionToken *rxn_token;
-
- count = 0;
- repeat = true;
- while (repeat)
- {
- count++;
- if (count > MAX_ADD_EQUATIONS)
- return false;
-
- repeat = false;
-
- count_rxn_orig = r_temp->token_list->Count();
- for (i = 1; i < count_rxn_orig; i++)
- {
- rxn_token = r_temp->token_list->Element(i);
-
- if (rxn_token->s->secondary == NULL)
- continue;
-
- if (rxn_token->s->secondary->rewrite)
- {
- repeat = true;
- r_temp->AddTRXN(rxn_token->s->secondary->rxn_secondary, rxn_token->coef, false);
- }
- }
-
- r_temp->TRXNCombine();
- }
-
- eos_list->Clear();
- parent_count = 0;
-
- ReactionToken *r;
- ElementOfSpecies *eos_p;
- for (i = 1; i < r_temp->token_list->Count(); i++)
- {
- r = r_temp->token_list->Element(i);
- j = eos_list->Count();
-
-
- r->s->name.Copy(token);
- ptr = &token[0];
- if (!GetElementsInSpecies(&ptr, r->coef))
- return false;
-
- for (k = j; k < eos_list->Count(); k++)
- {
- eos_p = (*eos_list)[k];
- if (r->s->secondary != NULL)
- master_ptr = r->s->secondary->e->primary;
- else
- master_ptr = r->s->primary;
-
- if (eos_p->e == master_ptr->e)
- {
- eos_p->coef = 0.0;
- break;
- }
- }
-
- if (r->s->secondary == NULL)
- r->s->primary->e->name.Copy(token);
- else
- r->s->secondary->e->name.Copy(token);
-
- ptr = &token[0];
- if (!GetSecondaryElementsInSpecies(&ptr, r->coef))
- return false;
- }
-
- CombineElements();
-
- return true;
- }
-
- bool ModelEngine::AddSurfaceChargeBalance()
- {
- //
- // Include charge balance in list for mass-balance equations
- //
-
- int i;
- char *ptr;
- char token[DEFAULT_STRING_LENGTH];
-
- Master *master_ptr;
- Unknown *unknown_ptr;
- ElementOfSpecies *eos_p;
-
- if (md->use.sur_p->type != DDL)
- return true;
-
- master_ptr = NULL;
-
- // Find master species
- for (i = 0; i < eos_list->Count(); i++)
- {
- eos_p = (*eos_list)[i];
-
- if (eos_p->e->primary->s->type == SURF)
- {
- master_ptr = eos_p->e->primary;
- break;
- }
- }
-
- if (i >= eos_list->Count())
- throw exception();
-
- // Find potential unknown for surface species
- unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI);
-
- if (unknown_ptr == NULL)
- throw exception();
-
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- master_ptr->e->name.Copy(token);
- ptr = &token[0];
-
- // Include charge balance in list for mass-balance equations
- LDBLE coef = 1.0;
- if (!GetSecondaryElementsInSpecies(&ptr, coef))
- return false;
-
- return true;
- }
-
- bool ModelEngine::AddCDMusicChargeBalances(int n)
- {
- //
- // Add the potential factor to surface mass-action equations.
- // Factor is essentially the activity coefficient, representing
- // the work required to bring charged ions to the surface
- //
-
- int i;
- char *ptr;
- char token[DEFAULT_STRING_LENGTH];
-
- Master *master_ptr;
- Unknown *unknown_ptr;
- ElementOfSpecies *eos_p;
-
- if (md->use.sur_p->type != CD_MUSIC)
- return true;
-
- master_ptr = NULL;
-
- // Find master species
- for (i = 0; i < eos_list->Count(); i++)
- {
- eos_p = (*eos_list)[i];
-
- if (eos_p->e->primary->s->type == SURF)
- {
- master_ptr = eos_p->e->primary;
- break;
- }
- }
-
- if (i >= eos_list->Count())
- throw exception();
-
- // Find potential unknown for plane 0
- unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI);
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- // Include charge balance in list for mass-balance equations
- master_ptr->e->name.Copy(token);
- ptr = &token[0];
- if (!GetSecondaryElementsInSpecies(&ptr, gd->species_list[n]->dz[0]))
- return false;
-
- // Find potential unknown for plane 1
- unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI1);
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- // Include charge balance in list for mass-balance equations
- master_ptr->e->name.Copy(token);
- ptr = &token[0];
- if (!GetSecondaryElementsInSpecies(&ptr, gd->species_list[n]->dz[1]))
- return false;
-
- // Find potential unknown for plane 2
- unknown_ptr = FindSurfaceChargeUnknown(master_ptr->e->name, SURF_PSI2);
- master_ptr = (*unknown_ptr->master)[0]; // potential for surface component
-
- // Include charge balance in list for mass-balance equations
- master_ptr->e->name.Copy(token);
- ptr = &token[0];
- if (!GetSecondaryElementsInSpecies(&ptr, gd->species_list[n]->dz[2]))
- return false;
-
- return true;
- }
-
- bool ModelEngine::MBForSpeciesAQ(int n)
- {
- //
- // Make list of mass balance and charge balance equations in which
- // to insert species n.
- //
- // count_mb_unknowns - number of equations and summation relations
- // mb_unknowns.unknown - pointer to unknown which contains row number
- // mb_unknowns.source - pointer to the LDBLE number to be multiplied
- // by coef, usually moles.
- // mb_unknowns.coef - coefficient of s[n] in equation or relation
-
- int i, j;
- Species *s;
- Master *master_ptr;
- Unknown *unknown_ptr;
-
- md->mb_unknowns->Clear();
-
- s = gd->species_list[n];
-
- // e- does not appear in any mass balances
- if (s->type == EMINUS)
- return true;
-
- // Do not include diffuse layer in cb, alk, ah2o, mu
- if (md->charge_balance_unknown != NULL && s->type < H2O)
- StoreMBUnknowns(md->charge_balance_unknown, &s->moles, s->z, &s->dg);
-
- if (md->alkalinity_unknown != NULL && s->type < H2O)
- StoreMBUnknowns (md->alkalinity_unknown, &s->moles, s->alk, &s->dg);
-
- if (md->ah2o_unknown != NULL && s->type < H2O)
- StoreMBUnknowns (md->ah2o_unknown, &s->moles, 1.0, &s->dg);
-
- if (md->mu_unknown != NULL && s->type < H2O)
- StoreMBUnknowns (md->mu_unknown, &s->moles, s->z * s->z, &s->dg);
-
- // Include diffuse layer in hydrogen and oxygen mass balance
- if (md->mass_hydrogen_unknown != NULL)
- {
- if (md->dl_type_x != NO_DL && md->state >= REACTION)
- StoreMBUnknowns (md->mass_hydrogen_unknown, &s->tot_g_moles, s->h - 2 * s->o, &s->dg_total_g);
- else
- StoreMBUnknowns (md->mass_hydrogen_unknown, &s->moles, s->h - 2 * s->o, &s->dg);
- }
-
- if (md->mass_oxygen_unknown != NULL)
- {
- if (md->dl_type_x != NO_DL && md->state >= REACTION)
- StoreMBUnknowns (md->mass_oxygen_unknown, &s->tot_g_moles, s->o, &s->dg_total_g);
- else
- StoreMBUnknowns (md->mass_oxygen_unknown, &s->moles, s->o, &s->dg);
- }
-
- // Sum diffuse layer charge into (surface + DL) charge balance
- SpeciesDiffLayer *sdl_p;
- if (md->use.sur_p != NULL && s->type < H2O && md->dl_type_x != NO_DL)
- {
- j = 0;
- for (i = 0; i < count_unknowns; i++)
- {
- unknown_ptr = (*unknown_list)[i];
-
- if (unknown_ptr->type == SURFACE_CB)
- {
- if (md->use.sur_p->type == CD_MUSIC)
- unknown_ptr = (*unknown_list)[i + 2];
-
- sdl_p = (*s->diff_layer)[j];
- StoreMBUnknowns(unknown_ptr, &sdl_p->g_moles, s->z, &sdl_p->dg_g_moles);
- j++;
- }
- }
- }
-
- // Other mass balances
- ElementOfSpecies *eos;
- for (i = 0; i < eos_list->Count(); i++)
- {
- eos = (*eos_list)[i];
-
- if (eos->e->master->s->type > AQ && eos->e->master->s->type < SOLID)
- continue;
-
- master_ptr = eos->e->master;
-
- if (master_ptr->primary && master_ptr->s->secondary != NULL)
- master_ptr = master_ptr->s->secondary;
-
- if (master_ptr->u == md->ph_unknown)
- continue;
- else if (master_ptr->u == md->pe_unknown)
- continue;
- else if (master_ptr->u == md->charge_balance_unknown)
- continue;
- else if (master_ptr->u == md->alkalinity_unknown)
- continue;
- else if (master_ptr->u->type == SOLUTION_PHASE_BOUNDARY)
- continue;
-
- if (md->dl_type_x != NO_DL && md->state >= REACTION)
- StoreMBUnknowns(master_ptr->u, &s->tot_g_moles, eos->coef * master_ptr->coef, &s->dg_total_g);
- else
- StoreMBUnknowns(master_ptr->u, &s->moles, eos->coef * master_ptr->coef, &s->dg);
- }
-
- return true;
- }
-
- bool ModelEngine::StoreMBUnknowns(Unknown *unknown_ptr, LDBLE *LDBLE_ptr, LDBLE coef, LDBLE *gamma_ptr)
- {
- //
- // Takes an unknown pointer and a coefficient and puts in
- // list of mb_unknowns
- //
-
- if (Equal(coef, (LDBLE)0.0, (LDBLE)TOLERANCE))
- return true;
-
- UnknownInfo *m = md->mb_unknowns->AddNew();
-
- m->u = unknown_ptr;
- m->source = LDBLE_ptr;
- m->gamma_source = gamma_ptr;
- m->coef = coef;
-
- return true;
- }
-
- bool ModelEngine::MBForSpeciesEX(int n)
- {
- //
- // Make list of mass balance and charge balance equations in which
- // to insert exchange species n.
- //
- // count_mb_unknowns - number of equations and summation relations
- // mb_unknowns.source - pointer to the LDBLE number to be multiplied
- // by coef, usually moles.
- // mb_unknowns.unknown - pointer to unknown which contains row number
- // mb_unknowns.coef - coefficient of s[n] in equation or relation
- //
-
- int i;
- Master *master_ptr;
- Species *s;
-
- md->mb_unknowns->Clear();
-
- // Master species for exchange do not appear in any mass balances
- s = gd->species_list[n];
- if (s->type == EX && s->primary != NULL)
- return true;
-
- // Include diffuse layer in hydrogen and oxygen mass balance
- if (md->charge_balance_unknown != NULL)
- StoreMBUnknowns(md->charge_balance_unknown, &s->moles, s->z, &s->dg);
-
- if (md->mass_hydrogen_unknown != NULL)
- StoreMBUnknowns(md->mass_hydrogen_unknown, &s->moles, s->h - 2 * s->o, &s->dg);
-
- if (md->mass_oxygen_unknown != NULL)
- StoreMBUnknowns(md->mass_oxygen_unknown, &s->moles, s->o, &s->dg);
-
- // Other mass balances
- ElementOfSpecies *eos_p;
- for (i = 0; i < eos_list->Count(); i++)
- {
- eos_p = (*eos_list)[i];
- if (eos_p->e->master->s->type > AQ && eos_p->e->master->s->type < SOLID)
- continue;
-
- master_ptr = eos_p->e->master;
-
- if (master_ptr->primary && master_ptr->s->secondary != NULL)
- master_ptr = master_ptr->s->secondary;
-
- // Special for ph_unknown, pe_unknown, and alkalinity_unknown
- if (master_ptr->u == md->ph_unknown)
- continue;
- else if (master_ptr->u == md->pe_unknown)
- continue;
- else if (master_ptr->u == md->alkalinity_unknown)
- continue;
-
- // EX, sum exchange species only into EXCH mass balance in initial calculation
- // into all mass balances in reaction calculation
- if (md->state >= REACTION || master_ptr->s->type == EX)
- StoreMBUnknowns(master_ptr->u, &s->moles, eos_p->coef * master_ptr->coef, &s->dg);
- }
-
- return true;
- }
-
- bool ModelEngine::MBForSpeciesSURF(int n)
- {
- //
- // Make list of mass balance and charge balance equations in which
- // to insert species n.
- //
- // count_mb_unknowns - number of equations and summation relations
- // mb_unknowns.source - pointer to the LDBLE number to be multiplied
- // by coef, usually moles.
- // mb_unknowns.unknown - pointer to unknown which contains row number
- // mb_unknowns.coef - coefficient of s[n] in equation or relation
- //
-
- int i;
- Master *master_ptr;
-
- md->mb_unknowns->Clear();
-
- // Include in charge balance, if diffuse_layer_x == false
- Species *s = gd->species_list[n];
- if (md->charge_balance_unknown != NULL && md->dl_type_x == NO_DL)
- StoreMBUnknowns(md->charge_balance_unknown, &s->moles, s->z, &s->dg);
-
- // Include diffuse layer in hydrogen and oxygen mass balance
- if (md->mass_hydrogen_unknown != NULL)
- StoreMBUnknowns (md->mass_hydrogen_unknown, &s->moles, s->h - 2 * s->o, &s->dg);
-
- if (md->mass_oxygen_unknown != NULL)
- StoreMBUnknowns (md->mass_oxygen_unknown, &s->moles, s->o, &s->dg);
-
- // Other mass balances
- ElementOfSpecies *eos_p;
- for (i = 0; i < eos_list->Count(); i++)
- {
- eos_p = (*eos_list)[i];
-
- // Skip H+, e-, and H2O
- if (eos_p->e->master->s->type > AQ && eos_p->e->master->s->type < SOLID)
- continue;
-
- master_ptr = eos_p->e->master;
- if (master_ptr->primary && master_ptr->s->secondary != NULL)
- master_ptr = master_ptr->s->secondary;
-
- // SURF_PSI, sum surface species in (surface + DL) charge balance
- if (master_ptr->s->type == SURF_PSI && md->use.sur_p->type != CD_MUSIC)
- {
- StoreMBUnknowns (master_ptr->u, &s->moles, s->z, &s->dg);
- continue;
- }
-
- if (master_ptr->s->type == SURF_PSI && md->use.sur_p->type == CD_MUSIC)
- {
- StoreMBUnknowns (master_ptr->u, &s->moles, s->dz[0], &s->dg);
- continue;
- }
-
- if (master_ptr->s->type == SURF_PSI1)
- {
- StoreMBUnknowns (master_ptr->u, &s->moles, s->dz[1], &s->dg);
- continue;
- }
-
- if (master_ptr->s->type == SURF_PSI2)
- {
- StoreMBUnknowns (master_ptr->u, &s->moles, s->dz[2], &s->dg);
- continue;
- }
-
- // Special for ph_unknown, pe_unknown, and alkalinity_unknown
- if (master_ptr->u == md->ph_unknown)
- continue;
- else if (master_ptr->u == md->pe_unknown)
- continue;
- else if (master_ptr->u == md->alkalinity_unknown)
- continue;
-
- // SURF, sum surface species only into SURFACE mass balance in initial calculation
- // into all mass balances in reaction calculation
- if (md->state >= REACTION || master_ptr->s->type == SURF)
- {
- StoreMBUnknowns (master_ptr->u, &s->moles, eos_p->coef * master_ptr->coef, &s->dg);
- }
- }
-
- return true;
- }
-
- bool ModelEngine::BuildMBSums()
- {
- //
- // Function builds lists sum_mb1 and sum_mb2 that describe how to sum molalities
- // to calculate mass balance sums, including activity of water, ionic strength,
- // charge balance, and alkalinity.
- //
-
- int i;
- //LDBLE *target;
- UnknownInfo *ui;
- int count = md->mb_unknowns->Count();
-
- for (i = 0; i < count; i++)
- {
- ui = (*md->mb_unknowns)[i];
- StoreMB(ui->source, &(ui->u->f), ui->coef);
-
- #ifdef DEBUG_MOHID
- if (d.debug_status)
- fprintf(d.buildjacobiansums_f, "BuildMBSums-1) i:%d %d %s %20.20e\n",
- i,
- ui->u->number,
- ui->u->name.CharPtr(),
- ui->coef);
- #endif
- }
-
- return true;
- }
-
- bool ModelEngine::BuildJacobianSums(int k)
- {
- //
- // Function builds lists sum_jacob1 and sum_jacob2 that describe how to sum molalities
- // to form jacobian.
- //
-
- int i, j, kk;
- int index;
- int count_g;
- LDBLE coef;
- //LDBLE *source, *target;
- UnknownInfo *mb_unknown;
- SpeciesDiffLayer *sdl_p;
- Unknown *u, *uprior, *u_kk;
-
- Species *s = gd->species_list[k];
-
- // Calculate jacobian coefficients for each mass balance equation
- for (i = 0; i < md->mb_unknowns->Count(); i++)
- {
- mb_unknown = (*md->mb_unknowns)[i];
-
- // Store d(moles) for a mass balance equation
-
- // initial solution only
- if (mb_unknown->u->type == SOLUTION_PHASE_BOUNDARY)
- continue;
-
- coef = mb_unknown->coef;
-
- StoreDN(k, mb_unknown->source, mb_unknown->u->number, coef, mb_unknown->gamma_source);
-
- #ifdef DEBUG_MOHID
- if (d.debug_status)
- fprintf(d.buildjacobiansums_f, "1) k:%d i:%d u_number:%d: %20.20e\n",
- k,
- i,
- mb_unknown->u->number,
- coef);
- #endif
-
- // Add extra terms for change in dg/dx in diffuse layer model
- if (s->type >= H2O || md->dl_type_x == NO_DL)
- continue;
- else if ((mb_unknown->u->type == MB || mb_unknown->u->type == MH || mb_unknown->u->type == MH2O) && md->state >= REACTION)
- {
- if (md->mass_oxygen_unknown != NULL)
- {
- // term for water, sum of all surfaces
- index = mb_unknown->u->number * (count_unknowns + 1) + md->mass_oxygen_unknown->number;
- StoreJacob(&s->tot_dh2o_moles, &arr[index], coef);
-
- #ifdef DEBUG_MOHID
- if (d.debug_status)
- fprintf(d.buildjacobiansums_f, "2) index:%d s_number:%d u_number:%d: %d %20.20e\n",
- index,
- s->number,
- md->mass_oxygen_unknown->number,
- coef);
- #endif
- }
-
- // terms for psi, one for each surface
- count_g = 0;
- for (j = 0; j < count_unknowns; j++)
- {
- u = (*unknown_list)[j];
- sdl_p = (*s->diff_layer)[count_g];
-
- if (u->type != SURFACE_CB)
- continue;
-
- index = mb_unknown->u->number * (count_unknowns + 1) + u->number;
- StoreJacob(&sdl_p->dx_moles, &arr[index], coef);
- count_g++;
-
- #ifdef DEBUG_MOHID
- if (d.debug_status)
- fprintf(d.buildjacobiansums_f, "3) index:%d %d %d %d %20.20e\n",
- index,
- mb_unknown->u->number,
- u->number,
- count_g,
- coef);
- #endif
-
- if (count_g >= md->use.sur_p->charge->Count())
- break;
- }
-
- // terms for related phases
- count_g = 0;
- for (j = 0; j < count_unknowns; j++)
- {
- u = (*unknown_list)[j];
- sdl_p = (*s->diff_layer)[count_g];
-
- if (u->type != SURFACE_CB)
- continue;
-
- uprior = (*unknown_list)[j - 1];
-
- // has related phase
- if (uprior->surface_comp->phase_name.IsEmpty())
- continue;
-
- // now find the related phase
- for (kk = count_unknowns - 1; kk >= 0; kk--)
- {
- u_kk = (*unknown_list)[kk];
-
- if (u_kk->type != PP)
- continue;
-
- if (u_kk->p->name == uprior->surface_comp->phase_name)
- break;
- }
-
- if (kk >= 0)
- {
- index = mb_unknown->u->number * (count_unknowns + 1) + (*unknown_list)[kk]->number;
- StoreJacob(&sdl_p->drelated_moles, &arr[index], coef);
-
- #ifdef DEBUG_MOHID
- if (d.debug_status)
- fprintf(d.buildjacobiansums_f, "4) index:%d %d %d %20.20e\n",
- index,
- mb_unknown->u->number,
- (*unknown_list)[kk]->number,
- coef);
- #endif
- }
-
- count_g++;
- if (count_g >= md->use.sur_p->charge->Count())
- break;
- }
- }
- else if (mb_unknown->u->type == SURFACE_CB)
- {
- count_g = 0;
-
- for (j = 0; j < count_unknowns; j++)
- {
- u = (*unknown_list)[j];
- sdl_p = (*s->diff_layer)[count_g];
-
- if (u->type != SURFACE_CB)
- continue;
-
- if (mb_unknown->u->number == u->number)
- {
- index = mb_unknown->u->number * (count_unknowns + 1) + u->number;
- StoreJacob(&sdl_p->dx_moles, &arr[index], coef);
-
- uprior = (*unknown_list)[j - 1];
-
- // term for related phase
- // has related phase
- if (!uprior->surface_comp->phase_name.IsEmpty())
- {
- // now find the related phase
- for (kk = count_unknowns - 1; kk >= 0; kk--)
- {
- u_kk = (*unknown_list)[kk];
-
- if (u_kk->type != PP)
- continue;
-
- if (u_kk->p->name == uprior->surface_comp->phase_name)
- break;
- }
-
- if (kk >= 0)
- {
- index = mb_unknown->u->number * (count_unknowns + 1) + u_kk->number;
- StoreJacob(&sdl_p->drelated_moles, &arr[index], coef);
- }
- }
-
- if (md->mass_oxygen_unknown != NULL)
- {
- // term for water, for same surfaces
- index = mb_unknown->u->number * (count_unknowns + 1) + md->mass_oxygen_unknown->number;
- StoreJacob(&sdl_p->dh2o_moles, &arr[index], coef);
- }
- break;
- }
-
- count_g++;
-
- if (count_g >= md->use.sur_p->charge->Count())
- break;
- }
- }
- }
-
- return true;
- }
-
- bool ModelEngine::WriteMBForSpeciesList(int n)
- {
- //
- // Sets up data to add to species_list
- // Original secondary redox states are retained
- //
-
- int i;
- char *ptr, token[DEFAULT_STRING_LENGTH];
- ElementOfSpecies *new_eos;
-
- r_temp->Reset();
- r_temp->AddTRXN(gd->species_list[n]->rxn_s, 1.0, false);
-
- //Copy to elt_list
- eos_list->Clear();
- parent_count = 0;
-
- ReactionToken *r;
- for (i = 1; i < r_temp->token_list->Count(); i++)
- {
- r = r_temp->token_list->Element(i);
-
- if (r->s->secondary == NULL)
- {
- r->s->primary->e->name.Copy(token);
- ptr = &token[0];
- GetSecondaryElementsInSpecies(&ptr, r->coef);
- }
- else
- {
- r->s->secondary->e->name.Copy(token);
- ptr = &token[0];
- GetSecondaryElementsInSpecies(&ptr, r->coef);
- }
- }
-
- ElementOfSpecies *eos;
- for (i = 0; i < eos_list->Count(); i++)
- {
- eos = (*eos_list)[i];
- if (eos->e->name == "O(-2)")
- {
- new_eos = eos_list->AddNew();
-
- new_eos->e = gd->e_h_one;
- new_eos->coef = eos->coef * 2;
- }
- }
-
- CombineElements();
-
- eos_list->CopyTo(gd->species_list[n]->eos_sys_total);
-
- return true;
- }
-
- bool ModelEngine::BuildSpeciesList(int n)
- {
- //
- // Builds a list that includes an entry for each master species in each
- // secondary reaction. Used for summing species of each element and
- // printing results.
- //
-
- int j;
- SpeciesInfo *new_s;
- ElementOfSpecies *eos_p;
- Species *s = gd->species_list[n];
-
- // Treat species made only with H+, e-, and H2O specially
- if (IsSpecial(s))
- {
- new_s = md->species_info_list->AddNew();
-
- new_s->master_s = gd->s_hplus;
- new_s->s = s;
- new_s->coef = 0.0;
-
- return true;
- }
-
- // Treat exchange species specially
- if (s->type == EX)
- {
- if (s->primary != NULL)
- return true; // master species has md->zero molality
-
- for (j = 0; j < eos_list->Count(); j++)
- {
- eos_p = (*eos_list)[j];
-
- if (eos_p->e->master->s->type != EX)
- continue;
-
- new_s = md->species_info_list->AddNew();
-
- new_s->master_s = eos_p->e->master->s;
- new_s->s = s;
- new_s->coef = eos_p->e->master->coef * eos_p->coef;
- }
-
- return true;
- }
-
- // Treat surface species specially
- if (s->type == SURF_PSI)
- return true;
-
- if (s->type == SURF)
- {
- for (j = 0; j < eos_list->Count(); j++)
- {
- eos_p = (*eos_list)[j];
-
- if (eos_p->e->master->s->type != SURF)
- continue;
-
- new_s = md->species_info_list->AddNew();
-
- new_s->master_s = eos_p->e->master->s;
- new_s->s = s;
- new_s->coef = eos_p->e->master->coef * eos_p->coef;
- }
-
- return true;
- }
-
- // Other aqueous species
- Master *master_ptr;
- for (j = 0; j < eos_list->Count(); j++)
- {
- eos_p = (*eos_list)[j];
-
- if (IsSpecial(eos_p->e->master->s))
- continue;
-
- if (eos_p->e->master->s->secondary != NULL)
- master_ptr = eos_p->e->master->s->secondary;
- else
- master_ptr = eos_p->e->master->s->primary;
-
- new_s = md->species_info_list->AddNew();
-
- new_s->master_s = master_ptr->s;
- new_s->s = s;
- new_s->coef = master_ptr->coef * eos_p->coef;
- }
-
- return true;
- }
-
- bool ModelEngine::StoreMB(LDBLE * source, LDBLE * target, LDBLE coef)
- {
- //
- // Adds item to list sum_mb1 or sum_mb2
- // If coef is 1.0, adds to sum_mb1, which does not require a multiply
- // else, adds to sum_mb2, which will multiply by coef
- //
-
- STCoef *new_item;
-
- if (Equal(coef, (LDBLE)1.0, (LDBLE)TOLERANCE))
- {
- new_item = md->sum_mb1->AddNew();
-
- new_item->source = source;
- new_item->target = target;
- }
- else
- {
- new_item = md->sum_mb2->AddNew();
-
- new_item->source = source;
- new_item->target = target;
- new_item->coef = coef;
- }
-
- return true;
- }
-
- bool ModelEngine::WritePhaseSysTotal(int n)
- {
- //
- // Sets up data to add to species_list
- // Original secondary redox states are retained
- //
-
- int i;
- ReactionToken *r;
- char *ptr, token[DEFAULT_STRING_LENGTH];
- Phase *p = gd->phase_list[n];
-
- // Start with secondary reaction
- r_temp->Reset();
- r_temp->AddTRXN(p->rxn_s, 1.0, false);
-
- // Copy to elt_list
- eos_list->Clear();
- parent_count = 0;
-
- for (i = 1; i < r_temp->token_list->Count(); i++)
- {
- r = r_temp->token_list->Element(i);
-
- if (r->s->secondary == NULL)
- r->s->primary->name.Copy(token); //master and element have the same name (change made in the load of database)
- else
- r->s->secondary->name.Copy(token);
-
- ptr = &token[0];
- GetSecondaryElementsInSpecies(&ptr, r->coef);
- }
-
- ElementOfSpecies *eos, *new_eos;
- for (i = 0; i < eos_list->Count(); i++)
- {
- eos = (*eos_list)[i];
-
- if (eos->e->name == "O(-2)")
- {
- new_eos = eos_list->AddNew();
-
- new_eos->e = gd->e_h_one;
- new_eos->coef = eos->coef * 2;
- }
- }
-
- CombineElements();
-
- eos_list->CopyTo(p->eos_sys_total);
-
- return true;
- }
-
- bool ModelEngine::BuildSolutionPhaseBoundaries()
- {
- //
- // Build into sums the logic to calculate inverse saturation indices for
- // solution phase boundaries
- //
-
- int i, j;
- Master *master_ptr;
- ReactionToken *rxn_ptr;
- Unknown *x;
-
- if (md->solution_phase_boundary_unknown == NULL)
- return true;
-
- // Calculate inverse saturation index
- for (i = 0; i < count_unknowns; i++)
- {
- x = (*unknown_list)[i];
-
- if (x->type != SOLUTION_PHASE_BOUNDARY)
- continue;
-
- StoreMB(&x->p->lk, &x->f, 1.0);
- StoreMB(&x->si, &x->f, 1.0);
-
- if (!x->p->in)
- return false;
-
- for (j = 1; j < x->p->rxn_x->token_list->Count(); j++)
- {
- rxn_ptr = x->p->rxn_x->token_list->Element(j);
- StoreMB(&rxn_ptr->s->la, &x->f, -rxn_ptr->coef);
- }
- }
-
- // Put coefficients into array
- for (i = 0; i < count_unknowns; i++)
- {
- x = (*unknown_list)[i];
-
- if (x->type != SOLUTION_PHASE_BOUNDARY)
- continue;
-
- for (j = 1; j < x->p->rxn_x->token_list->Count(); j++)
- {
- rxn_ptr = x->p->rxn_x->token_list->Element(j);
-
- if (rxn_ptr->s->secondary != NULL && rxn_ptr->s->secondary->in)
- master_ptr = rxn_ptr->s->secondary;
- else
- master_ptr = rxn_ptr->s->primary;
-
- if (master_ptr->u == NULL)
- continue;
-
- StoreJacob0(x->number, master_ptr->u->number, rxn_ptr->coef);
- }
- }
-
- return true;
- }
-
-
-
- bool ModelEngine::BuildMinExch()
- {
- //
- // Defines proportionality factor between mineral and exchanger to jacob0
- //
-
- int i, j, k, jj;
- int row;
- ExchComp *comp_ptr;
- Master *master_ptr, *m0;
- Unknown *unknown_ptr, *u_j, *u_k;
- ElementOfSpecies *eos_p;
- LDBLE coef;
-
- if (md->use.exc_p == NULL)
- return true;
-
- if (!md->use.exc_p->related_phases)
- return true;
-
- for (i = 0; i < md->use.exc_p->comps->Count(); i++)
- {
- comp_ptr = (*md->use.exc_p->comps)[i];
-
- if (comp_ptr->phase_name.IsEmpty())
- continue;
-
- // find unknown number
- for (j = count_unknowns - 1; j >= 0; j--)
- {
- u_j = (*unknown_list)[j];
- m0 = (*u_j->master)[0];
-
- if (u_j->type != EXCH)
- continue;
-
- if (m0 == comp_ptr->master)
- break;
- }
-
- for (k = count_unknowns - 1; k >= 0; k--)
- {
- u_k = (*unknown_list)[k];
-
- if (u_k->type != PP)
- continue;
-
- if (u_k->p->name == comp_ptr->phase_name)
- break;
- }
-
- if (j == -1)
- return false;
-
- if (k == -1)
- continue;
-
- // Build jacobian
-
- // charge balance
- StoreJacob0(md->charge_balance_unknown->number, u_k->number, comp_ptr->formula_z * comp_ptr->phase_proportion);
- StoreSumDeltas(&delta[k], &md->charge_balance_unknown->delta, -comp_ptr->formula_z * comp_ptr->phase_proportion);
-
- // mole balance balance
- eos_list->Clear();
- parent_count = 0;
-
- CopyToTempEOSList(comp_ptr->formula_totals, 1.0);
- ChangeHydrogenInEOSList(0);
-
- for (jj = 0; jj < eos_list->Count(); jj++)
- {
- eos_p = (*eos_list)[jj];
- master_ptr = (*eos_list)[jj]->e->primary;
-
- if (!master_ptr->in)
- master_ptr = master_ptr->s->secondary;
-
- if (master_ptr == NULL)
- return false;
-
- if (master_ptr->s->type == EX)
- {
- if (!Equal(u_j->moles, u_k->moles * eos_p->coef * comp_ptr->phase_proportion, (LDBLE)5.0 * md->convergence_tolerance))
- u_j->moles = u_k->moles * eos_p->coef * comp_ptr->phase_proportion;
- }
-
- coef = eos_p->coef;
-
- if (master_ptr->s == gd->s_hplus)
- {
- row = md->mass_hydrogen_unknown->number;
- unknown_ptr = md->mass_hydrogen_unknown;
- }
- else if (master_ptr->s == gd->s_h2o)
- {
- row = md->mass_oxygen_unknown->number;
- unknown_ptr = md->mass_oxygen_unknown;
- }
- else
- {
- row = master_ptr->u->number;
- unknown_ptr = master_ptr->u;
- }
-
- StoreJacob0 (row, u_k->number, coef * comp_ptr->phase_proportion);
- StoreSumDeltas(&delta[k], &unknown_ptr->delta, -coef * comp_ptr->phase_proportion);
- }
- }
-
- return true;
- }
-
- bool ModelEngine::BuildMinSurface()
- {
- //
- // Defines proportionality factor between mineral and surface to jacob0
- //
-
- int i, j, k, jj, row;
- List<ElementOfSpecies> *next_elt;
- ElementOfSpecies *eos_p;
- SurfaceComp *comp_ptr, *sc_p;
- Unknown *unknown_ptr, *u_j, *u_k, *u;
- Master *master_ptr, *m;
- LDBLE coef;
-
- if (md->use.sur_p == NULL)
- return true;
-
- if (!md->use.sur_p->related_phases)
- return true;
-
- for (i = 0; i < md->use.sur_p->comps->Count(); i++)
- {
- sc_p = (*md->use.sur_p->comps)[i];
-
- if (sc_p->phase_name.IsEmpty())
- continue;
-
- // find unknown number
- for (j = count_unknowns - 1; j >= 0; j--)
- {
- u_j = (*unknown_list)[j];
- m = (*u_j->master)[0];
-
- if (u_j->type != SURFACE)
- continue;
-
- if (m == sc_p->master)
- break;
- }
-
- for (k = count_unknowns - 1; k >= 0; k--)
- {
- u_k = (*unknown_list)[k];
-
- if (u_k->type != PP)
- continue;
-
- if (u_k->p->name == sc_p->phase_name)
- break;
- }
-
- if (j == -1)
- return false;
-
- if (k == -1)
- continue;
-
- comp_ptr = u_j->surface_comp;
-
- if (md->use.sur_p->type == CD_MUSIC)
- {
- // Add formula for CD_MUSIC
- next_elt = comp_ptr->formula_totals;
- }
- else
- {
- // Add master species for non CD_MUSIC
- next_elt = (*u_j->master)[0]->s->eos_list;
- }
-
- // update grams == moles in this case
- u = (*unknown_list)[j + 1];
- if (j < count_unknowns - 1 && u->type == SURFACE_CB)
- StoreSumDeltas (&delta[k], &u->related_moles, -1.0);
-
- // charge balance
- StoreJacob0 (md->charge_balance_unknown->number, u_k->number, comp_ptr->formula_z * comp_ptr->phase_proportion);
- StoreSumDeltas (&delta[k], &md->charge_balance_unknown->delta, -comp_ptr->formula_z * comp_ptr->phase_proportion);
-
- eos_list->Clear();
- parent_count = 0;
-
- CopyToTempEOSList(next_elt, 1.0);
- ChangeHydrogenInEOSList(0);
-
- for (jj = 0; jj < eos_list->Count(); jj++)
- {
- eos_p = (*eos_list)[jj];
-
- master_ptr = eos_p->e->primary;
- if (!master_ptr->in)
- master_ptr = master_ptr->s->secondary;
-
- if (master_ptr == NULL)
- return false;
-
- if (master_ptr->s->type == SURF)
- {
- if (!Equal(u_j->moles, u_k->moles * eos_p->coef * comp_ptr->phase_proportion, (LDBLE)5.0 * md->convergence_tolerance))
- u_j->moles = u_k->moles * eos_p->coef * comp_ptr->phase_proportion;
- }
-
- coef = eos_p->coef;
-
- if (master_ptr->s == gd->s_hplus)
- {
- row = md->mass_hydrogen_unknown->number;
- unknown_ptr = md->mass_hydrogen_unknown;
- }
- else if (master_ptr->s == gd->s_h2o)
- {
- row = md->mass_oxygen_unknown->number;
- unknown_ptr = md->mass_oxygen_unknown;
- }
- else
- {
- row = master_ptr->u->number;
- unknown_ptr = master_ptr->u;
- }
-
- StoreJacob0(row, u_k->number, coef * comp_ptr->phase_proportion);
- StoreSumDeltas(&delta[k], &unknown_ptr->delta, -coef * comp_ptr->phase_proportion);
- }
- }
-
- return true;
- }
-
- bool ModelEngine::BuildGasPhase()
- {
- //
- // Put coefficients into lists to sum iaps to test for equilibrium
- // Put coefficients into lists to build jacobian for
- // sum of partial pressures equation and
- // mass balance equations for elements contained in gases
- //
-
- int i, j, k;
- int row, col;
- Master *master_ptr;
- ReactionToken *rxn_ptr;
- GasComp *gas_comp_ptr;
- Phase *phase_ptr;
- Unknown *unknown_ptr;
- LDBLE coef, coef_elt;
- ElementOfSpecies *eos_p;
-
- if (md->gas_unknown == NULL)
- return true;
-
- for (i = 0; i < md->use.gas_p->comps->Count(); i++)
- {
- // Determine elements in gas component
- eos_list->Clear();
- parent_count = 0;
-
- gas_comp_ptr = (*md->use.gas_p->comps)[i];
- phase_ptr = gas_comp_ptr->phase;
-
- if (phase_ptr->rxn_x->token_list->Count() <= 0)
- continue;
-
- CopyToTempEOSList(phase_ptr->eos_list, 1.0);
- ChangeHydrogenInEOSList(0);
-
- // Build mass balance sums for each element in gas
-
- // All elements in gas
- for (j = 0; j < eos_list->Count(); j++)
- {
- eos_p = (*eos_list)[j];
-
- unknown_ptr = NULL;
- if (eos_p->e->name == "H")
- unknown_ptr = md->mass_hydrogen_unknown;
- else if (eos_p->e->name == "O")
- unknown_ptr = md->mass_oxygen_unknown;
- else
- {
- if (eos_p->e->primary->in)
- unknown_ptr = eos_p->e->primary->u;
- else if (eos_p->e->primary->s->secondary != NULL)
- unknown_ptr = eos_p->e->primary->s->secondary->u;
- }
-
- if (unknown_ptr != NULL)
- {
- coef = eos_p->coef;
- StoreMB(&(gas_comp_ptr->phase->moles_x), &(unknown_ptr->f), coef);
- }
- }
-
- if (md->use.gas_p->type == PRESSURE)
- {
- // Total pressure of gases
- StoreMB(&(gas_comp_ptr->phase->p_soln_x), &(md->gas_unknown->f), 1.0);
- }
-
- // Build jacobian sums for mass balance equations
- for (j = 0; j < eos_list->Count(); j++)
- {
- eos_p = (*eos_list)[j];
-
- unknown_ptr = NULL;
-
- if (eos_p->e->name == "H")
- unknown_ptr = md->mass_hydrogen_unknown;
- else if (eos_p->e->name == "O")
- unknown_ptr = md->mass_oxygen_unknown;
- else
- {
- if (eos_p->e->primary->in)
- unknown_ptr = eos_p->e->primary->u;
- else if (eos_p->e->primary->s->secondary != NULL)
- unknown_ptr = eos_p->e->primary->s->secondary->u;
- }
-
- if (unknown_ptr == NULL)
- continue;
-
- row = unknown_ptr->number * (count_unknowns + 1);
- coef_elt = eos_p->coef;
-
- for (k = 1; k < phase_ptr->rxn_x->token_list->Count(); k++)
- {
- rxn_ptr = phase_ptr->rxn_x->token_list->Element(k);
-
- if (rxn_ptr->s->secondary != NULL && rxn_ptr->s->secondary->in)
- master_ptr = rxn_ptr->s->secondary;
- else
- master_ptr = rxn_ptr->s->primary;
-
- if (master_ptr == NULL)
- return false;
-
- if (master_ptr->u == NULL)
- continue;
-
- if (!master_ptr->in)
- return false;
-
- col = master_ptr->u->number;
- coef = coef_elt * rxn_ptr->coef;
- StoreJacob (&(gas_comp_ptr->phase->moles_x), &(arr[row + col]), coef);
- }
-
- if (md->use.gas_p->type == PRESSURE)
- {
- // derivative wrt total moles of gas
- StoreJacob(&(gas_comp_ptr->phase->fraction_x), &(arr[row + md->gas_unknown->number]), coef_elt);
- }
- }
-
- // Build jacobian sums for sum of partial pressures equation
- if (md->use.gas_p->type != PRESSURE)
- continue;
-
- unknown_ptr = md->gas_unknown;
- row = unknown_ptr->number * (count_unknowns + 1);
-
- for (k = 1; k < phase_ptr->rxn_x->token_list->Count(); k++)
- {
- rxn_ptr = phase_ptr->rxn_x->token_list->Element(k);
-
- if (rxn_ptr->s->secondary != NULL && rxn_ptr->s->secondary->in)
- master_ptr = rxn_ptr->s->secondary;
- else
- master_ptr = rxn_ptr->s->primary;
-
- if (master_ptr->u == NULL)
- continue;
-
- if (!master_ptr->in)
- return false;
-
- col = master_ptr->u->number;
- coef = rxn_ptr->coef;
- StoreJacob (&(gas_comp_ptr->phase->p_soln_x), &(arr[row + col]), coef);
- }
- }
-
- return true;
- }
-
- bool ModelEngine::BuildSSAssemblage()
- {
- //
- // Put coefficients into lists to sum iaps to test for equilibrium
- // Put coefficients into lists to build jacobian for
- // mass action equation for component
- // mass balance equations for elements contained in solid solutions
- //
-
- int i, j, k, l, m;
- bool stop;
- int row, col;
- Master *master_ptr, *m_2, *m0;
- ReactionToken *rxn_ptr;
- Unknown *u, *u2;
- SS *s_s_ptr, *s_s_ptr_old;
- //SSComp *ssc_p;
- ElementOfSpecies *eos_p;
- char token[MAX_LENGTH];
- char *ptr;
-
- if (md->s_s_unknown == NULL)
- return true;
-
- s_s_ptr_old = NULL;
- col = 0;
-
- for (i = 0; i < count_unknowns; i++)
- {
- u = (*unknown_list)[i];
-
- if (u->type != S_S_MOLES)
- continue;
-
- s_s_ptr = u->s_s;
-
- if (s_s_ptr != s_s_ptr_old)
- {
- col = u->number;
- s_s_ptr_old = s_s_ptr;
- …
Large files files are truncated, but you can click here to view the full file