GeochemistrySpeciesSwapper Class Reference

Class to swap basis species with equilibrium species. More...

#include <GeochemistrySpeciesSwapper.h>

Public Member Functions
	GeochemistrySpeciesSwapper (unsigned basis_size, Real stoi_tol)
bool	operator== (const GeochemistrySpeciesSwapper &rhs) const
void	checkSwap (const ModelGeochemicalDatabase &mgd, const std::string &replace_this, const std::string &with_this)
	Check that replacing the named basis species with the named equilibrium species is valid. More...
void	checkSwap (const ModelGeochemicalDatabase &mgd, unsigned basis_index_to_replace, unsigned eqm_index_to_insert)
	Check that replacing the indexed basis species with the indexed equilibrium species is valid. More...
void	performSwap (ModelGeochemicalDatabase &mgd, const std::string &replace_this, const std::string &with_this)
	Check that replacing the named basis species with the named equilibrium species is valid, and then perform this swap by altering the mgd data structure. More...
void	performSwap (ModelGeochemicalDatabase &mgd, unsigned basis_index_to_replace, unsigned eqm_index_to_insert)
	Check that replacing the indexed basis species with the indexed equilibrium species is valid, and then perform this swap by altering the mgd data structure. More...
void	performSwap (ModelGeochemicalDatabase &mgd, DenseVector< Real > &bulk_composition, const std::string &replace_this, const std::string &with_this)
	Check that replacing the named basis species with the named equilibrium species is valid, and then perform this swap by altering the mgd data structure, and calculate the bulk composition of the new system. More...
void	performSwap (ModelGeochemicalDatabase &mgd, DenseVector< Real > &bulk_composition, unsigned basis_index_to_replace, unsigned eqm_index_to_insert)
	Check that replacing the indexed basis species with the indexed equilibrium species is valid, and then perform this swap by altering the mgd data structure, and calculate the bulk composition of the new system. More...
bool	findBestEqmSwap (unsigned basis_ind, const ModelGeochemicalDatabase &mgd, const std::vector< Real > &eqm_molality, bool minerals_allowed, bool gas_allowed, bool sorption_allowed, unsigned &best_eqm_species) const
	For the the given basis index, find the equilibrium index that it should be swapped with. More...

Private Member Functions
void	constructInverseMatrix (const ModelGeochemicalDatabase &mgd, unsigned basis_index_to_replace, unsigned eqm_index_to_insert)
	Construct the swap matrix, and its inverse, that describes the swap between the indexed basis species and the indexed equilibrium index. More...
void	alterMGD (ModelGeochemicalDatabase &mgd, unsigned basis_index_to_replace, unsigned eqm_index_to_insert)
	Modify the ModelGeochemicalDatabase mgd to swap the indexed basis species and the indexed equilibrium index. More...
void	alterBulkComposition (unsigned basis_size, DenseVector< Real > &bulk_composition) const
	Calculate the bulk composition in the new basis based on the bulk composition in the original basis. More...

Private Attributes
Real	_stoi_tol
	tolerance for matrix inversion and stoichiometric coefficients More...
DenseMatrix< Real >	_swap_matrix
	swap matrix More...
DenseMatrix< Real >	_true_swap_matrix
	Before _swap_matrix.svd (a non-const method) is called, we copy _true_swap_matrix = _swap_matrix. More...
DenseMatrix< Real >	_inv_swap_matrix
	inverse of swap matrix More...
DenseVector< Real >	_swap_sigma
	used in SVD decomposition of swap matrix More...
DenseMatrix< Real >	_swap_U
	used in SVD decomposition of swap matrix More...
DenseMatrix< Real >	_swap_VT
	used in SVD decomposition of swap matrix More...

Detailed Description

Class to swap basis species with equilibrium species.

Definition at line 19 of file GeochemistrySpeciesSwapper.h.

Constructor & Destructor Documentation

GeochemistrySpeciesSwapper()

GeochemistrySpeciesSwapper::GeochemistrySpeciesSwapper	(	unsigned	basis_size,
		Real	stoi_tol
	)

Parameters

basis_size	Number of basis species in the geochemical model
stoi_tol	Swapping involves inverting matrices via a singular value decomposition. During this process: (1) if abs(singular value) < stoi_tol * L1norm(singular values), then the matrix is deemed singular (so the basis swap is deemed invalid); (2) if abs(any stoichiometric coefficient) < stoi_tol then it is set to zero. A suggested reasonable value is stoi_tol=1E-6

Definition at line 12 of file GeochemistrySpeciesSwapper.C.

  : _stoi_tol(stoi_tol),
    _swap_matrix(basis_size, basis_size),
    _true_swap_matrix(basis_size, basis_size),
    _inv_swap_matrix(basis_size, basis_size),
    _swap_sigma(basis_size),
    _swap_U(basis_size, basis_size),
    _swap_VT(basis_size, basis_size)
{
}

Member Function Documentation

alterBulkComposition()

void GeochemistrySpeciesSwapper::alterBulkComposition ( unsigned  basis_size, DenseVector< Real > &  bulk_composition  ) const

private

Calculate the bulk composition in the new basis based on the bulk composition in the original basis.

Parameters

basis_size	Size of basis, which is only used for error checking
[in/out]	bulk_composition Upon entry, the bulk composition in the original basis. Upon exit, the bulk composition in the basis after the swap

Definition at line 355 of file GeochemistrySpeciesSwapper.C.

Referenced by performSwap().

{
  if (bulk_composition.size() != basis_size)
    mooseError("GeochemistrySpeciesSwapper: bulk_composition has size ",
               bulk_composition.size(),
               " which differs from the basis size");
  DenseVector<Real> result;
  _inv_swap_matrix.vector_mult_transpose(result, bulk_composition);
  bulk_composition = result;
}

alterMGD()

void GeochemistrySpeciesSwapper::alterMGD ( ModelGeochemicalDatabase &  mgd, unsigned  basis_index_to_replace, unsigned  eqm_index_to_insert  )

private

Modify the ModelGeochemicalDatabase mgd to swap the indexed basis species and the indexed equilibrium index.

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
basis_index_to_replace	The index of the basis component that will be removed from the basis and added to the equilibrium species list
eqm_index_to_insert	The index of the equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 162 of file GeochemistrySpeciesSwapper.C.

Referenced by performSwap().

{
  const unsigned num_cols = mgd.basis_species_index.size();
  const unsigned num_rows = mgd.eqm_species_index.size();
  const unsigned num_temperatures = mgd.eqm_log10K.n();
  const unsigned num_redox = mgd.redox_stoichiometry.m();
  const unsigned kin_rows = mgd.kin_stoichiometry.m();
  const unsigned num_rate = mgd.kin_rate.size();
  const unsigned pro_ind_eqm =
      num_cols + eqm_index_to_insert; // index of the (eqm species that is being swapped into the
                                      // basis) in the kinetic-rate promoting_indices vectors

  // change names
  const std::string basis_name = mgd.basis_species_name[basis_index_to_replace];
  const std::string eqm_name = mgd.eqm_species_name[eqm_index_to_insert];
  mgd.basis_species_name[basis_index_to_replace] = eqm_name;
  mgd.basis_species_index.erase(basis_name);
  mgd.basis_species_index[eqm_name] = basis_index_to_replace;
  mgd.eqm_species_name[eqm_index_to_insert] = basis_name;
  mgd.eqm_species_index.erase(eqm_name);
  mgd.eqm_species_index[basis_name] = eqm_index_to_insert;

  // flag indicating whether the redox_lhs is the equilibrium species that is being put into the
  // basis
  const bool redox_lhs_going_to_basis = (eqm_name == mgd.redox_lhs);
  mooseAssert(
      !redox_lhs_going_to_basis,
      "Should not be able to swap the redox LHS into the basis"); // At present, it is impossible to
                                                                  // swap redox_lhs into the basis,
                                                                  // since redox_lhs cannot be an
                                                                  // equilibrium species, because it
                                                                  // is explicitly excluded in
                                                                  // PertinentGeochemicalSystem when
                                                                  // building the secondary species.
  /*
  In past versions of the code, it was possible to swap the redox LHS into the basis.
  If you want to make it possible again then mgd needs to be modified as follows
  if (redox_lhs_going_to_basis)
  {
    // need to make the redox_lhs equal to the species that is being taken from the basis
    mgd.redox_lhs = basis_name;
    for (unsigned red = 0; red < num_redox; ++red)
    {
      const Real alpha = mgd.eqm_stoichiometry(
          eqm_index_to_insert, basis_index_to_replace); // must be nonzero due to valid swap
      const Real alpha_r = mgd.redox_stoichiometry(red, basis_index_to_replace);
      mooseAssert(alpha != alpha_r,
                  "Cannot swap equilibrium species "
                      << eqm_name << " with basis species " << basis_name
                      << " because a redox reaction would result in 0 <-> 1");
      const Real coef = 1.0 / (alpha - alpha_r);
      for (unsigned i = 0; i < num_cols; ++i)
        mgd.redox_stoichiometry(red, i) = coef * (mgd.redox_stoichiometry(red, i) -
                                                  mgd.eqm_stoichiometry(eqm_index_to_insert, i));
      mgd.redox_stoichiometry(red, basis_index_to_replace) = 0.0;
      for (unsigned t = 0; t < num_temperatures; ++t)
        mgd.redox_log10K(red, t) =
            coef * (mgd.redox_log10K(red, t) - mgd.eqm_log10K(eqm_index_to_insert, t));
    }
  }
  */

  // record that the swap has occurred
  if (mgd.swap_to_original_basis.n() == 0)
    mgd.swap_to_original_basis = _true_swap_matrix;
  else
    mgd.swap_to_original_basis.left_multiply(_true_swap_matrix);
  mgd.have_swapped_out_of_basis.push_back(basis_index_to_replace);
  mgd.have_swapped_into_basis.push_back(eqm_index_to_insert);

  // express stoichiometry in new basis
  for (unsigned i = 0; i < num_cols; ++i)
    mgd.eqm_stoichiometry(eqm_index_to_insert, i) = 0.0;
  mgd.eqm_stoichiometry(eqm_index_to_insert, basis_index_to_replace) =
      1.0; // 1 * replace_this <-> 1 * replace_this
  mgd.eqm_stoichiometry.right_multiply(_inv_swap_matrix);
  for (unsigned i = 0; i < num_rows; ++i)
    for (unsigned j = 0; j < num_cols; ++j)
      if (std::abs(mgd.eqm_stoichiometry(i, j)) < _stoi_tol)
        mgd.eqm_stoichiometry(i, j) = 0.0;

  // if the redox_lhs is not changed by the swap, alter the redox stoichiometry
  if (!redox_lhs_going_to_basis)
    mgd.redox_stoichiometry.right_multiply(_inv_swap_matrix);
  for (unsigned red = 0; red < num_redox; ++red)
    for (unsigned j = 0; j < num_cols; ++j)
      if (std::abs(mgd.redox_stoichiometry(red, j)) < _stoi_tol)
        mgd.redox_stoichiometry(red, j) = 0.0;

  // express kinetic stoichiometry in new basis
  mgd.kin_stoichiometry.right_multiply(_inv_swap_matrix);
  for (unsigned i = 0; i < kin_rows; ++i)
    for (unsigned j = 0; j < num_cols; ++j)
      if (std::abs(mgd.kin_stoichiometry(i, j)) < _stoi_tol)
        mgd.kin_stoichiometry(i, j) = 0.0;

  // alter equilibrium constants for each temperature point
  for (unsigned t = 0; t < num_temperatures; ++t)
  {
    const Real log10k_eqm_species = mgd.eqm_log10K(eqm_index_to_insert, t);
    mgd.eqm_log10K(eqm_index_to_insert, t) =
        0.0; // 1 * replace_this <-> 1 * replace_this with log10K = 0
    for (unsigned row = 0; row < num_rows; ++row)
      mgd.eqm_log10K(row, t) -=
          mgd.eqm_stoichiometry(row, basis_index_to_replace) * log10k_eqm_species;

    // similar for the redox equations
    if (!redox_lhs_going_to_basis)
      for (unsigned red = 0; red < num_redox; ++red)
        mgd.redox_log10K(red, t) -=
            mgd.redox_stoichiometry(red, basis_index_to_replace) * log10k_eqm_species;

    // similar for kinetic
    for (unsigned kin = 0; kin < kin_rows; ++kin)
      mgd.kin_log10K(kin, t) -=
          mgd.kin_stoichiometry(kin, basis_index_to_replace) * log10k_eqm_species;
  }

  // swap the "is mineral" information
  const bool basis_was_mineral = mgd.basis_species_mineral[basis_index_to_replace];
  mgd.basis_species_mineral[basis_index_to_replace] = mgd.eqm_species_mineral[eqm_index_to_insert];
  mgd.eqm_species_mineral[eqm_index_to_insert] = basis_was_mineral;

  // swap the "is gas" information
  const bool basis_was_gas = mgd.basis_species_gas[basis_index_to_replace];
  mgd.basis_species_gas[basis_index_to_replace] = mgd.eqm_species_gas[eqm_index_to_insert];
  mgd.eqm_species_gas[eqm_index_to_insert] = basis_was_gas;

  // swap the "is transported" information
  const bool basis_was_transported = mgd.basis_species_transported[basis_index_to_replace];
  mgd.basis_species_transported[basis_index_to_replace] =
      mgd.eqm_species_transported[eqm_index_to_insert];
  mgd.eqm_species_transported[eqm_index_to_insert] = basis_was_transported;

  // swap the "charge" information
  const Real basis_charge = mgd.basis_species_charge[basis_index_to_replace];
  mgd.basis_species_charge[basis_index_to_replace] = mgd.eqm_species_charge[eqm_index_to_insert];
  mgd.eqm_species_charge[eqm_index_to_insert] = basis_charge;

  // swap the "radius" information
  const Real basis_radius = mgd.basis_species_radius[basis_index_to_replace];
  mgd.basis_species_radius[basis_index_to_replace] = mgd.eqm_species_radius[eqm_index_to_insert];
  mgd.eqm_species_radius[eqm_index_to_insert] = basis_radius;

  // swap the "molecular weight" information
  const Real basis_molecular_weight = mgd.basis_species_molecular_weight[basis_index_to_replace];
  mgd.basis_species_molecular_weight[basis_index_to_replace] =
      mgd.eqm_species_molecular_weight[eqm_index_to_insert];
  mgd.eqm_species_molecular_weight[eqm_index_to_insert] = basis_molecular_weight;

  // swap the "molecular volume" information
  const Real basis_molecular_volume = mgd.basis_species_molecular_volume[basis_index_to_replace];
  mgd.basis_species_molecular_volume[basis_index_to_replace] =
      mgd.eqm_species_molecular_volume[eqm_index_to_insert];
  mgd.eqm_species_molecular_volume[eqm_index_to_insert] = basis_molecular_volume;

  // No need to swap surface_complexation_info or gas_chi because they are not
  // bound to basis or eqm species

  // swap promoting indices in the rates
  for (unsigned r = 0; r < num_rate; ++r)
  {
    const Real promoting_index_of_original_basis =
        mgd.kin_rate[r].promoting_indices[basis_index_to_replace];
    mgd.kin_rate[r].promoting_indices[basis_index_to_replace] =
        mgd.kin_rate[r].promoting_indices[pro_ind_eqm];
    mgd.kin_rate[r].promoting_indices[pro_ind_eqm] = promoting_index_of_original_basis;
    const Real promoting_monod_index_of_original_basis =
        mgd.kin_rate[r].promoting_monod_indices[basis_index_to_replace];
    mgd.kin_rate[r].promoting_monod_indices[basis_index_to_replace] =
        mgd.kin_rate[r].promoting_monod_indices[pro_ind_eqm];
    mgd.kin_rate[r].promoting_monod_indices[pro_ind_eqm] = promoting_monod_index_of_original_basis;
    const Real promoting_half_saturation_of_original_basis =
        mgd.kin_rate[r].promoting_half_saturation[basis_index_to_replace];
    mgd.kin_rate[r].promoting_half_saturation[basis_index_to_replace] =
        mgd.kin_rate[r].promoting_half_saturation[pro_ind_eqm];
    mgd.kin_rate[r].promoting_half_saturation[pro_ind_eqm] =
        promoting_half_saturation_of_original_basis;
  }

  // swap progeny, if needed
  for (unsigned r = 0; r < num_rate; ++r)
  {
    if (mgd.kin_rate[r].progeny_index == pro_ind_eqm)
      mgd.kin_rate[r].progeny_index = basis_index_to_replace;
    else if (mgd.kin_rate[r].progeny_index == basis_index_to_replace)
      mgd.kin_rate[r].progeny_index = pro_ind_eqm;
  }
}

checkSwap() [1/2]

void GeochemistrySpeciesSwapper::checkSwap	(	const ModelGeochemicalDatabase &	mgd,
		const std::string &	replace_this,
		const std::string &	with_this
	)

Check that replacing the named basis species with the named equilibrium species is valid.

In performing this check, the "swap" matrix is inverted, so the check is somewhat expensive. Therefore, if the swap will be performed, assuming it's valid, the performSwap method should be used instead. A mooseError or mooseException will result if the swap is invalid

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
replace_this	The basis component that will be removed from the basis and added to the equilibrium species list
with_this	The equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 24 of file GeochemistrySpeciesSwapper.C.

Referenced by performSwap(), and TEST().

{
  if (replace_this == "H2O")
    mooseError("Cannot remove H2O from the basis");
  if (mgd.basis_species_index.count(replace_this) == 0)
    mooseError(replace_this + " is not in the basis, so cannot be removed from the basis");
  if (mgd.eqm_species_index.count(with_this) == 0)
    mooseError(with_this + " is not an equilibrium species, so cannot be "
                           "removed from the equilibrium species list");

  checkSwap(mgd, mgd.basis_species_index.at(replace_this), mgd.eqm_species_index.at(with_this));
}

checkSwap() [2/2]

void GeochemistrySpeciesSwapper::checkSwap	(	const ModelGeochemicalDatabase &	mgd,
		unsigned	basis_index_to_replace,
		unsigned	eqm_index_to_insert
	)

Check that replacing the indexed basis species with the indexed equilibrium species is valid.

In performing this check, the "swap" matrix is inverted, so the check is somewhat expensive. Therefore, if the swap will be performed, assuming it's valid, the performSwap method should be used instead. A mooseError or mooseException will result if the swap is invalid

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
basis_index_to_replace	The index of the basis component that will be removed from the basis and added to the equilibrium species list
eqm_index_to_insert	The index of the equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 40 of file GeochemistrySpeciesSwapper.C.

{
  const unsigned num_cols = mgd.basis_species_index.size();
  const unsigned num_rows = mgd.eqm_species_index.size();
  if (basis_index_to_replace == 0)
    mooseError("Cannot remove H2O from the basis");
  if (basis_index_to_replace >= num_cols)
    mooseError(basis_index_to_replace, " exceeds the number of basis species in the problem");
  if (eqm_index_to_insert >= num_rows)
    mooseError(eqm_index_to_insert, " exceeds the number of equilibrium species in the problem");
  if (mgd.surface_sorption_related[eqm_index_to_insert])
    mooseError(
        "Equilibrium species ",
        mgd.eqm_species_name[eqm_index_to_insert],
        " is involved in surface sorption so cannot be swapped into the basis.  If this is truly "
        "necessary, code enhancements will need to be made including: recording whether basis "
        "species are involved in surface sorption, including them in the surface-potential "
        "calculations, and carefully swapping surface-potential-modified equilibrium constants");

  constructInverseMatrix(mgd, basis_index_to_replace, eqm_index_to_insert);
}

constructInverseMatrix()

void GeochemistrySpeciesSwapper::constructInverseMatrix ( const ModelGeochemicalDatabase &  mgd, unsigned  basis_index_to_replace, unsigned  eqm_index_to_insert  )

private

Construct the swap matrix, and its inverse, that describes the swap between the indexed basis species and the indexed equilibrium index.

The inverse of the swap matrix is used in alterMGD to modify its datastructures to implement the swap. A mooseError or mooseException will result if the swap is invalid

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
basis_index_to_replace	The index of the basis component that will be removed from the basis and added to the equilibrium species list
eqm_index_to_insert	The index of the equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 65 of file GeochemistrySpeciesSwapper.C.

Referenced by checkSwap().

{
  const unsigned num_cols = mgd.basis_species_index.size();

  if (_swap_matrix.n() != num_cols)
    mooseError("GeochemistrySpeciesSwapper constructed with incorrect "
               "basis_species size");

  // This is a private method, called from checkSwap, so we know that all the inputs have been
  // sanity-checked.  The only way the swap could be invalid is that the swap matrix is not
  // invertible.  This could be due to the solve algorithm attempting to perform an invalid swap, so
  // only mooseExceptions are thrown below, instead of mooseError, to allow the solve algorithm a
  // chance to attempt another swap.

  // construct the swap matrix.  new_basis = swap_matrix * old_basis
  // We shove the equilibrium species into the "slot" of the basis species it is
  // replacing
  _swap_matrix.zero();
  for (unsigned i = 0; i < num_cols; ++i)
    _swap_matrix(i, i) = 1.0;
  for (unsigned i = 0; i < num_cols; ++i)
    _swap_matrix(basis_index_to_replace, i) = mgd.eqm_stoichiometry(eqm_index_to_insert, i);

  // record the swap matrix since the SVD trashes it
  _true_swap_matrix = _swap_matrix;

  // In order to invert _swap_matrix, perform the SVD: A = U * D * VT
  // Although every matrix has a SVD, this process might fail due to numerical errors.  Therefore
  // the following could be wrapped in a try-catch block, but I have never been able to trigger a
  // failure since _swap_matrix is always so well-formed due to all the checks on the database while
  // parsing it
  _swap_matrix.svd(_swap_sigma, _swap_U, _swap_VT);
  const Real l1norm = _swap_sigma.l1_norm();
  for (unsigned i = 0; i < num_cols; ++i)
    if (std::abs(_swap_sigma(i) / l1norm) < _stoi_tol)
      mooseException("Matrix is not invertible, which signals an invalid basis swap");

  // Find the inverse, which is VT^-1 * D^-1 * U^-1 = V * D^-1 * UT
  for (unsigned i = 0; i < num_cols; ++i)
    _swap_U.scale_column(i, 1.0 / _swap_sigma(i)); // (scale the columns of U)^T = D^-1 * UT
  for (unsigned i = 0; i < num_cols; ++i)
    for (unsigned j = 0; j < num_cols; ++j)
      _inv_swap_matrix(i, j) = _swap_VT(j, i);        // _inv_swap_matrix = V
  _inv_swap_matrix.right_multiply_transpose(_swap_U); // V * (scale the columns of U)^T
}

findBestEqmSwap()

bool GeochemistrySpeciesSwapper::findBestEqmSwap	(	unsigned	basis_ind,
		const ModelGeochemicalDatabase &	mgd,
		const std::vector< Real > &	eqm_molality,
		bool	minerals_allowed,
		bool	gas_allowed,
		bool	sorption_allowed,
		unsigned &	best_eqm_species
	)		const

For the the given basis index, find the equilibrium index that it should be swapped with.

Parameters

	basis_ind	index of basis species in mgd that we'd like to swap out of the basis
	mgd	the model's geochemical databgase
	eqm_molality	the molality of all the equilibrium species
	minerals_allowed	if false, best_eqm_species will not correspond to a mineral. Since equilibrium minerals have molality=0 it is very unlikely that best_eqm_species will correspond to a mineral, even if this flag is true.
	gases_allowed	if false, best_eqm_species will not correspond to a gas. Since equilibrium gases have molality=0 it is very unlikely that best_eqm_species will correspond to a gas, even if this flag is true.
	sorption_allowed	if false, best_eqm_species will not be involved in sorption
[out]	best_eqm_species	the index of the equilibrium species that should be swapped with basis_ind

Returns

true if a valid swap was found. if false, then best_eqm_species will be garbage

Definition at line 368 of file GeochemistrySpeciesSwapper.C.

Referenced by GeochemistryTimeDependentReactor::preSolveDump(), GeochemicalSolver::swapNeeded(), and TEST().

{
  const unsigned num_eqm = mgd.eqm_species_name.size();
  if (eqm_molality.size() != num_eqm)
    mooseError("Size of eqm_molality is ",
               eqm_molality.size(),
               " which is not equal to the number of equilibrium species ",
               num_eqm);
  if (basis_ind >= mgd.basis_species_name.size())
    mooseError("basis index ", basis_ind, " must be less than ", mgd.basis_species_name.size());
  best_eqm_species = 0;
  bool legitimate_swap_found = false;
  Real best_stoi = 0.0;
  for (unsigned j = 0; j < num_eqm; ++j)
  {
    if (mgd.eqm_stoichiometry(j, basis_ind) == 0.0)
      continue;
    if (!minerals_allowed && mgd.eqm_species_mineral[j])
      continue;
    if (!gas_allowed && mgd.eqm_species_gas[j])
      continue;
    if (!sorption_allowed && mgd.surface_sorption_related[j])
      continue;
    const Real stoi = std::abs(mgd.eqm_stoichiometry(j, basis_ind)) * eqm_molality[j];
    if (stoi >= best_stoi)
    {
      best_stoi = stoi;
      best_eqm_species = j;
      legitimate_swap_found = true;
    }
  }
  return legitimate_swap_found;
}

operator==()

bool GeochemistrySpeciesSwapper::operator== ( const GeochemistrySpeciesSwapper &  rhs ) const

inline

Definition at line 31 of file GeochemistrySpeciesSwapper.h.

  {
    return (_stoi_tol == rhs._stoi_tol) && (_swap_matrix.m() == rhs._swap_matrix.m());
  };

performSwap() [1/4]

void GeochemistrySpeciesSwapper::performSwap	(	ModelGeochemicalDatabase &	mgd,
		const std::string &	replace_this,
		const std::string &	with_this
	)

Check that replacing the named basis species with the named equilibrium species is valid, and then perform this swap by altering the mgd data structure.

A mooseError or mooseException will result if the swap is invalid

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
replace_this	The basis component that will be removed from the basis and added to the equilibrium species list
with_this	The equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 114 of file GeochemistrySpeciesSwapper.C.

Referenced by GeochemicalSystem::checkAndInitialize(), GeochemicalModelInterrogator::output(), performSwap(), GeochemicalSystem::performSwapNoCheck(), and TEST().

{
  // check the swap is valid (if not then mooseError or mooseException)
  // and if it's valid, create the inverse swap matrix
  checkSwap(mgd, replace_this, with_this);

  // perform the swap inside the MGD datastructure
  alterMGD(mgd, mgd.basis_species_index.at(replace_this), mgd.eqm_species_index.at(with_this));
}

performSwap() [2/4]

void GeochemistrySpeciesSwapper::performSwap	(	ModelGeochemicalDatabase &	mgd,
		unsigned	basis_index_to_replace,
		unsigned	eqm_index_to_insert
	)

Check that replacing the indexed basis species with the indexed equilibrium species is valid, and then perform this swap by altering the mgd data structure.

A mooseError or mooseException will result if the swap is invalid

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
basis_index_to_replace	The index of the basis component that will be removed from the basis and added to the equilibrium species list
eqm_index_to_insert	The index of the equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 127 of file GeochemistrySpeciesSwapper.C.

{
  // check the swap is valid (if not then mooseError or mooseException)
  // and if it's valid, create the inverse swap matrix
  checkSwap(mgd, basis_index_to_replace, eqm_index_to_insert);

  // perform the swap inside the MGD datastructure
  alterMGD(mgd, basis_index_to_replace, eqm_index_to_insert);
}

performSwap() [3/4]

void GeochemistrySpeciesSwapper::performSwap	(	ModelGeochemicalDatabase &	mgd,
		DenseVector< Real > &	bulk_composition,
		const std::string &	replace_this,
		const std::string &	with_this
	)

Check that replacing the named basis species with the named equilibrium species is valid, and then perform this swap by altering the mgd data structure, and calculate the bulk composition of the new system.

A mooseError or mooseException will result if the swap is invalid

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
[in/out]	bulk_composition Upon entry, the bulk composition in the original basis. Upon exit, the bulk composition in the basis after the swap
replace_this	The basis component that will be removed from the basis and added to the equilibrium species list
with_this	The equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 140 of file GeochemistrySpeciesSwapper.C.

{
  performSwap(mgd, replace_this, with_this);
  // compute the bulk composition expressed in the new basis
  alterBulkComposition(mgd.basis_species_index.size(), bulk_composition);
}

performSwap() [4/4]

void GeochemistrySpeciesSwapper::performSwap	(	ModelGeochemicalDatabase &	mgd,
		DenseVector< Real > &	bulk_composition,
		unsigned	basis_index_to_replace,
		unsigned	eqm_index_to_insert
	)

Check that replacing the indexed basis species with the indexed equilibrium species is valid, and then perform this swap by altering the mgd data structure, and calculate the bulk composition of the new system.

A mooseError or mooseException will result if the swap is invalid

Parameters

mgd	The ModelGeochemicalDatabase that holds all information regarding basis components, equilibrium species, stoichiometries, etc
[in/out]	bulk_composition Upon entry, the bulk composition in the original basis. Upon exit, the bulk composition in the basis after the swap
basis_index_to_replace	The index of the basis component that will be removed from the basis and added to the equilibrium species list
eqm_index_to_insert	The index of the equilibrium species that will be removed from the equilibrium species list and added to the basis component list

Definition at line 151 of file GeochemistrySpeciesSwapper.C.

{
  performSwap(mgd, basis_index_to_replace, eqm_index_to_insert);
  // compute the bulk composition expressed in the new basis
  alterBulkComposition(mgd.basis_species_index.size(), bulk_composition);
}

Member Data Documentation

_inv_swap_matrix

DenseMatrix<Real> GeochemistrySpeciesSwapper::_inv_swap_matrix

private

inverse of swap matrix

Definition at line 216 of file GeochemistrySpeciesSwapper.h.

Referenced by alterBulkComposition(), alterMGD(), and constructInverseMatrix().

_stoi_tol

Real GeochemistrySpeciesSwapper::_stoi_tol

private

tolerance for matrix inversion and stoichiometric coefficients

Definition at line 203 of file GeochemistrySpeciesSwapper.h.

Referenced by alterMGD(), constructInverseMatrix(), and operator==().

_swap_matrix

DenseMatrix<Real> GeochemistrySpeciesSwapper::_swap_matrix

private

swap matrix

Definition at line 206 of file GeochemistrySpeciesSwapper.h.

Referenced by constructInverseMatrix(), and operator==().

_swap_sigma

DenseVector<Real> GeochemistrySpeciesSwapper::_swap_sigma

private

used in SVD decomposition of swap matrix

Definition at line 219 of file GeochemistrySpeciesSwapper.h.

Referenced by constructInverseMatrix().

_swap_U

DenseMatrix<Real> GeochemistrySpeciesSwapper::_swap_U

private

used in SVD decomposition of swap matrix

Definition at line 222 of file GeochemistrySpeciesSwapper.h.

Referenced by constructInverseMatrix().

_swap_VT

DenseMatrix<Real> GeochemistrySpeciesSwapper::_swap_VT

private

used in SVD decomposition of swap matrix

Definition at line 225 of file GeochemistrySpeciesSwapper.h.

Referenced by constructInverseMatrix().

_true_swap_matrix

DenseMatrix<Real> GeochemistrySpeciesSwapper::_true_swap_matrix

private

Before _swap_matrix.svd (a non-const method) is called, we copy _true_swap_matrix = _swap_matrix.

It is necessary to record _swap_matrix in this fashion so that mgd.swap_to_original_basis can be modified using it

Definition at line 213 of file GeochemistrySpeciesSwapper.h.

Referenced by alterMGD(), and constructInverseMatrix().

---

The documentation for this class was generated from the following files:

• geochemistry/include/utils/GeochemistrySpeciesSwapper.h
• geochemistry/src/utils/GeochemistrySpeciesSwapper.C