GeochemicalSolver Class Reference

This class contains methods to solve the algebraic system in GeochemicalSystem. More...

#include <GeochemicalSolver.h>

Public Member Functions
	GeochemicalSolver (unsigned num_basis, unsigned num_kin, GeochemistryIonicStrength &is, Real abs_tol, Real rel_tol, unsigned max_iter, Real max_initial_residual, Real swap_threshold, unsigned max_swaps_allowed, const std::vector< std::string > &prevent_precipitation, Real max_ionic_strength, unsigned ramp_max_ionic_strength, bool evaluate_kin_always)
	Construct and check for sensible arguments. More...
void	solveSystem (GeochemicalSystem &egs, std::stringstream &ss, unsigned &tot_iter, Real &abs_residual, Real dt, DenseVector< Real > &mole_additions, DenseMatrix< Real > &dmole_additions)
	Solve the system. More...
void	setMaxInitialResidual (Real max_initial_residual)
	Set value for max_initial_residual. More...
Real	getMaxInitialResidual () const
	Get value for max_initial_residual. More...
void	setRampMaxIonicStrength (unsigned ramp_max_ionic_strength)
	Sets the value of _ramp_max_ionic_strength. More...
unsigned	getRampMaxIonicStrength () const
	Gets the value of _ramp_max_ionic_strength. More...

Private Member Functions
Real	computeResidual (const GeochemicalSystem &egs, DenseVector< Real > &residual, const DenseVector< Real > &mole_additions) const
	Builds the residual of the algebraic system. More...
void	solveAndUnderrelax (const GeochemicalSystem &egs, DenseMatrix< Real > &jacobian, DenseVector< Real > &new_mol) const
	Solves _jacobian * neg_change_mol = _residual for neg_change_mol, then performs an underrelaxation to get new_mol. More...
bool	swapNeeded (const GeochemicalSystem &egs, unsigned &swap_out_of_basis, unsigned &swap_into_basis, std::stringstream &ss) const
	Check if a basis swap is needed. More...
bool	reduceInitialResidual (GeochemicalSystem &egs, Real dt, DenseVector< Real > &mole_additions, DenseMatrix< Real > &dmole_additions)
	Progressively alter the initial-guess molalities for the algebraic system to attempt to reduce the residual. More...

Private Attributes
GeochemistryIonicStrength &	_is
	The ionic-strength calculator. More...
const unsigned	_num_basis
	Number of species in the basis. More...
const unsigned	_num_kin
	Number of kinetic species. More...
unsigned	_num_basis_in_algebraic_system
	Number of basis molalities (and potentially solvent water mass) in the algebraic system. More...
unsigned	_num_in_algebraic_system
	Number of unknowns (molalities and surface potentials) in the algebraic system. More...
DenseVector< Real >	_residual
	residual of the algebraic system we wish to solve More...
Real	_abs_residual
	L1 norm of residual. More...
DenseMatrix< Real >	_jacobian
	jacobian of the algebraic system More...
DenseVector< Real >	_new_mol
	the new molality after finding the solution of _jacobian * neg_change_mol = _residual More...
const Real	_abs_tol
	If the residual of the algebraic system falls below this value, the Newton process has converged. More...
const Real	_rel_tol
	If the residual of the algebraic system falls below this value times the initial residual, the Newton process has converged. More...
Real	_res0_times_rel
	_res0_times_rel = _rel_tol * initial residual More...
const unsigned	_max_iter
	maximum number of iterations allowed during an inner solve More...
Real	_max_initial_residual
	maximum desired initial residual More...
const Real	_swap_threshold
	If a basis molality < swap_threshold, we attempt to swap it out of the basis. More...
const unsigned	_max_swaps_allowed
	Maximum number of swaps allowed before the solve aborts. More...
const std::vector< std::string >	_prevent_precipitation
	The minerals named in this list can have positive saturation indices and will not precipitate. More...
const Real	_max_ionic_strength
	Maximum ionic strength allowed. More...
unsigned	_ramp_max_ionic_strength
	Number of iterations over which to increase the maximum ionic strength to _max_ionic_strength. More...
bool	_evaluate_kin_always
	When to compute the kinetic rates: if true then evaluate before every residual calculation, otherwise evaluate only at the start of the solve process (using, eg, the old molality values) More...
DenseVector< Real >	_input_mole_additions
	the mole_additions for the basis and kinetic species as specified in the argument of solveSystem More...
DenseMatrix< Real >	_input_dmole_additions
	d(mole_additions)/d(species) as specified in the argument of solveSystem More...

Detailed Description

This class contains methods to solve the algebraic system in GeochemicalSystem.

Definition at line 18 of file GeochemicalSolver.h.

Constructor & Destructor Documentation

GeochemicalSolver()

GeochemicalSolver::GeochemicalSolver	(	unsigned	num_basis,
		unsigned	num_kin,
		GeochemistryIonicStrength &	is,
		Real	abs_tol,
		Real	rel_tol,
		unsigned	max_iter,
		Real	max_initial_residual,
		Real	swap_threshold,
		unsigned	max_swaps_allowed,
		const std::vector< std::string > &	prevent_precipitation,
		Real	max_ionic_strength,
		unsigned	ramp_max_ionic_strength,
		bool	evaluate_kin_always
	)

Construct and check for sensible arguments.

Parameters

num_basis	Number of basis species in the system
mum_kin	Number of kinetic species in the system
is	The object to compute ionic strengths. GeochemicalSolver changes the maximum value of ionic strength as it progresses towards the solution
abs_tol	The Newton solution process is deemed to have converged if the L1 norm of the residual < abs_tol
rel_tol	The Newton solution process is deemed to have converged if the L1 norm of the residual < rel_tol * residual_0
max_iter	Only max_iter Newton iterations are allowed in solving the system
max_initial_residual	egs guesses the initial configuration, but sometimes this results in a very large initial residual. GeochemicalSolver will attempt to alter the initial guesses for molalities so that the initial residual is less than max_initial_residual. This parameter may be varied during runtime which might be useful in time-dependent simulations if the first residual is poor but the remainder are deemed likely to be OK
swap_threshold	If a basis molality < swap_threshold at the end of the Newton process, GeochemicalSolver attempts to swap it out of the basis
max_swaps_allowed	Maximum number of swaps allowed before the solve process aborts
prevent_precipitation	The minerals named in this vector will not be allowed to precipitate, even if their saturation indices are positive
max_ionic_strength	Maximum ionic strength ever allowed
ramp_max_ionic_strength	The maximum ionic strength is ramped from 0 to max_ionic_strength over the course of the first ramp_max_ionic_strength Newton iterations. This can help with convergence.
evaluate_kin_always	If true then evaluate the kinetic rates before each residual calculation, otherwise evaluate them only at the start of the solve process (which will be using the molalities, etc from the previous time-step)

Definition at line 13 of file GeochemicalSolver.C.

  : _is(is),
    _num_basis(num_basis),
    _num_kin(num_kin),
    _num_basis_in_algebraic_system(_num_basis),
    _num_in_algebraic_system(_num_basis + _num_kin),
    _residual(_num_in_algebraic_system),
    _abs_residual(0.0),
    _jacobian(_num_in_algebraic_system, _num_in_algebraic_system),
    _new_mol(_num_in_algebraic_system),
    _abs_tol(abs_tol),
    _rel_tol(rel_tol),
    _res0_times_rel(0.0),
    _max_iter(max_iter),
    _max_initial_residual(max_initial_residual),
    _swap_threshold(swap_threshold),
    _max_swaps_allowed(max_swaps_allowed),
    _prevent_precipitation(prevent_precipitation),
    _max_ionic_strength(max_ionic_strength),
    _ramp_max_ionic_strength(ramp_max_ionic_strength),
    _evaluate_kin_always(evaluate_kin_always),
    _input_mole_additions(_num_basis + _num_kin),
    _input_dmole_additions(_num_basis + _num_kin, _num_basis + _num_kin)
{
  if (_ramp_max_ionic_strength > _max_iter)
    mooseError("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter");
  if (max_initial_residual <= 0.0)
    mooseError("GeochemicalSolver: max_initial_residual must be positive");
  if (max_ionic_strength < 0.0)
    mooseError("GeochemicalSolver: max_ionic_strength must not be negative");
  if (abs_tol < 0.0)
    mooseError("GeochemicalSolver: abs_tol must not be negative");
  if (rel_tol < 0.0)
    mooseError("GeochemicalSolver: rel_tol must not be negative");
  if (rel_tol == 0.0 && abs_tol == 0.0)
    mooseError("GeochemicalSolver: either rel_tol or abs_tol must be positive");
}

Member Function Documentation

computeResidual()

Real GeochemicalSolver::computeResidual ( const GeochemicalSystem &  egs, DenseVector< Real > &  residual, const DenseVector< Real > &  mole_additions  ) const

private

Builds the residual of the algebraic system.

Parameters

egs	The geochemical system to use to compute the residual
residual	The residual components will be placed here
mole_additions	the increment of mole number of each basis species and kinetic species since the last timestep.

Returns

the L1 norm of residual

Definition at line 78 of file GeochemicalSolver.C.

Referenced by reduceInitialResidual(), and solveSystem().

{
  for (unsigned a_ind = 0; a_ind < _num_in_algebraic_system; ++a_ind)
    residual(a_ind) = egs.getResidualComponent(a_ind, mole_additions);
  return residual.l1_norm();
}

getMaxInitialResidual()

Real GeochemicalSolver::getMaxInitialResidual

(

)

const

Get value for max_initial_residual.

Definition at line 72 of file GeochemicalSolver.C.

{
  return _max_initial_residual;
}

getRampMaxIonicStrength()

unsigned GeochemicalSolver::getRampMaxIonicStrength

(

)

const

Gets the value of _ramp_max_ionic_strength.

Definition at line 470 of file GeochemicalSolver.C.

{
  return _ramp_max_ionic_strength;
}

reduceInitialResidual()

bool GeochemicalSolver::reduceInitialResidual ( GeochemicalSystem &  egs, Real  dt, DenseVector< Real > &  mole_additions, DenseMatrix< Real > &  dmole_additions  )

private

Progressively alter the initial-guess molalities for the algebraic system to attempt to reduce the residual.

Parameters

egs	The GeochemicalSystem we're trying to solve
dt	time-step size (used for determining kinetic rates)
mole_additions	the increment of mole number of each basis species and kinetic species since the last timestep. This may change during this function as kinetic rates depend on molalities
dmole_additions	d(mole_additions)/(molality and kinetic moles)

Definition at line 241 of file GeochemicalSolver.C.

Referenced by solveSystem().

{
  const Real initial_r = _abs_residual;

  // to get an indication of whether we should increase or decrease molalities in the algorithm
  // below, find the median of the original molalities
  const std::vector<Real> & original_molality = egs.getAlgebraicBasisValues();
  const std::vector<unsigned> mol_order =
      GeochemistrySortedIndices::sortedIndices(original_molality, false);
  const Real median_molality = original_molality[mol_order[_num_basis_in_algebraic_system / 2]];

  // get the index order of the residual vector (largest first): ignore residuals for surface
  // potentials (we're only using _num_basis_in_algebraic_system, not _num_in_algebraic_system)
  unsigned ind = 0;
  std::vector<unsigned> res_order(_num_basis_in_algebraic_system);
  std::iota(res_order.begin(), res_order.end(), ind++);
  std::sort(res_order.begin(),
            res_order.end(),
            [&](int i, int j) { return std::abs(_residual(i)) > std::abs(_residual(j)); });

  const std::vector<Real> & original_molality_and_pot = egs.getAlgebraicVariableValues();
  DenseVector<Real> new_molality_and_pot(original_molality_and_pot);
  for (const auto & a : res_order)
  {
    if (std::abs(_residual(a)) < _max_initial_residual)
      return false; // haven't managed to find a suitable new molality, and all remaining residual
                    // components are less than the cutoff, so cannot appreciably reduce from now
                    // on

    const Real multiplier = (original_molality_and_pot[a] > median_molality) ? 0.5 : 2.0;
    // try using the multiplier
    new_molality_and_pot(a) = multiplier * original_molality_and_pot[a];
    egs.setAlgebraicVariables(new_molality_and_pot);
    if (_evaluate_kin_always)
    {
      mole_additions = _input_mole_additions;
      dmole_additions = _input_dmole_additions;
      egs.addKineticRates(dt, mole_additions, dmole_additions);
    }
    _abs_residual = computeResidual(egs, _residual, mole_additions);
    if (_abs_residual < initial_r)
      return true; // success: found a new molality that decreases the initial |R|

    // the above approach did not decrease |R|, so try using 1/multiplier
    new_molality_and_pot(a) = original_molality_and_pot[a] / multiplier;
    egs.setAlgebraicVariables(new_molality_and_pot);
    if (_evaluate_kin_always)
    {
      mole_additions = _input_mole_additions;
      dmole_additions = _input_dmole_additions;
      egs.addKineticRates(dt, mole_additions, dmole_additions);
    }
    _abs_residual = computeResidual(egs, _residual, mole_additions);
    if (_abs_residual < initial_r)
      return true; // success: found a new molality that decreases the initial |R|

    // the new molalities did not decrease |R|, so revert to the original molality, and move to
    // next-biggest residual component
    new_molality_and_pot(a) = original_molality_and_pot[a];
    egs.setAlgebraicVariables(new_molality_and_pot);
    if (_evaluate_kin_always)
    {
      mole_additions = _input_mole_additions;
      dmole_additions = _input_dmole_additions;
      egs.addKineticRates(dt, mole_additions, dmole_additions);
    }
    _abs_residual = computeResidual(egs, _residual, mole_additions);
  }
  return false;
}

setMaxInitialResidual()

void GeochemicalSolver::setMaxInitialResidual

(

Real

max_initial_residual

)

Set value for max_initial_residual.

Definition at line 64 of file GeochemicalSolver.C.

Referenced by TEST().

{
  if (max_initial_residual <= 0.0)
    mooseError("GeochemicalSolver: max_initial_residual must be positive");
  _max_initial_residual = max_initial_residual;
}

setRampMaxIonicStrength()

void GeochemicalSolver::setRampMaxIonicStrength

(

unsigned

ramp_max_ionic_strength

)

Sets the value of _ramp_max_ionic_strength.

Definition at line 462 of file GeochemicalSolver.C.

Referenced by GeochemistryTimeDependentReactor::execute(), and GeochemistrySpatialReactor::execute().

{
  if (ramp_max_ionic_strength > _max_iter)
    mooseError("GeochemicalSolver: ramp_max_ionic_strength must be less than max_iter");
  _ramp_max_ionic_strength = ramp_max_ionic_strength;
}

solveAndUnderrelax()

void GeochemicalSolver::solveAndUnderrelax ( const GeochemicalSystem &  egs, DenseMatrix< Real > &  jacobian, DenseVector< Real > &  new_mol  ) const

private

Solves _jacobian * neg_change_mol = _residual for neg_change_mol, then performs an underrelaxation to get new_mol.

Parameters

egs	The geochemical system that we're trying to solve
jacobian	the jacobian of the system
new_mol	upon exit, this will be the new molality (and surface potential values, if any) values according to the underrelaxed Newton process

Definition at line 88 of file GeochemicalSolver.C.

Referenced by solveSystem().

{
  jacobian.lu_solve(_residual, new_mol);

  // at this point we want to do molality = molality - new_mol, but
  // Bethke recommends underrelaxation - probably want to do PETSc variational bounds in the
  // future
  Real one_over_delta = 1.0;
  const std::vector<Real> current_molality_and_pot = egs.getAlgebraicVariableValues();
  for (unsigned a_ind = 0; a_ind < _num_in_algebraic_system; ++a_ind)
    one_over_delta =
        std::max(one_over_delta, new_mol(a_ind) * 2.0 / current_molality_and_pot[a_ind]);
  for (unsigned a_ind = 0; a_ind < _num_in_algebraic_system; ++a_ind)
    new_mol(a_ind) = current_molality_and_pot[a_ind] - new_mol(a_ind) / one_over_delta;
}

solveSystem()

void GeochemicalSolver::solveSystem	(	GeochemicalSystem &	egs,
		std::stringstream &	ss,
		unsigned &	tot_iter,
		Real &	abs_residual,
		Real	dt,
		DenseVector< Real > &	mole_additions,
		DenseMatrix< Real > &	dmole_additions
	)

Solve the system.

Parameters

egs	The geochemical system that needs to be solved. This gets modified as the solve progresses
ss	Textual information such as (iteration, residual) and swap information is written to this stringstream
tot_iter	The total number of iterations used in the solve
abs_residual	The residual of the algebraic system as this method exits
mole_additions	the increment of mole number of each basis species and kinetic species since the last timestep. This must have size _num_basis + _num_kin. For the basis species, this is the amount of each species being injected into the system over the timestep. For the kinetic species, this is -dt*reaction_rate. Please note: do not decompose the kinetic-species additions into basis components and add them to the first slots of mole_additions! This method does that decomposition automatically. The first _num_basis slots of mole_additions contain kinetic-independent additions, while the last _num_kin slots contain kinetic-rate contributions. This may change during the solution process as kinetic rates depend on molalities.
dmole_additions	dmole_additions(a, b) = d(mole_additions(a))/d(basis_molality(b)) for b < _num_basis and d(mole_additions(a))/d(kinetic_moles(b - _num_basis)) otherwise.

Definition at line 316 of file GeochemicalSolver.C.

Referenced by GeochemistryTimeDependentReactor::execute(), GeochemistrySpatialReactor::execute(), GeochemistryTimeIndependentReactor::initialSetup(), GeochemistryTimeDependentReactor::initialSetup(), GeochemistrySpatialReactor::initialSetup(), GeochemistryTimeDependentReactor::postSolveExchanger(), and TEST().

{
  ss.str("");
  tot_iter = 0;
  unsigned num_swaps = 0;

  // capture the inputs expressed in the current basis
  _input_mole_additions = mole_additions;
  _input_dmole_additions = dmole_additions;

  const ModelGeochemicalDatabase & mgd = egs.getModelGeochemicalDatabase();

  bool still_swapping = true;
  while (still_swapping && num_swaps <= _max_swaps_allowed)
  {
    _num_basis_in_algebraic_system = egs.getNumBasisInAlgebraicSystem();
    _num_in_algebraic_system = egs.getNumInAlgebraicSystem();
    _residual = DenseVector<Real>(_num_in_algebraic_system);
    _jacobian = DenseMatrix<Real>(_num_in_algebraic_system, _num_in_algebraic_system);
    _new_mol = DenseVector<Real>(_num_in_algebraic_system);

    unsigned iter = 0;
    const Real max_is0 =
        std::min(1.0, (iter + 1.0) / (_ramp_max_ionic_strength + 1.0)) * _max_ionic_strength;
    _is.setMaxIonicStrength(max_is0);
    _is.setMaxStoichiometricIonicStrength(max_is0);

    mole_additions = _input_mole_additions;
    dmole_additions = _input_dmole_additions;
    egs.addKineticRates(dt, mole_additions, dmole_additions);

    _abs_residual = computeResidual(egs, _residual, mole_additions);
    bool reducing_initial_molalities = (_abs_residual > _max_initial_residual);
    const unsigned max_tries =
        _num_basis *
        unsigned(
            std::log(_abs_residual / _max_initial_residual) /
            std::log(1.1)); // assume each successful call to reduceInitialResidual reduces residual
                            // by about 1.1, but that the correct basis species has to be chosen
    unsigned tries = 0;
    while (reducing_initial_molalities && ++tries <= max_tries)
      reducing_initial_molalities = reduceInitialResidual(egs, dt, mole_additions, dmole_additions);

    ss << "iter = " << iter << " |R| = " << _abs_residual << std::endl;
    _res0_times_rel = _abs_residual * _rel_tol;
    while ((_abs_residual >= _res0_times_rel && _abs_residual >= _abs_tol && iter < _max_iter) ||
           iter < _ramp_max_ionic_strength)
    {
      iter += 1;
      tot_iter += 1;
      egs.computeJacobian(_residual, _jacobian, mole_additions, dmole_additions);
      solveAndUnderrelax(egs, _jacobian, _new_mol);
      const Real max_is =
          std::min(1.0, (iter + 1.0) / (_ramp_max_ionic_strength + 1.0)) * _max_ionic_strength;
      _is.setMaxIonicStrength(max_is);
      _is.setMaxStoichiometricIonicStrength(max_is);
      egs.setAlgebraicVariables(_new_mol);
      if (egs.alterChargeBalanceSpecies(_swap_threshold))
        ss << "Changed change balance species to "
           << mgd.basis_species_name.at(egs.getChargeBalanceBasisIndex()) << std::endl;
      if (_evaluate_kin_always)
      {
        mole_additions = _input_mole_additions;
        dmole_additions = _input_dmole_additions;
        egs.addKineticRates(dt, mole_additions, dmole_additions);
      }
      _abs_residual = computeResidual(egs, _residual, mole_additions);
      ss << "iter = " << iter << " |R| = " << _abs_residual << std::endl;
    }

    abs_residual = _abs_residual; // record the final residual

    if (iter >= _max_iter)
      ss << std::endl << "Warning: Number of iterations exceeds " << _max_iter << std::endl;

    unsigned swap_out_of_basis = 0;
    unsigned swap_into_basis = 0;
    try
    {
      // to ensure basis molality is correct:
      for (unsigned i = 0; i < _num_basis; ++i)
        egs.addToBulkMoles(i, mole_additions(i));
      // now check for small basis molalities, minerals consumed and minerals precipitated
      still_swapping = swapNeeded(egs, swap_out_of_basis, swap_into_basis, ss);
    }
    catch (const MooseException & e)
    {
      mooseException(e.what());
    }
    if (still_swapping)
    {
      // need to do a swap and re-solve
      try
      {
        // before swapping, remove any basis mole_additions that came from kinetics.  The
        // following loop over i, combined with the above egs.addToBulkMoles(mole_additions), where
        // mole_additions = _input_mole_additions + kinetic_additions, means that the bulk
        // moles in egs will have only been incremented by _input_mole_additions.  Hence,
        // _input_mole_additions can be set to zero in preparation for the next solve in the new
        // basis.
        // NOTE: if _input_mole_additions depend on molalities, this approach introduces error,
        // because the following call to addToBulkMoles and subsequent _input_mole_additions = 0
        // means the mole additions are FIXED
        for (unsigned i = 0; i < _num_basis; ++i)
        {
          egs.addToBulkMoles(i, _input_mole_additions(i) - mole_additions(i));
          _input_mole_additions(i) = 0.0;
          for (unsigned j = 0; j < _num_basis + _num_kin; ++j)
            _input_dmole_additions(i, j) = 0.0;
        }
        egs.performSwap(swap_out_of_basis, swap_into_basis);
        num_swaps += 1;
      }
      catch (const MooseException & e)
      {
        mooseException(e.what());
      }
    }
  }

  if (num_swaps > _max_swaps_allowed)
  {
    mooseException("Maximum number of swaps performed during solve");
  }
  else if (_abs_residual >= _res0_times_rel && _abs_residual >= _abs_tol)
  {
    mooseException("Failed to converge");
  }
  else
  {
    // the basis species additions have been added above with addtoBulkMoles: now add the kinetic
    // additions:
    for (unsigned i = 0; i < _num_basis; ++i)
      mole_additions(i) = 0.0;
    egs.updateOldWithCurrent(mole_additions);
    egs.enforceChargeBalance();
  }
}

swapNeeded()

bool GeochemicalSolver::swapNeeded ( const GeochemicalSystem &  egs, unsigned &  swap_out_of_basis, unsigned &  swap_into_basis, std::stringstream &  ss  ) const

private

Check if a basis swap is needed.

It is needed if:

• the Newton process did not converge and some non-minerals have small molality
• the free number of moles of a basis mineral is negative
• the saturation index of an equilibrium mineral is positive (and it is not in the prevent_precipitation list)

  Parameters

  egs	Geochemical system that we're trying to solve
  swap_out_of_basis	the index of the species in the basis that will be removed from the basis
  swap_into_basis	the index of the equilibrium mineneral that will be added to the basis

  Returns

  true if a swap is needed, false otherwise

Definition at line 107 of file GeochemicalSolver.C.

Referenced by solveSystem().

{
  bool swap_needed = false;

  const ModelGeochemicalDatabase & mgd = egs.getModelGeochemicalDatabase();

  // check if any basis minerals have negative free number of moles
  const std::vector<Real> & basis_molality = egs.getSolventMassAndFreeMolalityAndMineralMoles();
  const std::vector<unsigned> molality_order =
      GeochemistrySortedIndices::sortedIndices(basis_molality, true);

  // if the Newton did not converge then check for small molalities in the non-minerals
  if (_abs_residual > _res0_times_rel && _abs_residual > _abs_tol)
  {
    for (const auto & i : molality_order)
    {
      if (basis_molality[i] >= _swap_threshold)
      {
        // since we're going through basis_molality in ascending order, as soon as we get >=
        // _swap_threshold, we're safe
        break;
      }
      else if (!mgd.basis_species_mineral[i] && egs.getInAlgebraicSystem()[i] &&
               egs.getChargeBalanceBasisIndex() != i)
      {
        // a non-mineral in the algebraic system has super low molality: try to find a legitimate
        // swap
        swap_out_of_basis = i;
        bool legitimate_swap_found = egs.getSwapper().findBestEqmSwap(swap_out_of_basis,
                                                                      mgd,
                                                                      egs.getEquilibriumMolality(),
                                                                      true,
                                                                      false,
                                                                      false,
                                                                      swap_into_basis);
        if (legitimate_swap_found)
        {
          ss << "Basis species " << mgd.basis_species_name[swap_out_of_basis]
             << " has very low molality of " << basis_molality[i] << " compared to "
             << _swap_threshold << ".  Swapping it with equilibrium species "
             << mgd.eqm_species_name[swap_into_basis] << std::endl;
          swap_needed = true;
          break;
        }
        // if no legitimate swap is found, then loop around to the next basis species
      }
    }
  }
  if (swap_needed)
    return swap_needed;

  // now look through the molalities for minerals that are consumed
  for (const auto & i : molality_order)
  {
    if (basis_molality[i] > 0.0)
    {
      // since we're going through basis_molality in ascending order, as soon as we get >0, we're
      // safe
      break;
    }
    else if (mgd.basis_species_mineral[i])
    {
      swap_needed = true;
      swap_out_of_basis = i;
      bool legitimate_swap_found = egs.getSwapper().findBestEqmSwap(swap_out_of_basis,
                                                                    mgd,
                                                                    egs.getEquilibriumMolality(),
                                                                    false,
                                                                    false,
                                                                    false,
                                                                    swap_into_basis);
      if (!legitimate_swap_found)
        mooseException("Cannot find a legitimate swap for mineral ",
                       mgd.basis_species_name[swap_out_of_basis]);
      ss << "Mineral " << mgd.basis_species_name[swap_out_of_basis]
         << " consumed.  Swapping it with equilibrium species "
         << mgd.eqm_species_name[swap_into_basis] << std::endl;
      break;
    }
  }
  if (swap_needed)
    return swap_needed;

  // check maximum saturation index is not positive
  const std::vector<Real> & eqm_SI = egs.getSaturationIndices();
  const std::vector<unsigned> mineral_SI_order =
      GeochemistrySortedIndices::sortedIndices(eqm_SI, false);

  for (const auto & j : mineral_SI_order)
    if (eqm_SI[j] > 0.0)
      if (mgd.eqm_species_mineral[j])
        if (std::find(_prevent_precipitation.begin(),
                      _prevent_precipitation.end(),
                      mgd.eqm_species_name[j]) == _prevent_precipitation.end())
        {
          // mineral has positive saturation index and user is not preventing its precipitation.
          // determine the basis species to swap out
          swap_needed = true;
          swap_into_basis = j;
          bool legitimate_swap_found = false;
          swap_out_of_basis = 0;
          Real best_stoi = 0.0;
          for (unsigned i = 1; i < _num_basis; ++i) // never swap water (i=0)
          {
            if (basis_molality[i] > 0.0 && i != egs.getChargeBalanceBasisIndex() &&
                !mgd.basis_species_gas[i] && !egs.getBasisActivityKnown()[i])
            {
              // don't want to swap out the charge-balance species or any gases of fixed fugacity
              // or any species with fixed activity
              const Real stoi = std::abs(mgd.eqm_stoichiometry(j, i)) / basis_molality[i];
              if (stoi > best_stoi)
              {
                legitimate_swap_found = true;
                best_stoi = stoi;
                swap_out_of_basis = i;
              }
            }
          }
          if (!legitimate_swap_found)
            mooseException(
                "Cannot find a legitimate swap for the supersaturated equilibrium species ",
                mgd.eqm_species_name[j]);
          ss << "Mineral " << mgd.eqm_species_name[j]
             << " supersaturated.  Swapping it with basis species "
             << mgd.basis_species_name[swap_out_of_basis] << std::endl;
          break;
        }
  return swap_needed;
}

Member Data Documentation

_abs_residual

Real GeochemicalSolver::_abs_residual

private

L1 norm of residual.

Definition at line 118 of file GeochemicalSolver.h.

Referenced by reduceInitialResidual(), solveSystem(), and swapNeeded().

_abs_tol

const Real GeochemicalSolver::_abs_tol

private

If the residual of the algebraic system falls below this value, the Newton process has converged.

Definition at line 124 of file GeochemicalSolver.h.

Referenced by solveSystem(), and swapNeeded().

_evaluate_kin_always

bool GeochemicalSolver::_evaluate_kin_always

private

When to compute the kinetic rates: if true then evaluate before every residual calculation, otherwise evaluate only at the start of the solve process (using, eg, the old molality values)

Definition at line 147 of file GeochemicalSolver.h.

Referenced by reduceInitialResidual(), and solveSystem().

_input_dmole_additions

DenseMatrix<Real> GeochemicalSolver::_input_dmole_additions

private

d(mole_additions)/d(species) as specified in the argument of solveSystem

Definition at line 151 of file GeochemicalSolver.h.

Referenced by reduceInitialResidual(), and solveSystem().

_input_mole_additions

DenseVector<Real> GeochemicalSolver::_input_mole_additions

private

the mole_additions for the basis and kinetic species as specified in the argument of solveSystem

Definition at line 149 of file GeochemicalSolver.h.

Referenced by reduceInitialResidual(), and solveSystem().

_is

GeochemistryIonicStrength& GeochemicalSolver::_is

private

The ionic-strength calculator.

Definition at line 106 of file GeochemicalSolver.h.

Referenced by solveSystem().

_jacobian

DenseMatrix<Real> GeochemicalSolver::_jacobian

private

jacobian of the algebraic system

Definition at line 120 of file GeochemicalSolver.h.

Referenced by solveSystem().

_max_initial_residual

Real GeochemicalSolver::_max_initial_residual

private

maximum desired initial residual

Definition at line 132 of file GeochemicalSolver.h.

Referenced by getMaxInitialResidual(), reduceInitialResidual(), setMaxInitialResidual(), and solveSystem().

_max_ionic_strength

const Real GeochemicalSolver::_max_ionic_strength

private

Maximum ionic strength allowed.

Definition at line 140 of file GeochemicalSolver.h.

Referenced by solveSystem().

_max_iter

const unsigned GeochemicalSolver::_max_iter

private

maximum number of iterations allowed during an inner solve

Definition at line 130 of file GeochemicalSolver.h.

Referenced by GeochemicalSolver(), setRampMaxIonicStrength(), and solveSystem().

_max_swaps_allowed

const unsigned GeochemicalSolver::_max_swaps_allowed

private

Maximum number of swaps allowed before the solve aborts.

Definition at line 136 of file GeochemicalSolver.h.

Referenced by solveSystem().

_new_mol

DenseVector<Real> GeochemicalSolver::_new_mol

private

the new molality after finding the solution of _jacobian * neg_change_mol = _residual

Definition at line 122 of file GeochemicalSolver.h.

Referenced by solveSystem().

_num_basis

const unsigned GeochemicalSolver::_num_basis

private

Number of species in the basis.

Definition at line 108 of file GeochemicalSolver.h.

Referenced by solveSystem(), and swapNeeded().

_num_basis_in_algebraic_system

unsigned GeochemicalSolver::_num_basis_in_algebraic_system

private

Number of basis molalities (and potentially solvent water mass) in the algebraic system.

Definition at line 112 of file GeochemicalSolver.h.

Referenced by reduceInitialResidual(), and solveSystem().

_num_in_algebraic_system

unsigned GeochemicalSolver::_num_in_algebraic_system

private

Number of unknowns (molalities and surface potentials) in the algebraic system.

Definition at line 114 of file GeochemicalSolver.h.

Referenced by computeResidual(), solveAndUnderrelax(), and solveSystem().

_num_kin

const unsigned GeochemicalSolver::_num_kin

private

Number of kinetic species.

Definition at line 110 of file GeochemicalSolver.h.

Referenced by solveSystem().

_prevent_precipitation

const std::vector<std::string> GeochemicalSolver::_prevent_precipitation

private

The minerals named in this list can have positive saturation indices and will not precipitate.

Definition at line 138 of file GeochemicalSolver.h.

Referenced by swapNeeded().

_ramp_max_ionic_strength

unsigned GeochemicalSolver::_ramp_max_ionic_strength

private

Number of iterations over which to increase the maximum ionic strength to _max_ionic_strength.

Definition at line 142 of file GeochemicalSolver.h.

Referenced by GeochemicalSolver(), getRampMaxIonicStrength(), setRampMaxIonicStrength(), and solveSystem().

_rel_tol

const Real GeochemicalSolver::_rel_tol

private

If the residual of the algebraic system falls below this value times the initial residual, the Newton process has converged.

Definition at line 126 of file GeochemicalSolver.h.

Referenced by solveSystem().

_res0_times_rel

Real GeochemicalSolver::_res0_times_rel

private

_res0_times_rel = _rel_tol * initial residual

Definition at line 128 of file GeochemicalSolver.h.

Referenced by solveSystem(), and swapNeeded().

_residual

DenseVector<Real> GeochemicalSolver::_residual

private

residual of the algebraic system we wish to solve

Definition at line 116 of file GeochemicalSolver.h.

Referenced by reduceInitialResidual(), solveAndUnderrelax(), and solveSystem().

_swap_threshold

const Real GeochemicalSolver::_swap_threshold

private

If a basis molality < swap_threshold, we attempt to swap it out of the basis.

Definition at line 134 of file GeochemicalSolver.h.

Referenced by solveSystem(), and swapNeeded().

---

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

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