This class holds information about bulk composition, molalities, activities, activity coefficients, etc of the user-defined geochemical system in PertinentGeochemicalSystem. More...

#include <GeochemicalSystem.h>

Public Types
enum class	ConstraintMeaningEnum { MOLES_BULK_WATER , KG_SOLVENT_WATER , MOLES_BULK_SPECIES , FREE_MOLALITY , FREE_MOLES_MINERAL_SPECIES , FUGACITY , ACTIVITY }
	Each basis species has one of the following constraints. More...

enum class	ConstraintUserMeaningEnum { KG_SOLVENT_WATER , BULK_COMPOSITION , BULK_COMPOSITION_WITH_KINETIC , FREE_CONCENTRATION , FREE_MINERAL , ACTIVITY , LOG10ACTIVITY , FUGACITY , LOG10FUGACITY }
	Each basis species must be provided with a constraint, that is chosen by the user from the following enum. More...

Public Member Functions
	GeochemicalSystem (ModelGeochemicalDatabase &mgd, GeochemistryActivityCoefficients &gac, GeochemistryIonicStrength &is, GeochemistrySpeciesSwapper &swapper, const std::vector< std::string > &swap_out_of_basis, const std::vector< std::string > &swap_into_basis, const std::string &charge_balance_species, const std::vector< std::string > &constrained_species, const std::vector< Real > &constraint_value, const MultiMooseEnum &constraint_unit, const MultiMooseEnum &constraint_user_meaning, Real initial_temperature, unsigned iters_to_make_consistent, Real min_initial_molality, const std::vector< std::string > &kin_name, const std::vector< Real > &kin_initial_moles, const MultiMooseEnum &kin_unit)
	Construct the geochemical system, check consistency of the constructor's arguments, and initialize all internal variables (molalities, bulk compositions, equilibrium constants, activities, activity coefficients, surface potentials and kinetic mole numbers) and set up the algebraic system.

	GeochemicalSystem (ModelGeochemicalDatabase &mgd, GeochemistryActivityCoefficients &gac, GeochemistryIonicStrength &is, GeochemistrySpeciesSwapper &swapper, const std::vector< std::string > &swap_out_of_basis, const std::vector< std::string > &swap_into_basis, const std::string &charge_balance_species, const std::vector< std::string > &constrained_species, const std::vector< Real > &constraint_value, const std::vector< GeochemistryUnitConverter::GeochemistryUnit > &constraint_unit, const std::vector< ConstraintUserMeaningEnum > &constraint_user_meaning, Real initial_temperature, unsigned iters_to_make_consistent, Real min_initial_molality, const std::vector< std::string > &kin_name, const std::vector< Real > &kin_initial, const std::vector< GeochemistryUnitConverter::GeochemistryUnit > &kin_unit)

GeochemicalSystem &	operator= (const GeochemicalSystem &src)
	Copy assignment operator.

	GeochemicalSystem (const GeochemicalSystem &src)=default
	Copy constructor.

bool	operator== (const GeochemicalSystem &rhs) const

unsigned	getNumInBasis () const
	returns the number of species in the basis

unsigned	getNumInEquilibrium () const
	returns the number of species in equilibrium with the basis components

unsigned	getNumKinetic () const
	returns the number of kinetic species

void	initialize ()
	initialize all variables, ready for a Newton solve of the algebraic system

unsigned	getChargeBalanceBasisIndex () const
	return the index of the charge-balance species in the basis list

bool	alterChargeBalanceSpecies (Real threshold_molality)
	If the molality of the charge-balance basis species is less than threshold molality, then attempt to change the charge-balance species, preferably to another component with opposite charge and big molality.

Real	getLog10K (unsigned j) const

unsigned	getNumRedox () const

Real	getRedoxLog10K (unsigned red) const

Real	log10RedoxActivityProduct (unsigned red) const

Real	getKineticLog10K (unsigned kin) const

const std::vector< Real > &	getKineticLog10K () const

Real	log10KineticActivityProduct (unsigned kin) const

Real	getKineticMoles (unsigned kin) const

void	setKineticMoles (unsigned kin, Real moles)
	Sets the current AND old mole number for a kinetic species.

const std::vector< Real > &	getKineticMoles () const

void	updateOldWithCurrent (const DenseVector< Real > &mole_additions)
	Usually used at the end of a solve, to provide correct initial conditions to the next time-step, this method:

unsigned	getNumInAlgebraicSystem () const

unsigned	getNumBasisInAlgebraicSystem () const
	return the number of basis species in the algebraic system

unsigned	getNumSurfacePotentials () const
	return the number of surface potentials

const std::vector< bool > &	getInAlgebraicSystem () const

const std::vector< unsigned > &	getBasisIndexOfAlgebraicSystem () const

const std::vector< unsigned > &	getAlgebraicIndexOfBasisSystem () const

std::vector< Real >	getAlgebraicVariableValues () const

std::vector< Real >	getAlgebraicBasisValues () const

DenseVector< Real >	getAlgebraicVariableDenseValues () const

Real	getSolventWaterMass () const
	Returns the mass of solvent water.

const std::vector< Real > &	getBulkMolesOld () const

DenseVector< Real >	getBulkOldInOriginalBasis () const

DenseVector< Real >	getTransportedBulkInOriginalBasis () const

void	computeTransportedBulkFromMolalities (std::vector< Real > &transported_bulk) const
	Computes the value of transported bulk moles for all basis species using the existing molalities.

const std::vector< Real > &	getSolventMassAndFreeMolalityAndMineralMoles () const

const std::vector< bool > &	getBasisActivityKnown () const

Real	getBasisActivity (unsigned i) const

const std::vector< Real > &	getBasisActivity () const

Real	getEquilibriumMolality (unsigned j) const

const std::vector< Real > &	getEquilibriumMolality () const

Real	getBasisActivityCoefficient (unsigned i) const

const std::vector< Real > &	getBasisActivityCoefficient () const

Real	getEquilibriumActivityCoefficient (unsigned j) const

const std::vector< Real > &	getEquilibriumActivityCoefficient () const

Real	getTotalChargeOld () const

Real	getResidualComponent (unsigned algebraic_ind, const DenseVector< Real > &mole_additions) const
	Return the residual of the algebraic system for the given algebraic index.

const ModelGeochemicalDatabase &	getModelGeochemicalDatabase () const

void	setModelGeochemicalDatabase (const ModelGeochemicalDatabase &mgd)
	Copies a ModelGeochemicalDatabase into our _mgd structure.

void	computeJacobian (const DenseVector< Real > &res, DenseMatrix< Real > &jac, const DenseVector< Real > &mole_additions, const DenseMatrix< Real > &dmole_additions) const
	Compute the Jacobian for the algebraic system and put it in jac.

void	setAlgebraicVariables (const DenseVector< Real > &algebraic_var)
	Set the variables in the algebraic system (molalities and potentially the mass of solvent water, and surface potentials, and kinetic mole numbers if any) to algebraic_var.

void	enforceChargeBalance ()
	Enforces charge balance by altering the constraint_value and bulk_moles_old of the charge-balance species.

bool	revertToOriginalChargeBalanceSpecies ()
	Changes the charge-balance species to the original that is specified in the constructor (this might not be possible since the original charge-balance species might have been swapped out)

std::vector< Real >	getSaturationIndices () const

void	performSwap (unsigned swap_out_of_basis, unsigned swap_into_basis)
	Perform the basis swap, and ensure that the resulting system is consistent.

Real	getIonicStrength () const
	Get the ionic strength.

Real	getStoichiometricIonicStrength () const
	Get the stoichiometric ionic strength.

Real	getSurfacePotential (unsigned sp) const

Real	getSurfaceCharge (unsigned sp) const

const std::vector< Real > &	getSorbingSurfaceArea () const

Real	getTemperature () const

void	setTemperature (Real temperature)
	Sets the temperature to the given quantity.

void	setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles (const std::vector< std::string > &names, const std::vector< Real > &values, const std::vector< bool > &constraints_from_molalities)
	Sets:

const std::vector< GeochemicalSystem::ConstraintMeaningEnum > &	getConstraintMeaning () const

void	closeSystem ()
	Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and all FREE_MOLALITY and FREE_MOLES_MINERAL_SPECIES to MOLES_BULK_SPECIES using the current values of _bulk_moles_old.

void	changeConstraintToBulk (unsigned basis_ind)
	Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis index.

void	changeConstraintToBulk (unsigned basis_ind, Real value)
	Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis index.

void	addToBulkMoles (unsigned basis_ind, Real value)
	Add to the MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis species.

void	setConstraintValue (unsigned basis_ind, Real value)
	Set the constraint value for the basis species.

void	computeConsistentConfiguration ()
	Compute a reasonably consistent configuration that can be used as an initial condition in a Newton process to solve the algebraic system.

const std::string &	getOriginalRedoxLHS () const

void	performSwapNoCheck (unsigned swap_out_of_basis, unsigned swap_into_basis)
	Perform the basis swap, and ensure that the resulting system is consistent.

const GeochemistrySpeciesSwapper &	getSwapper () const

void	setMineralRelatedFreeMoles (Real value)
	Set the free mole number of mineral-related species to the value provided.

const std::vector< Real > &	computeAndGetEquilibriumActivity ()
	Computes activity for all equilibrium species (_eqm_activity) and returns a reference to the vector.

Real	getEquilibriumActivity (unsigned eqm_ind) const
	Returns the value of activity for the equilibrium species with index eqm_index.

void	addKineticRates (Real dt, DenseVector< Real > &mole_additions, DenseMatrix< Real > &dmole_additions)
	Computes the kinetic rates and their derivatives based on the current values of molality, mole number, etc, and then, using dt, calculates the appropriate mole_additions and dmole_additions.

Private Member Functions
void	buildTemperatureDependentQuantities (Real temperature)
	Using the provided value of temperature, build _eqm_log10K for each eqm species and redox species.

void	buildAlgebraicInfo (std::vector< bool > &in_algebraic_system, unsigned &num_basis_in_algebraic_system, unsigned &num_in_algebraic_system, std::vector< unsigned > &algebraic_index, std::vector< unsigned > &basis_index) const
	Builds in_algebraic_system, algebraic_index and basis_index, and sets num_basis_in_algebraic_system appropriately.

void	initBulkAndFree (std::vector< Real > &bulk_moles_old, std::vector< Real > &basis_molality) const
	based on _constrained_value and _constraint_meaning, populate nw, bulk_moles_old and basis_molality with reasonable initial conditions that may be used during the Newton solve of the algebraic system

void	buildKnownBasisActivities (std::vector< bool > &basis_activity_known, std::vector< Real > &basis_activity) const
	Fully populates basis_activity_known, which is true if activity has been set by the user, or the fugacity has been set by the user, or the basis species is a mineral.

void	computeRemainingBasisActivities (std::vector< Real > &basis_activity) const
	Compute the activity of water and put in basis_activity[0], and use the _basis_activity_coef and _basis_molality to compute the remaining basis activities (for those that are not minerals, gases or aqueous species provided with an activity)

void	updateBasisMolalityForKnownActivity (std::vector< Real > &basis_molality) const
	For basis aqueous species (not water, minerals or gases) with activity set by the user, set basis_molality = activity / activity_coefficient.

Real	log10ActivityProduct (unsigned eqm_j) const

void	computeEqmMolalities (std::vector< Real > &eqm_molality) const
	compute the equilibrium molalities (not for minerals or gases)

void	enforceChargeBalanceIfSimple (std::vector< Real > &constraint_value, std::vector< Real > &bulk_moles_old) const
	If all non-charged basis species are provided with a bulk number of moles, alter the bulk number of moles (old) of the charge_balance_species so that charges are balanced.

void	computeBulk (std::vector< Real > &bulk_moles_old) const
	Computes the value of bulk_moles_old for all basis species.

void	enforceChargeBalance (std::vector< Real > &constraint_value, std::vector< Real > &bulk_moles_old) const
	Enforces charge balance by altering the constraint_value and bulk_moles_old of the charge-balance species.

void	computeFreeMineralMoles (std::vector< Real > &basis_molality) const
	Compute the free mineral moles (ie, basis_molality for basis species that are minerals), using the bulk_moles_old.

void	checkAndInitialize (const std::vector< std::string > &kin_name, const std::vector< Real > &kin_initial, const std::vector< GeochemistryUnitConverter::GeochemistryUnit > &kin_unit)
	Used during construction: checks for sane inputs and initializes molalities, etc, using the initialize() method.

void	setChargeBalanceSpecies (unsigned new_charge_balance_index)
	Set the charge-balance species to the basis index provided.

Real	surfaceSorptionModifier (unsigned eqm_j) const

Real	surfacePotPrefactor (unsigned sp) const

void	computeSorbingSurfaceArea (std::vector< Real > &sorbing_surface_area) const
	Compute sorbing_surface_area depending on the current molality of the sorbing minerals.

Real	computeBulkFromMolalities (unsigned basis_ind) const

void	alterSystemBecauseBulkChanged ()
	Alter the GeochemicalSystem to reflect changes in bulk composition constraints that occur through, for instance, setConstraintValue or addToBulkMoles or changeConstraintToBulk (the latter because charge neutrality might be easily enforced, changing constraints and hence bulk moles).

Private Attributes
ModelGeochemicalDatabase &	_mgd
	The minimal geochemical database.

const unsigned	_num_basis
	number of basis species

const unsigned	_num_eqm
	number of equilibrium species

const unsigned	_num_redox
	number of redox species

const unsigned	_num_surface_pot
	number of surface potentials

const unsigned	_num_kin
	number of kinetic species

GeochemistrySpeciesSwapper &	_swapper
	swapper that can swap species into and out of the basis

const std::vector< std::string >	_swap_out
	Species to immediately remove from the basis in favor of _swap_in.

const std::vector< std::string >	_swap_in
	Species to immediately add to the basis in favor of _swap_out.

GeochemistryActivityCoefficients &	_gac
	Object to compute the activity coefficients and activity of water.

GeochemistryIonicStrength &	_is
	Object that provides the ionic strengths.

std::string	_charge_balance_species
	The species used to enforce charge balance.

const std::string	_original_charge_balance_species
	The species used to enforce charge balance, as provided in the constructor.

unsigned	_charge_balance_basis_index
	The index in the list of basis species corresponding to the charge-balance species. This gets altered as the simulation progresses, due to swaps, and alterChargeBalanceSpecies.

std::vector< std::string >	_constrained_species
	Names of the species in that have their values fixed to _constraint_value with meanings _constraint_meaning. In the constructor, this is ordered to have the same ordering as the basis species.

std::vector< Real >	_constraint_value
	Numerical values of the constraints on _constraint_species. In the constructor, this is ordered to have the same ordering as the basis species.

std::vector< Real >	_original_constraint_value
	Numerical values of the constraints on _constraint_species. In the constructor, this is ordered to have the same ordering as the basis species. Since values can change due to charge-balance, this holds the original values set by the user.

std::vector< GeochemistryUnitConverter::GeochemistryUnit >	_constraint_unit
	Units of the _constraint_value when the GeochemicalSystem is constructed. This is used during the constructor to convert the constraint_values into mole units.

std::vector< ConstraintUserMeaningEnum >	_constraint_user_meaning
	The user-defined meaning of the values in _constraint_value. In the constructor, this is ordered to have the same ordering as the basis species. During the process of unit conversion, from the user-supplied _constraint_unit, to mole-based units, _constraint_user_meaning is used to populate _constraint_meaning, and henceforth usually only _constraint_meaning is used in the code.

std::vector< ConstraintMeaningEnum >	_constraint_meaning
	The meaning of the values in _constraint_value. In the constructor, this is ordered to have the same ordering as the basis species.

std::vector< Real >	_eqm_log10K
	equilibrium constant of the equilibrium species

std::vector< Real >	_redox_log10K
	equilibrium constant of the redox species

std::vector< Real >	_kin_log10K
	equilibrium constant of the kinetic species

unsigned	_num_basis_in_algebraic_system
	number of unknown molalities in the algebraic system. Note: surface potentials (if any) are extra quantities in the algebraic system

unsigned	_num_in_algebraic_system
	number of unknowns in the algebraic system (includes molalities and surface potentials)

std::vector< bool >	_in_algebraic_system
	if _in_algebraic_system[i] == true then we need to solve for the i^th basis-species molality

std::vector< unsigned >	_algebraic_index
	_algebraic_index[i] = index in the algebraic system of the basis species i. _basis_index[_algebraic_index[i]] = i

std::vector< unsigned >	_basis_index
	_basis_index[i] = index in the basis of the algebraic system species i, for i<num_basis_in_algebraic_system. _basis_index[_algebraic_index[i]] == i

std::vector< Real >	_bulk_moles_old
	Number of bulk moles of basis species (this includes contributions from kinetic species, in contrast to Bethke)

std::vector< Real >	_basis_molality
	IMPORTANT: this is.

std::vector< bool >	_basis_activity_known
	whether basis_activity[i] is fixed by the user

std::vector< Real >	_basis_activity
	values of activity (for water, minerals and aqueous basis species) or fugacity (for gases)

std::vector< Real >	_eqm_molality
	molality of equilibrium species

std::vector< Real >	_basis_activity_coef
	basis activity coefficients

std::vector< Real >	_eqm_activity_coef
	equilibrium activity coefficients

std::vector< Real >	_eqm_activity
	equilibrium activities.

std::vector< Real >	_surface_pot_expr
	surface potential expressions.

std::vector< Real >	_sorbing_surface_area
	surface areas of the sorbing minerals

std::vector< Real >	_kin_moles
	mole number of kinetic species

std::vector< Real >	_kin_moles_old
	old mole number of kinetic species

unsigned	_iters_to_make_consistent
	Iterations to make the initial configuration consistent.

Real	_temperature
	The temperature in degC.

Real	_min_initial_molality
	Minimum molality ever used in an initial guess.

std::string	_original_redox_lhs
	The left-hand-side of the redox equations in _mgd after initial swaps.

Detailed Description

This class holds information about bulk composition, molalities, activities, activity coefficients, etc of the user-defined geochemical system in PertinentGeochemicalSystem.

Hence, it extends the generic information in PertinentGeochemicalSystem to a system that is specific to one spatio-temporal location. It offers methods to manipulate these molalities, bulk compositions, etc, as well as performing basis swaps.

Related to this, this class also builds the so-called "algebraic system" which is the nonlinear system of algebraic equations for the unknown basis molalities, surface potentials and mole number of kinetic species, and computes the residual vector and jacobian of this algebraic system. This class initialises the algebraic variables (the unknown molalities, surface potentials and kinetic moles) to reaponable values, but it does not solve the algebraic system: instead, it has a "setter" method, setAlgebraicValues that can be used to set the unknowns, hence another class can use whatever method it desires to solve the system.

WARNING: because GeochemicalSystem performs swaps, it will change the ModelGeochemicalDatabase object. WARNING: because GeochemicalSystem sets internal parameters in the activity-coefficient object, it will change the GeochemistryActivityCoefficient object.

Definition at line 39 of file GeochemicalSystem.h.

Member Enumeration Documentation

◆ ConstraintMeaningEnum

enum class GeochemicalSystem::ConstraintMeaningEnum

strong

Each basis species has one of the following constraints.

During the process of units conversion (from user-prescribed units to mole-based units) each ConstraintMeaningUserEnum supplied by the user is translated to the appropriate ConstraintMeaningEnum

Enumerator
MOLES_BULK_WATER
KG_SOLVENT_WATER
MOLES_BULK_SPECIES
FREE_MOLALITY
FREE_MOLES_MINERAL_SPECIES
FUGACITY
ACTIVITY

Definition at line 47 of file GeochemicalSystem.h.

  {
    MOLES_BULK_WATER,
    KG_SOLVENT_WATER,
    MOLES_BULK_SPECIES,
    FREE_MOLALITY,
    FREE_MOLES_MINERAL_SPECIES,
    FUGACITY,
    ACTIVITY
  };

◆ ConstraintUserMeaningEnum

enum class GeochemicalSystem::ConstraintUserMeaningEnum

strong

Each basis species must be provided with a constraint, that is chosen by the user from the following enum.

Enumerator
KG_SOLVENT_WATER
BULK_COMPOSITION
BULK_COMPOSITION_WITH_KINETIC
FREE_CONCENTRATION
FREE_MINERAL
ACTIVITY
LOG10ACTIVITY
FUGACITY
LOG10FUGACITY

Definition at line 59 of file GeochemicalSystem.h.

  {
    KG_SOLVENT_WATER,
    BULK_COMPOSITION,
    BULK_COMPOSITION_WITH_KINETIC,
    FREE_CONCENTRATION,
    FREE_MINERAL,
    ACTIVITY,
    LOG10ACTIVITY,
    FUGACITY,
    LOG10FUGACITY
  };

Constructor & Destructor Documentation

◆ GeochemicalSystem() [1/3]

GeochemicalSystem::GeochemicalSystem	(	ModelGeochemicalDatabase &	mgd,
		GeochemistryActivityCoefficients &	gac,
		GeochemistryIonicStrength &	is,
		GeochemistrySpeciesSwapper &	swapper,
		const std::vector< std::string > &	swap_out_of_basis,
		const std::vector< std::string > &	swap_into_basis,
		const std::string &	charge_balance_species,
		const std::vector< std::string > &	constrained_species,
		const std::vector< Real > &	constraint_value,
		const MultiMooseEnum &	constraint_unit,
		const MultiMooseEnum &	constraint_user_meaning,
		Real	initial_temperature,
		unsigned	iters_to_make_consistent,
		Real	min_initial_molality,
		const std::vector< std::string > &	kin_name,
		const std::vector< Real > &	kin_initial_moles,
		const MultiMooseEnum &	kin_unit
	)

Construct the geochemical system, check consistency of the constructor's arguments, and initialize all internal variables (molalities, bulk compositions, equilibrium constants, activities, activity coefficients, surface potentials and kinetic mole numbers) and set up the algebraic system.

Parameters

mgd	the model's geochemical database, which is a model-specific version of the full geochemical database. WARNING: because GeochemicalSystem performs swaps, it will change mgd
gac	the object that computes activity coefficients. WARNING: because GeochemicalSystem sets internal parameters in gac, it will change it.
is	object to compute ionic strengths
swapper	object to perform swaps on mgd
swap_out_of_basis	A list of basis species that should be swapped out of the basis and replaced with swap_into_basis, prior to any other computations
swap_into_basis	A list of equilibrium species that should be swapped into the basis in place of swap_out_of_basis
charge_balance_species	The charge-balance species, which must be a basis species after the swaps are performed
constrained_species	A list of the basis species after the swaps have been performed
constraint_values	Numerical values of the constraints placed on the basis species. Each basis species must have exactly one constraint
constraint_unit	Units of the constraint_values. Each constraint_value must have exactly one constraint_unit. This is only used during construction, to convert the constraint_values to mole units
constraint_user_meaning	A list the provides physical meaning to the constraint_values
initial_temperature	Initial temperature
iters_to_make_consistent	The initial equilibrium molalities depend on the activity coefficients, which depend on the basis and equilibrium molalities. This circular dependence means it is usually impossible to define an exactly consistent initial configuration for molalities. An iterative process is used to approach the consistent initial configuration, using iters_to_make_consistent iterations. Usually iters_to_make_consistent=0 is reasonable, because there are so many approximations used in the solution process that a slightly inconsistent initial configuration is fine
min_initial_molality	Minimum value of equilibrium molality used in the initial condition
kin_name	Names of the kinetic species that are provided with initial conditions in kin_initial_moles. All kinetic species must be provided with an initial condition
kin_initial	Values of the initial mole number or mass or volume (depending on kin_unit) for the kinetic species
kin_unit	The units of the numbers given in kin_initial

Definition at line 15 of file GeochemicalSystem.C.

  : _mgd(mgd),
    _num_basis(mgd.basis_species_index.size()),
    _num_eqm(mgd.eqm_species_index.size()),
    _num_redox(mgd.redox_stoichiometry.m()),
    _num_surface_pot(mgd.surface_sorption_name.size()),
    _num_kin(mgd.kin_species_index.size()),
    _swapper(swapper),
    _swap_out(swap_out_of_basis),
    _swap_in(swap_into_basis),
    _gac(gac),
    _is(is),
    _charge_balance_species(charge_balance_species),
    _original_charge_balance_species(charge_balance_species),
    _charge_balance_basis_index(0),
    _constrained_species(constrained_species),
    _constraint_value(constraint_value),
    _original_constraint_value(constraint_value),
    _constraint_unit(constraint_unit.size()),
    _constraint_user_meaning(constraint_user_meaning.size()),
    _constraint_meaning(constraint_user_meaning.size()),
    _eqm_log10K(_num_eqm),
    _redox_log10K(_num_redox),
    _kin_log10K(_num_kin),
    _num_basis_in_algebraic_system(0),
    _num_in_algebraic_system(0),
    _in_algebraic_system(_num_basis),
    _algebraic_index(_num_basis),
    _basis_index(_num_basis),
    _bulk_moles_old(_num_basis),
    _basis_molality(_num_basis),
    _basis_activity_known(_num_basis),
    _basis_activity(_num_basis),
    _eqm_molality(_num_eqm),
    _basis_activity_coef(_num_basis),
    _eqm_activity_coef(_num_eqm),
    _eqm_activity(_num_eqm),
    _surface_pot_expr(_num_surface_pot),
    _sorbing_surface_area(_num_surface_pot),
    _kin_moles(_num_kin),
    _kin_moles_old(_num_kin),
    _iters_to_make_consistent(iters_to_make_consistent),
    _temperature(initial_temperature),
    _min_initial_molality(min_initial_molality),
    _original_redox_lhs()
{
  for (unsigned i = 0; i < constraint_user_meaning.size(); ++i)
    _constraint_user_meaning[i] =
        static_cast<ConstraintUserMeaningEnum>(constraint_user_meaning.get(i));
  for (unsigned i = 0; i < constraint_unit.size(); ++i)
    _constraint_unit[i] =
        static_cast<GeochemistryUnitConverter::GeochemistryUnit>(constraint_unit.get(i));
  const unsigned ku_size = kin_unit.size();
  std::vector<GeochemistryUnitConverter::GeochemistryUnit> k_unit(ku_size);
  for (unsigned i = 0; i < ku_size; ++i)
    k_unit[i] = static_cast<GeochemistryUnitConverter::GeochemistryUnit>(kin_unit.get(i));
  checkAndInitialize(kin_name, kin_initial, k_unit);
}

◆ GeochemicalSystem() [2/3]

GeochemicalSystem::GeochemicalSystem	(	ModelGeochemicalDatabase &	mgd,
		GeochemistryActivityCoefficients &	gac,
		GeochemistryIonicStrength &	is,
		GeochemistrySpeciesSwapper &	swapper,
		const std::vector< std::string > &	swap_out_of_basis,
		const std::vector< std::string > &	swap_into_basis,
		const std::string &	charge_balance_species,
		const std::vector< std::string > &	constrained_species,
		const std::vector< Real > &	constraint_value,
		const std::vector< GeochemistryUnitConverter::GeochemistryUnit > &	constraint_unit,
		const std::vector< ConstraintUserMeaningEnum > &	constraint_user_meaning,
		Real	initial_temperature,
		unsigned	iters_to_make_consistent,
		Real	min_initial_molality,
		const std::vector< std::string > &	kin_name,
		const std::vector< Real > &	kin_initial,
		const std::vector< GeochemistryUnitConverter::GeochemistryUnit > &	kin_unit
	)

Definition at line 90 of file GeochemicalSystem.C.

  : _mgd(mgd),
    _num_basis(mgd.basis_species_index.size()),
    _num_eqm(mgd.eqm_species_index.size()),
    _num_redox(mgd.redox_stoichiometry.m()),
    _num_surface_pot(mgd.surface_sorption_name.size()),
    _num_kin(mgd.kin_species_index.size()),
    _swapper(swapper),
    _swap_out(swap_out_of_basis),
    _swap_in(swap_into_basis),
    _gac(gac),
    _is(is),
    _charge_balance_species(charge_balance_species),
    _original_charge_balance_species(charge_balance_species),
    _charge_balance_basis_index(0),
    _constrained_species(constrained_species),
    _constraint_value(constraint_value),
    _original_constraint_value(constraint_value),
    _constraint_unit(constraint_unit),
    _constraint_user_meaning(constraint_user_meaning),
    _constraint_meaning(constraint_user_meaning.size()),
    _eqm_log10K(_num_eqm),
    _redox_log10K(_num_redox),
    _kin_log10K(_num_kin),
    _num_basis_in_algebraic_system(0),
    _num_in_algebraic_system(0),
    _in_algebraic_system(_num_basis),
    _algebraic_index(_num_basis),
    _basis_index(_num_basis),
    _bulk_moles_old(_num_basis),
    _basis_molality(_num_basis),
    _basis_activity_known(_num_basis),
    _basis_activity(_num_basis),
    _eqm_molality(_num_eqm),
    _basis_activity_coef(_num_basis),
    _eqm_activity_coef(_num_eqm),
    _eqm_activity(_num_eqm),
    _surface_pot_expr(_num_surface_pot),
    _sorbing_surface_area(_num_surface_pot),
    _kin_moles(_num_kin),
    _kin_moles_old(_num_kin),
    _iters_to_make_consistent(iters_to_make_consistent),
    _temperature(initial_temperature),
    _min_initial_molality(min_initial_molality),
    _original_redox_lhs()
{
  checkAndInitialize(kin_name, kin_initial, kin_unit);
}

◆ GeochemicalSystem() [3/3]

GeochemicalSystem::GeochemicalSystem ( const GeochemicalSystem & src )

default

Copy constructor.

Member Function Documentation

◆ addKineticRates()

void GeochemicalSystem::addKineticRates	(	Real	dt,
		DenseVector< Real > &	mole_additions,
		DenseMatrix< Real > &	dmole_additions
	)

Computes the kinetic rates and their derivatives based on the current values of molality, mole number, etc, and then, using dt, calculates the appropriate mole_additions and dmole_additions.

NOTE: this adds the kinetic contributes to mole_additions and dmole_additions. This method is not const because it modifies _eqm_activity for equilibrium gases and eqm species H+ and OH- (if there are any).

Parameters

dt	time-step size
mole_additions	The kinetic rates multiplied by dt get placed in the last num_kin slots
dmole_additions	d(mole_additions)/d(molality or kinetic_moles)

Definition at line 2383 of file GeochemicalSystem.C.

{
  if (_num_kin == 0)
    return;
 
  // check sizes
  const unsigned tot = mole_additions.size();
  if (!(tot == _num_kin + _num_basis && dmole_additions.m() == tot && dmole_additions.n() == tot))
    mooseError("addKineticRates: incorrectly sized additions: ",
               tot,
               " ",
               dmole_additions.m(),
               " ",
               dmole_additions.n());
 
  // construct eqm_activity for species that we need
  for (unsigned j = 0; j < _num_eqm; ++j)
    if (_mgd.eqm_species_gas[j] || _mgd.eqm_species_name[j] == "H+" ||
        _mgd.eqm_species_name[j] == "OH-")
      _eqm_activity[j] = getEquilibriumActivity(j);
 
  // calculate the rates and put into appropriate slots
  Real rate;
  Real drate_dkin;
  std::vector<Real> drate_dmol(_num_basis);
  for (const auto & krd : _mgd.kin_rate)
  {
    const unsigned kin = krd.kinetic_species_index;
    GeochemistryKineticRateCalculator::calculateRate(krd.promoting_indices,
                                                     krd.promoting_monod_indices,
                                                     krd.promoting_half_saturation,
                                                     krd.description,
                                                     _mgd.basis_species_name,
                                                     _mgd.basis_species_gas,
                                                     _basis_molality,
                                                     _basis_activity,
                                                     _basis_activity_known,
                                                     _mgd.eqm_species_name,
                                                     _mgd.eqm_species_gas,
                                                     _eqm_molality,
                                                     _eqm_activity,
                                                     _mgd.eqm_stoichiometry,
                                                     _kin_moles[kin],
                                                     _mgd.kin_species_molecular_weight[kin],
                                                     _kin_log10K[kin],
                                                     log10KineticActivityProduct(kin),
                                                     _mgd.kin_stoichiometry,
                                                     kin,
                                                     _temperature,
                                                     rate,
                                                     drate_dkin,
                                                     drate_dmol);
 
    /* the following block of code may be confusing at first sight.
     * (1) The usual case is that the kinetic reaction is written
     * kinetic -> basis_species, and direction = BOTH, and kinetic_bio_efficiency = -1
     * In this case, the following does mole_additions(kinetic_species) += -1 * rate * dt
     * If rate > 0 then the kinetic species decreases in mass.
     * The solver will then detect that the kinetic_species has changed, and adjust basis_species
     * accordingly, viz it will add rate * dt of basis species to the aqueous solution.
     * (2) The common biologically-catalysed case is the reaction is written
     * reactants -> products, catalysed by microbe.  In the database file this is written as
     * microbe -> -reactants + products
     * The input file will set direction = BOTH, and kinetic_bio_efficiency != -1
     * Suppose that rate > 0, corresponding to reactants being converted to products.
     * With kinetic_bio_efficiency = -1) this would correspond to microbe decreasing in mass,
     * however with kinetic_bio_efficiency > 0, this is not true.
     * The following block of code sets
     * mole_additions(microbe) = kinetic_bio_efficiency * rate * dt > 0
     * that is, the microbe mass is increasing.
     * However, this means the solver will decrease products and increase reactants,
     * by kinetic_bio_efficiency * rate * dt because it sees the reaction microbe -> -reactants +
     * products in the database file.
     * This is clearly wrong, because we actually want the end result to be that the products have
     * increased by rate * dt, and the reactants to have decreased by rate * dt.
     * So, to counter the solver, we add (-1 - kinetic_bio_efficiency) * rate * dt of reactants and
     * add (1 + kinetic_bio_efficienty) * rate * dt of products.
     * (3) The other situation is the biological death, where direction = DEATH.
     * In this case the database file contains microbe -> -reactants + products
     * However, independent of rate, we don't want to change the molality of reactants or products.
     * Following the logic of (2), to counter the solver, we add -kinetic_bio_efficiency * rate * dt
     * of reactants and kinetic_bio_efficiency * rate * dt of products
     */
    const unsigned ind = kin + _num_basis;
    mole_additions(ind) += krd.description.kinetic_bio_efficiency * rate * dt;
    dmole_additions(ind, ind) += krd.description.kinetic_bio_efficiency * drate_dkin * dt;
    const Real extra_additions = (krd.description.direction == DirectionChoiceEnum::DEATH)
                                     ? krd.description.kinetic_bio_efficiency
                                     : krd.description.kinetic_bio_efficiency + 1.0;
    for (unsigned i = 0; i < _num_basis; ++i)
    {
      dmole_additions(ind, i) += krd.description.kinetic_bio_efficiency * drate_dmol[i] * dt;
      const Real stoi_fac = _mgd.kin_stoichiometry(kin, i) * extra_additions * dt;
      mole_additions(i) += stoi_fac * rate;
      dmole_additions(i, ind) += stoi_fac * drate_dkin;
      for (unsigned j = 0; j < _num_basis; ++j)
        dmole_additions(i, j) += stoi_fac * drate_dmol[j];
    }
 
    const Real eff = krd.description.progeny_efficiency;
    if (eff != 0.0)
    {
      const unsigned bio_i = krd.progeny_index;
      if (bio_i < _num_basis)
      {
        mole_additions(bio_i) += eff * rate * dt;
        dmole_additions(bio_i, ind) += eff * drate_dkin * dt;
        for (unsigned i = 0; i < _num_basis; ++i)
          dmole_additions(bio_i, i) += eff * drate_dmol[i] * dt;
      }
      else
      {
        for (unsigned i = 0; i < _num_basis; ++i)
        {
          mole_additions(i) += _mgd.eqm_stoichiometry(bio_i - _num_basis, i) * eff * rate * dt;
          dmole_additions(i, ind) +=
              _mgd.eqm_stoichiometry(bio_i - _num_basis, i) * eff * drate_dkin * dt;
          for (unsigned j = 0; j < _num_basis; ++j)
            dmole_additions(i, j) +=
                _mgd.eqm_stoichiometry(bio_i - _num_basis, i) * eff * drate_dmol[j] * dt;
        }
      }
    }
  }
}

Referenced by GeochemicalSolver::reduceInitialResidual(), GeochemicalSolver::solveSystem(), TEST_F(), TEST_F(), TEST_F(), and TEST_F().

◆ addToBulkMoles()

void GeochemicalSystem::addToBulkMoles	(	unsigned	basis_ind,
		Real	value
	)

Add to the MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis species.

Note that if the constraint on the basis species is not MOLES_BULK_SPECIES (or MOLES_BULK_WATER) then this will have no impact. Note that charge-balance is performed if all constraints are BULK-type.

Parameters

basis_ind	the index of the basis species
value	the value to add to the bulk number of moles

Definition at line 2182 of file GeochemicalSystem.C.

{
  if (basis_ind >= _num_basis)
    mooseError("Cannot addToBulkMoles for species ",
               basis_ind,
               " because there are only ",
               _num_basis,
               " basis species");
  if (!(_constraint_meaning[basis_ind] ==
            GeochemicalSystem::ConstraintMeaningEnum::MOLES_BULK_WATER ||
        _constraint_meaning[basis_ind] ==
            GeochemicalSystem::ConstraintMeaningEnum::MOLES_BULK_SPECIES))
    return;
  setConstraintValue(basis_ind, _constraint_value[basis_ind] + value);
}

Referenced by GeochemistryTimeDependentReactor::postSolveFlowThrough(), GeochemistryTimeDependentReactor::preSolveDump(), GeochemicalSolver::solveSystem(), TEST_F(), TEST_F(), and updateOldWithCurrent().

◆ alterChargeBalanceSpecies()

bool GeochemicalSystem::alterChargeBalanceSpecies ( Real threshold_molality )

If the molality of the charge-balance basis species is less than threshold molality, then attempt to change the charge-balance species, preferably to another component with opposite charge and big molality.

Parameters

threshold_molality if the charge-balance basis species has molality greater than this, this routine does nothing

Returns: true if the charge-balance species was changed

Definition at line 1720 of file GeochemicalSystem.C.

{
  if (_basis_molality[_charge_balance_basis_index] > threshold_molality)
    return false;
  unsigned best_species_opp_sign = 0;
  unsigned best_species_same_sign = 0;
  Real best_molality_opp_sign = 0.0;
  Real best_molality_same_sign = 0.0;
  for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
  {
    if (basis_i == _charge_balance_basis_index)
      continue;
    else if (_constraint_meaning[basis_i] != ConstraintMeaningEnum::MOLES_BULK_SPECIES)
      continue;
    else if (_mgd.basis_species_charge[basis_i] == 0.0)
      continue;
    else if (_basis_molality[basis_i] <= threshold_molality)
      continue;
    // we know basis_i is a charged species with set bulk moles and high molality
    if (_mgd.basis_species_charge[basis_i] *
            _mgd.basis_species_charge[_charge_balance_basis_index] <
        0.0)
    {
      // charge of opposite sign
      if (_basis_molality[basis_i] > best_molality_opp_sign)
      {
        best_species_opp_sign = basis_i;
        best_molality_opp_sign = _basis_molality[basis_i];
      }
    }
    else
    {
      // charge of same sign
      if (_basis_molality[basis_i] > best_molality_same_sign)
      {
        best_species_same_sign = basis_i;
        best_molality_same_sign = _basis_molality[basis_i];
      }
    }
  }
  if (best_species_opp_sign != 0)
  {
    // this is preferred over the same-sign version
    setChargeBalanceSpecies(best_species_opp_sign);
    return true;
  }
  else if (best_species_same_sign != 0)
  {
    // this is not preferred, but is better than no charge-balance species change
    setChargeBalanceSpecies(best_species_same_sign);
    return true;
  }
  return false;
}

Referenced by GeochemicalSolver::solveSystem().

◆ alterSystemBecauseBulkChanged()

void GeochemicalSystem::alterSystemBecauseBulkChanged ( )

private

Alter the GeochemicalSystem to reflect changes in bulk composition constraints that occur through, for instance, setConstraintValue or addToBulkMoles or changeConstraintToBulk (the latter because charge neutrality might be easily enforced, changing constraints and hence bulk moles).

This method enforces charge neutrality if simple (ie, all constraints are BULK-type), then sets _bulk_moles_old, and free mineral moles appropriately

Definition at line 2246 of file GeochemicalSystem.C.

{
  // Altering the constraints on bulk number of moles impacts various other things.
  // Here, we follow the initialize() and computeConsistentConfiguration() methods, picking out
  // what might have changed.
  // Because the constraint meanings have changed, enforcing charge-neutrality might be easy
  enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
  // Because the constraint values have changed, either through charge-neutrality or directly
  // through changing the constraint values, the bulk moles must be updated:
  for (unsigned i = 0; i < _num_basis; ++i)
    if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES ||
        _constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_WATER)
      _bulk_moles_old[i] = _constraint_value[i];
  // Because the bulk mineral moles might have changed, the free mineral moles might have
  // changed
  computeFreeMineralMoles(_basis_molality); // free mineral moles might be changed
}

Referenced by setConstraintValue().

◆ buildAlgebraicInfo()

void GeochemicalSystem::buildAlgebraicInfo	(	std::vector< bool > &	in_algebraic_system,
		unsigned &	num_basis_in_algebraic_system,
		unsigned &	num_in_algebraic_system,
		std::vector< unsigned > &	algebraic_index,
		std::vector< unsigned > &	basis_index
	)		const

private

Builds in_algebraic_system, algebraic_index and basis_index, and sets num_basis_in_algebraic_system appropriately.

Definition at line 609 of file GeochemicalSystem.C.

{
  // build in_algebraic_system
  for (const auto & name_index : _mgd.basis_species_index)
  {
    const std::string name = name_index.first;
    const unsigned basis_ind = name_index.second;
    const ConstraintMeaningEnum meaning = _constraint_meaning[basis_ind];
    if (name == "H2O")
      in_algebraic_system[basis_ind] = (meaning == ConstraintMeaningEnum::MOLES_BULK_WATER);
    else if (_mgd.basis_species_gas[basis_ind])
      in_algebraic_system[basis_ind] = false;
    else if (_mgd.basis_species_mineral[basis_ind])
      in_algebraic_system[basis_ind] = false;
    else
      in_algebraic_system[basis_ind] = (meaning == ConstraintMeaningEnum::MOLES_BULK_SPECIES);
  }
 
  // build algebraic_index and basis_index
  num_basis_in_algebraic_system = 0;
  algebraic_index.resize(_num_basis, 0);
  basis_index.resize(_num_basis, 0);
  for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
    if (in_algebraic_system[basis_ind])
    {
      algebraic_index[basis_ind] = _num_basis_in_algebraic_system;
      basis_index[_num_basis_in_algebraic_system] = basis_ind;
      num_basis_in_algebraic_system += 1;
    }
 
  num_in_algebraic_system = num_basis_in_algebraic_system + _num_surface_pot + _num_kin;
}

Referenced by changeConstraintToBulk(), initialize(), and performSwapNoCheck().

◆ buildKnownBasisActivities()

void GeochemicalSystem::buildKnownBasisActivities	(	std::vector< bool > &	basis_activity_known,
		std::vector< Real > &	basis_activity
	)		const

private

Fully populates basis_activity_known, which is true if activity has been set by the user, or the fugacity has been set by the user, or the basis species is a mineral.

Populates only those slots in basis_activity for which basis_activity_known == true

Definition at line 837 of file GeochemicalSystem.C.

{
  basis_activity_known = std::vector<bool>(_num_basis, false);
  basis_activity.resize(_num_basis);
 
  // all aqueous species with provided activity, and all gases with fugacity have known activity
  for (unsigned i = 0; i < _num_basis; ++i)
  {
    const ConstraintMeaningEnum meaning = _constraint_meaning[i];
    if (meaning == ConstraintMeaningEnum::ACTIVITY || meaning == ConstraintMeaningEnum::FUGACITY)
    {
      basis_activity_known[i] = true;
      basis_activity[i] = _constraint_value[i];
    }
  }
 
  // all minerals have activity = 1.0
  for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
    if (_mgd.basis_species_mineral[basis_ind])
    {
      basis_activity_known[basis_ind] = true;
      basis_activity[basis_ind] = 1.0;
    }
}

Referenced by changeConstraintToBulk(), initialize(), performSwapNoCheck(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ buildTemperatureDependentQuantities()

void GeochemicalSystem::buildTemperatureDependentQuantities ( Real temperature )

private

Using the provided value of temperature, build _eqm_log10K for each eqm species and redox species.

Definition at line 579 of file GeochemicalSystem.C.

{
  const std::vector<Real> temps = _mgd.original_database->getTemperatures();
  const unsigned numT = temps.size();
  const std::string model_type = _mgd.original_database->getLogKModel();
 
  for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
  {
    EquilibriumConstantInterpolator interp(
        temps, _mgd.eqm_log10K.sub_matrix(eqm_j, 1, 0, numT).get_values(), model_type);
    interp.generate();
    _eqm_log10K[eqm_j] = interp.sample(temperature);
  }
  for (unsigned red = 0; red < _num_redox; ++red)
  {
    EquilibriumConstantInterpolator interp(
        temps, _mgd.redox_log10K.sub_matrix(red, 1, 0, numT).get_values(), model_type);
    interp.generate();
    _redox_log10K[red] = interp.sample(temperature);
  }
  for (unsigned kin = 0; kin < _num_kin; ++kin)
  {
    EquilibriumConstantInterpolator interp(
        temps, _mgd.kin_log10K.sub_matrix(kin, 1, 0, numT).get_values(), model_type);
    interp.generate();
    _kin_log10K[kin] = interp.sample(temperature);
  }
}

Referenced by initialize(), performSwapNoCheck(), and setTemperature().

◆ changeConstraintToBulk() [1/2]

void GeochemicalSystem::changeConstraintToBulk ( unsigned basis_ind )

Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis index.

The constraint value (the bulk number of moles of the species) is computed from the current molality values. Note that if basis_ind corresponds to a gas, then changing the constraint involves performing a basis swap with a non-gaseous equilibrium species (eg O2(g) is removed from the basis and O2(aq) enters in its place). If, after performing the change to MOLES_BULK_SPECIES, all constraints are MOLES_BULK_WATER or MOLES_BULK_SPECIES then charge-balance is performed by altering the bulk_moles_old for the charge-balance species.

Parameters

basis_ind the index of the basis species

Definition at line 2110 of file GeochemicalSystem.C.

{
  if (basis_ind >= _num_basis)
    mooseError("Cannot changeConstraintToBulk for species ",
               basis_ind,
               " because there are only ",
               _num_basis,
               " basis species");
  if (_mgd.basis_species_gas[basis_ind])
  {
    // this is a special case where we have to swap out the gas in favour of an equilibrium
    // species
    unsigned swap_into_basis = 0;
    bool legitimate_swap_found = false;
    Real best_stoi = 0.0;
    for (unsigned j = 0; j < _num_eqm; ++j)
    {
      if (_mgd.eqm_species_gas[j] || _mgd.eqm_stoichiometry(j, basis_ind) == 0.0 ||
          _mgd.surface_sorption_related[j])
        continue;
      const Real stoi = std::abs(_mgd.eqm_stoichiometry(j, basis_ind));
      if (stoi > best_stoi)
      {
        best_stoi = stoi;
        swap_into_basis = j;
        legitimate_swap_found = true;
      }
    }
    if (legitimate_swap_found)
      performSwapNoCheck(basis_ind, swap_into_basis);
    else
      mooseError("Attempting to change constraint of gas ",
                 _mgd.basis_species_name[basis_ind],
                 " to MOLES_BULK_SPECIES, which involves a search for a suitable non-gas species "
                 "to swap with.  No such species was found");
  }
  else
    changeConstraintToBulk(basis_ind, computeBulkFromMolalities(basis_ind));
}

Referenced by changeConstraintToBulk(), closeSystem(), GeochemistryTimeDependentReactor::execute(), and TEST_F().

◆ changeConstraintToBulk() [2/2]

void GeochemicalSystem::changeConstraintToBulk	(	unsigned	basis_ind,
		Real	value
	)

Changes the constraint to MOLES_BULK_SPECIES (or MOLES_BULK_WATER if basis_ind = 0) for the basis index.

The constraint value (the bulk number of moles of the species) is set to value. It is an error to call this function when basis_ind corresponds to a gas, because changing the constraint involves performing a basis swap with a non-gaseous equilibrium species (eg O2(g) is removed from the basis and O2(aq) enters in its place), so the bulk number of moles will be computed by the swapper: use changeConstraintToBulk(basis_ind) instead. If, after performing the change to MOLES_BULK_SPECIES, all constraints are MOLES_BULK_WATER or MOLES_BULK_SPECIES then charge-balance is performed by altering the bulk_moles_old for the charge-balance species.

Parameters

basis_ind	the index of the basis species
value	the value to set the bulk number of moles to

Definition at line 2151 of file GeochemicalSystem.C.

{
  if (basis_ind >= _num_basis)
    mooseError("Cannot changeConstraintToBulk for species ",
               basis_ind,
               " because there are only ",
               _num_basis,
               " basis species");
  if (_mgd.basis_species_gas[basis_ind])
    mooseError("Attempting to changeConstraintToBulk for a gas species.  Since a swap is involved, "
               "you cannot specify a value for the bulk number of moles.  You must use "
               "changeConstraintToBulk(basis_ind) method instead of "
               "changeConstraintToBulk(basis_ind, value)");
  if (basis_ind == 0)
    _constraint_meaning[basis_ind] = ConstraintMeaningEnum::MOLES_BULK_WATER;
  else
    _constraint_meaning[basis_ind] = ConstraintMeaningEnum::MOLES_BULK_SPECIES;
  setConstraintValue(basis_ind, value);
 
  // it is possible that FIXED_ACTIVITY just became MOLES_BULK_SPECIES
  buildKnownBasisActivities(_basis_activity_known, _basis_activity);
 
  // it is likely the algebraic system has changed
  buildAlgebraicInfo(_in_algebraic_system,
                     _num_basis_in_algebraic_system,
                     _num_in_algebraic_system,
                     _algebraic_index,
                     _basis_index);
}

◆ checkAndInitialize()

void GeochemicalSystem::checkAndInitialize	(	const std::vector< std::string > &	kin_name,
		const std::vector< Real > &	kin_initial,
		const std::vector< GeochemistryUnitConverter::GeochemistryUnit > &	kin_unit
	)

private

Used during construction: checks for sane inputs and initializes molalities, etc, using the initialize() method.

Parameters

kin_name	names of kinetic species
kin_initial	initial mole numbers (or mass or volume, depending on kin_unit) of the species named in kin_name
kin_unit	units of the numbers provided in kin_initial

Definition at line 157 of file GeochemicalSystem.C.

{
  // initialize every the kinetic species
  const unsigned num_kin_name = kin_name.size();
  if (!(_num_kin == num_kin_name && _num_kin == kin_initial.size() && _num_kin == kin_unit.size()))
    mooseError("Initial mole number (or mass or volume) and a unit must be provided for each "
               "kinetic species ",
               _num_kin,
               " ",
               num_kin_name,
               " ",
               kin_initial.size(),
               " ",
               kin_unit.size());
  for (const auto & name_index : _mgd.kin_species_index)
  {
    const unsigned ind = std::distance(
        kin_name.begin(), std::find(kin_name.begin(), kin_name.end(), name_index.first));
    if (ind < num_kin_name)
    {
      if (!(kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::MOLES ||
            kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::KG ||
            kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::G ||
            kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::MG ||
            kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::UG ||
            kin_unit[ind] == GeochemistryUnitConverter::GeochemistryUnit::CM3))
        mooseError("Kinetic species ",
                   name_index.first,
                   ": units must be moles or mass, or volume in the case of minerals");
      const Real moles = GeochemistryUnitConverter::toMoles(
          kin_initial[ind], kin_unit[ind], name_index.first, _mgd);
      setKineticMoles(name_index.second, moles);
    }
    else
      mooseError("Initial mole number or mass or volume for kinetic species ",
                 name_index.first,
                 " must be provided");
  }
 
  // check sanity of swaps
  if (_swap_out.size() != _swap_in.size())
    mooseError("swap_out_of_basis must have same length as swap_into_basis");
  for (unsigned i = 0; i < _swap_out.size(); ++i)
    if (_swap_out[i] == _charge_balance_species)
      mooseError("Cannot swap out ",
                 _charge_balance_species,
                 " because it is the charge-balance species\n");
 
  // do swaps desired by user.  any exception here is an error
  for (unsigned i = 0; i < _swap_out.size(); ++i)
    try
    {
      _swapper.performSwap(_mgd, _swap_out[i], _swap_in[i]);
    }
    catch (const MooseException & e)
    {
      mooseError(e.what());
    }
 
  // check charge-balance species is in the basis and has a charge
  if (_mgd.basis_species_index.count(_charge_balance_species) == 0)
    mooseError("Cannot enforce charge balance using ",
               _charge_balance_species,
               " because it is not in the basis");
  _charge_balance_basis_index = _mgd.basis_species_index.at(_charge_balance_species);
  if (_mgd.basis_species_charge[_charge_balance_basis_index] == 0.0)
    mooseError("Cannot enforce charge balance using ",
               _charge_balance_species,
               " because it has zero charge");
 
  // check that constraint vectors have appropriate sizes
  if (_constrained_species.size() != _constraint_value.size())
    mooseError("Constrained species names must have same length as constraint values");
  if (_constrained_species.size() != _constraint_unit.size())
    mooseError("Constrained species names must have same length as constraint units");
  if (_constrained_species.size() != _constraint_user_meaning.size())
    mooseError("Constrained species names must have same length as constraint meanings");
  if (_constrained_species.size() != _num_basis)
    mooseError("Constrained species names must have same length as the number of species in the "
               "basis (each component must be provided with a single constraint");
 
  // check that each _constrained_species name appears in the basis
  for (const auto & name : _mgd.basis_species_name)
    if (std::find(_constrained_species.begin(), _constrained_species.end(), name) ==
        _constrained_species.end())
      mooseError("The basis species ", name, " must appear in the constrained species list");
 
  // order the constraints in the same way as the basis species.  This makes the remainder of the
  // code much cleaner.
  std::vector<std::string> c_s(_num_basis);
  std::vector<Real> c_v(_num_basis);
  std::vector<GeochemistryUnitConverter::GeochemistryUnit> c_u(_num_basis);
  std::vector<ConstraintUserMeaningEnum> c_m(_num_basis);
  for (unsigned i = 0; i < _num_basis; ++i)
  {
    const unsigned basis_ind = _mgd.basis_species_index.at(_constrained_species[i]);
    c_s[basis_ind] = _constrained_species[i];
    c_v[basis_ind] = _constraint_value[i];
    c_u[basis_ind] = _constraint_unit[i];
    c_m[basis_ind] = _constraint_user_meaning[i];
  }
  _constrained_species = c_s;
  _constraint_value = c_v;
  _constraint_unit = c_u;
  _original_constraint_value = c_v;
  _constraint_user_meaning = c_m;
 
  // run through the constraints, checking physical and chemical consistency, converting to mole
  // units, and building constraint_meaning
  for (unsigned i = 0; i < _constrained_species.size(); ++i)
  {
    const std::string name = _constrained_species[i];
 
    switch (_constraint_user_meaning[i])
    {
      case ConstraintUserMeaningEnum::KG_SOLVENT_WATER:
      {
        // if the mass of solvent water is provided, check it is positive
        if (_constraint_value[i] <= 0.0)
          mooseError("Specified mass of solvent water must be positive: you entered ",
                     _constraint_value[i]);
        if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::KG)
          mooseError("Units for kg_solvent_water must be kg");
        _constraint_meaning[i] = ConstraintMeaningEnum::KG_SOLVENT_WATER;
        break;
      }
      case ConstraintUserMeaningEnum::BULK_COMPOSITION:
      case ConstraintUserMeaningEnum::BULK_COMPOSITION_WITH_KINETIC:
      {
        // convert to mole units and specify correct constraint_meaning
        _constraint_value[i] = GeochemistryUnitConverter::toMoles(
            _constraint_value[i], _constraint_unit[i], name, _mgd);
        // add contributions from kinetic moles, if necessary
        if (_constraint_user_meaning[i] == ConstraintUserMeaningEnum::BULK_COMPOSITION)
          for (unsigned k = 0; k < _mgd.kin_species_name.size(); ++k)
            _constraint_value[i] += _kin_moles[k] * _mgd.kin_stoichiometry(k, i);
        if (!(_constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MOLES ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::KG ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::G ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MG ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::UG))
          mooseError("Species ", name, ": units for bulk composition must be moles or mass");
        if (name == "H2O")
          _constraint_meaning[i] = ConstraintMeaningEnum::MOLES_BULK_WATER;
        else
          _constraint_meaning[i] = ConstraintMeaningEnum::MOLES_BULK_SPECIES;
        break;
      }
      case ConstraintUserMeaningEnum::FREE_CONCENTRATION:
      {
        // if free concentration, check it is positive and perform the translation to mole units
        if (_constraint_value[i] <= 0.0)
          mooseError("Specified free concentration values must be positive: you entered ",
                     _constraint_value[i]);
        if (!(_constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MOLAL ||
              _constraint_unit[i] ==
                  GeochemistryUnitConverter::GeochemistryUnit::KG_PER_KG_SOLVENT ||
              _constraint_unit[i] ==
                  GeochemistryUnitConverter::GeochemistryUnit::G_PER_KG_SOLVENT ||
              _constraint_unit[i] ==
                  GeochemistryUnitConverter::GeochemistryUnit::MG_PER_KG_SOLVENT ||
              _constraint_unit[i] ==
                  GeochemistryUnitConverter::GeochemistryUnit::UG_PER_KG_SOLVENT))
          mooseError(
              "Species ",
              name,
              ": units for free concentration quantities must be molal or mass_per_kg_solvent");
        _constraint_value[i] = GeochemistryUnitConverter::toMoles(
            _constraint_value[i], _constraint_unit[i], name, _mgd);
        _constraint_meaning[i] = ConstraintMeaningEnum::FREE_MOLALITY;
        break;
      }
      case ConstraintUserMeaningEnum::FREE_MINERAL:
      {
        // if free mineral, check it is positive and perform the translation to mole units
        if (_constraint_value[i] <= 0.0)
          mooseError("Specified free mineral values must be positive: you entered ",
                     _constraint_value[i]);
        if (!(_constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MOLES ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::KG ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::G ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::MG ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::UG ||
              _constraint_unit[i] == GeochemistryUnitConverter::GeochemistryUnit::CM3))
          mooseError("Species ",
                     name,
                     ": units for free mineral quantities must be moles, mass or volume");
        _constraint_value[i] = GeochemistryUnitConverter::toMoles(
            _constraint_value[i], _constraint_unit[i], name, _mgd);
        _constraint_meaning[i] = ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES;
        break;
      }
      case ConstraintUserMeaningEnum::ACTIVITY:
      {
        if (_constraint_value[i] <= 0.0)
          mooseError("Specified activity values must be positive: you entered ",
                     _constraint_value[i]);
        if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
          mooseError(
              "Species ", name, ": dimensionless units must be used when specifying activity");
        _constraint_meaning[i] = ConstraintMeaningEnum::ACTIVITY;
        break;
      }
      case ConstraintUserMeaningEnum::LOG10ACTIVITY:
      {
        // if log10activity is provided, convert it to activity
        if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
          mooseError("Species ",
                     name,
                     ": dimensionless units must be used when specifying log10activity\n");
        _constraint_value[i] = std::pow(10.0, _constraint_value[i]);
        _constraint_meaning[i] = ConstraintMeaningEnum::ACTIVITY;
        break;
      }
      case ConstraintUserMeaningEnum::FUGACITY:
      {
        // if fugacity is provided, check it is positive
        if (_constraint_value[i] <= 0.0)
          mooseError("Specified fugacity values must be positive: you entered ",
                     _constraint_value[i]);
        if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
          mooseError(
              "Species ", name, ": dimensionless units must be used when specifying fugacity\n");
        _constraint_meaning[i] = ConstraintMeaningEnum::FUGACITY;
        break;
      }
      case ConstraintUserMeaningEnum::LOG10FUGACITY:
      {
        // if log10fugacity is provided, convert it to fugacity
        if (_constraint_unit[i] != GeochemistryUnitConverter::GeochemistryUnit::DIMENSIONLESS)
          mooseError("Species ",
                     name,
                     ": dimensionless units must be used when specifying log10fugacity\n");
        _constraint_value[i] = std::pow(10.0, _constraint_value[i]);
        _constraint_meaning[i] = ConstraintMeaningEnum::FUGACITY;
        break;
      }
    }
 
    // check that water is provided with correct meaning
    if (name == "H2O")
      if (!(_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_WATER ||
            _constraint_meaning[i] == ConstraintMeaningEnum::KG_SOLVENT_WATER ||
            _constraint_meaning[i] == ConstraintMeaningEnum::ACTIVITY))
        mooseError("H2O must be provided with either a mass of solvent water, a bulk composition "
                   "in moles or mass, or an activity");
 
    // check that gases are provided with the correct meaning
    if (_mgd.basis_species_gas[i])
      if (_constraint_meaning[i] != ConstraintMeaningEnum::FUGACITY)
        mooseError("The gas ", name, " must be provided with a fugacity");
 
    // check that minerals are provided with the correct meaning
    if (_mgd.basis_species_mineral[i])
      if (!(_constraint_meaning[i] == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES ||
            _constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES))
        mooseError("The mineral ",
                   name,
                   " must be provided with either: a free number of moles, a free mass or a free "
                   "volume; or a bulk composition of moles or mass");
 
    // check that non-water, non-minerals, non-gases are provided with the correct meaning
    if (name != "H2O" && !_mgd.basis_species_gas[i] && !_mgd.basis_species_mineral[i])
      if (!(_constraint_meaning[i] == ConstraintMeaningEnum::FREE_MOLALITY ||
            _constraint_meaning[i] == ConstraintMeaningEnum::ACTIVITY ||
            _constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES))
        mooseError("The basis species ",
                   name,
                   " must be provided with a free concentration, bulk composition or an activity");
 
    // check that the charge-balance species has been provided MOLES_BULK_SPECIES
    if (name == _charge_balance_species)
      if (_constraint_meaning[i] != ConstraintMeaningEnum::MOLES_BULK_SPECIES)
        mooseError("For code consistency, the species ",
                   name,
                   " must be provided with a bulk composition because it is the charge-balance "
                   "species.  The value provided should be a reasonable estimate of the mole "
                   "number, but will be overridden as the solve progresses");
  }
  _original_redox_lhs = _mgd.redox_lhs;
 
  initialize();
}

Referenced by GeochemicalSystem(), and GeochemicalSystem().

◆ closeSystem()

void GeochemicalSystem::closeSystem ( )

Changes a KG_SOLVENT_WATER constraint to MOLES_BULK_WATER (if there is such a constraint) and all FREE_MOLALITY and FREE_MOLES_MINERAL_SPECIES to MOLES_BULK_SPECIES using the current values of _bulk_moles_old.

If, after performing these changes, all constraints are MOLES_BULK_WATER or MOLES_BULK_SPECIES then charge-balance is performed by altering the bulk_moles_old for the charge-balance species.

Definition at line 2100 of file GeochemicalSystem.C.

{
  for (unsigned basis_ind = 0; basis_ind < _num_basis; ++basis_ind)
    if (_constraint_meaning[basis_ind] == ConstraintMeaningEnum::KG_SOLVENT_WATER ||
        _constraint_meaning[basis_ind] == ConstraintMeaningEnum::FREE_MOLALITY ||
        _constraint_meaning[basis_ind] == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES)
      changeConstraintToBulk(basis_ind);
}

Referenced by GeochemistryTimeDependentReactor::execute(), TEST(), TEST(), and TEST_F().

◆ computeAndGetEquilibriumActivity()

const std::vector< Real > & GeochemicalSystem::computeAndGetEquilibriumActivity ( )

Computes activity for all equilibrium species (_eqm_activity) and returns a reference to the vector.

The vector _eqm_activity is only used for output purposes (such as recording the fugacity of equilibrium gases) so it is only computed by this method and not internally during computeConsistentConfiguration.

Returns: equilibrium activities

Definition at line 2342 of file GeochemicalSystem.C.

{
  for (unsigned j = 0; j < _num_eqm; ++j)
    _eqm_activity[j] = getEquilibriumActivity(j);
  return _eqm_activity;
}

Referenced by TEST_F().

◆ computeBulk()

void GeochemicalSystem::computeBulk ( std::vector< Real > & bulk_moles_old ) const

private

Computes the value of bulk_moles_old for all basis species.

For species with BULK-type constraints bulk_moles_old gets set to the constraint value, while for other species bulk_moles_old is computed from the molalities, and, in contrast to Bethke, the kinetic species

Definition at line 1105 of file GeochemicalSystem.C.

{
  for (unsigned i = 0; i < _num_basis; ++i)
  {
    const Real value = _constraint_value[i];
    const ConstraintMeaningEnum meaning = _constraint_meaning[i];
    switch (meaning)
    {
      case ConstraintMeaningEnum::MOLES_BULK_SPECIES:
      case ConstraintMeaningEnum::MOLES_BULK_WATER:
      {
        bulk_moles_old[i] = value;
        break;
      }
      case ConstraintMeaningEnum::KG_SOLVENT_WATER:
      case ConstraintMeaningEnum::FREE_MOLALITY:
      case ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES:
      case ConstraintMeaningEnum::FUGACITY:
      case ConstraintMeaningEnum::ACTIVITY:
      {
        bulk_moles_old[i] = computeBulkFromMolalities(i);
        break;
      }
    }
  }
}

Referenced by computeConsistentConfiguration(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and updateOldWithCurrent().

◆ computeBulkFromMolalities()

Real GeochemicalSystem::computeBulkFromMolalities ( unsigned basis_ind ) const

private

Parameters

basis_ind the index of the basis species in mgd

Returns: the number of bulk moles of the species basis_ind. This depends on the current molality of the basis_ind species, as well as the current molalities of the equilibrium species that depend on this basis species, and, in contrast to the Bethke approach, the current mole number of kinetic species that depend on this basis species.

Definition at line 1133 of file GeochemicalSystem.C.

{
  const Real nw = _basis_molality[0];
  Real bulk = 0.0;
  if (basis_ind == 0)
    bulk = GeochemistryConstants::MOLES_PER_KG_WATER;
  else if (_mgd.basis_species_mineral[basis_ind])
    bulk = _basis_molality[basis_ind] / nw; // because of multiplication by nw, below
  else
    bulk = _basis_molality[basis_ind];
  for (unsigned j = 0; j < _num_eqm; ++j)
    bulk += _mgd.eqm_stoichiometry(j, basis_ind) * _eqm_molality[j];
  bulk *= nw;
  for (unsigned kin = 0; kin < _num_kin; ++kin)
    bulk += _mgd.kin_stoichiometry(kin, basis_ind) * _kin_moles[kin];
  return bulk;
}

Referenced by changeConstraintToBulk(), computeBulk(), getResidualComponent(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ computeConsistentConfiguration()

void GeochemicalSystem::computeConsistentConfiguration ( )

Compute a reasonably consistent configuration that can be used as an initial condition in a Newton process to solve the algebraic system.

Note that unless _iters_to_make_consistent is very large, the resulting configuration will not be completely consistent, but this is conventional in geochemical modelling: there are so many unknowns and approximations used in the Newton process, that starting from a completely consistent configuration is unimportant.

Before entering this method, it is assumed that::

basis molality has been set for all basis species. For species where this is unknown (viz, it will be solved for in the Newton process) a reasonable initial condition should be set. Use the initBulkAndFree method.
basis_activity_known has been set to true for basis species whose activities are specified by the user, for basis gases and for basis minerals. Use the buildKnownBasisActivities method.
basis_activity has been set for basis species with basis_activity_known = true. Use the buildKnownBasisActivities method.
equilibrium molality has been pre-initialised (eg, to zero) for all equilibrium species prior to entering this method.

This method computes:

ionic strength and stoichiometric ionic strength
basis activity coefficients (not for water, minerals or gases)
equilibrium activity coefficients (not for minerals or gases)
basis molality for basis species where the user has specified the activity
basis molality for mineral basis species
basis activity for all basis species
equilibrium molality for all equilibrium species (molality for minerals and gases are set to zero)
bulk_moles_old for all basis species (set to the constraint_value for BULK-type constraints, otherwise computed from the current values of molality)
sorbing surface area

It does not compute activity for equilibrium species: use computeAndGetEquilibriumActivity for that

Definition at line 464 of file GeochemicalSystem.C.

{
  // the steps 1 and 2 below could be iterated for a long time (or a Newton process could even be
  // followed) to provide better estimates of activities and molalities, but this is not done in the
  // conventional geochemistry approach: there are just too many unknowns and approximations
  // employed during the algebraic-system solve to justify iterating towards the perfectly
  // consistent initial condition
  for (unsigned picard = 0; picard < _iters_to_make_consistent + 1; ++picard)
  {
    // Step 1: compute ionic strengths and activities using the eqm molalities
    _gac.setInternalParameters(_temperature, _mgd, _basis_molality, _eqm_molality, _kin_moles);
    _gac.buildActivityCoefficients(_mgd, _basis_activity_coef, _eqm_activity_coef);
    updateBasisMolalityForKnownActivity(_basis_molality);
    computeRemainingBasisActivities(_basis_activity);
 
    // Step 2: compute equilibrium molality based on the activities just computed
    computeEqmMolalities(_eqm_molality);
  }
 
  computeBulk(_bulk_moles_old);
  computeFreeMineralMoles(_basis_molality);
  computeSorbingSurfaceArea(_sorbing_surface_area);
}

Referenced by GeochemistryTimeDependentReactor::execute(), initialize(), performSwapNoCheck(), GeochemistryTimeDependentReactor::postSolveExchanger(), setAlgebraicVariables(), and setConstraintValue().

◆ computeEqmMolalities()

void GeochemicalSystem::computeEqmMolalities ( std::vector< Real > & eqm_molality ) const

private

compute the equilibrium molalities (not for minerals or gases)

Definition at line 1007 of file GeochemicalSystem.C.

{
  for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
  {
    if (_mgd.eqm_species_mineral[eqm_j] || _mgd.eqm_species_gas[eqm_j])
      eqm_molality[eqm_j] = 0.0;
    else
    {
      // compute log10 version first, in an attempt to eliminate overflow and underflow problems
      // such as 10^(1000)
      const Real log10m = log10ActivityProduct(eqm_j) - _eqm_log10K[eqm_j];
      eqm_molality[eqm_j] =
          std::pow(10.0, log10m) / _eqm_activity_coef[eqm_j] * surfaceSorptionModifier(eqm_j);
    }
  }
}

Referenced by computeConsistentConfiguration().

◆ computeFreeMineralMoles()

void GeochemicalSystem::computeFreeMineralMoles ( std::vector< Real > & basis_molality ) const

private

Compute the free mineral moles (ie, basis_molality for basis species that are minerals), using the bulk_moles_old.

Note that usually computeBulk should have been called immediately preceding this in order to correctly set bulk_moles_old.

Definition at line 1577 of file GeochemicalSystem.C.

{
 
  const Real nw = _basis_molality[0];
  for (unsigned i = 0; i < _num_basis; ++i)
    if (_mgd.basis_species_mineral[i])
    {
      if (_constraint_meaning[i] == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES)
        basis_molality[i] = _constraint_value[i];
      else
      {
        basis_molality[i] = _bulk_moles_old[i];
        for (unsigned j = 0; j < _num_eqm; ++j)
          basis_molality[i] -= nw * _mgd.eqm_stoichiometry(j, i) * _eqm_molality[j];
        for (unsigned kin = 0; kin < _num_kin; ++kin)
          basis_molality[i] -= _mgd.kin_stoichiometry(kin, i) * _kin_moles[kin];
      }
    }
}

Referenced by alterSystemBecauseBulkChanged(), and computeConsistentConfiguration().

◆ computeJacobian()

void GeochemicalSystem::computeJacobian	(	const DenseVector< Real > &	res,
		DenseMatrix< Real > &	jac,
		const DenseVector< Real > &	mole_additions,
		const DenseMatrix< Real > &	dmole_additions
	)		const

Compute the Jacobian for the algebraic system and put it in jac.

Parameters

res	The residual of the algebraic system. This is used to speed up computations of the jacobian
jac	The jacobian entries are placed here
mole_additions	The molar additions to the basis species and the kinetic species
dmole_additions	d(mole_additions)/d(molality or kinetic_moles)

Definition at line 1287 of file GeochemicalSystem.C.

{
  /* To the reader: yes, this is awfully complicated.  The best way of understanding it is to very
   * slowly go through the residual calculation and compute the derivatives by hand.  It is quite
   * well tested, but it'd be great if you could add more tests!
   */
  if (res.size() != _num_in_algebraic_system)
    mooseError(
        "Jacobian: residual size must be ", _num_in_algebraic_system, " but it is ", res.size());
  if (mole_additions.size() != _num_basis + _num_kin)
    mooseError("Jacobian: the increment in mole numbers (mole_additions) needs to be of size ",
               _num_basis,
               " + ",
               _num_kin,
               " but it is of size ",
               mole_additions.size());
  if (!(dmole_additions.n() == _num_basis + _num_kin &&
        dmole_additions.m() == _num_basis + _num_kin))
    mooseError("Jacobian: the derivative of mole additions (dmole_additions) needs to be of size ",
               _num_basis + _num_kin,
               "x",
               _num_basis + _num_kin,
               " but it is of size ",
               dmole_additions.m(),
               "x",
               dmole_additions.n());
  /*
Note that in this function we make use of the fact that species can only be in the algebraic
system if their molality is unknown (or in the case of water, the mass of solvent mass is
unknown).  This means that the molality of the equilibrium species depends on
(activity_coefficient * basis molality)^stoi, so that the derivative is quite simple.  Eg, we
don't have to worry about derivatives with respect to fixed-activity things, or fixed fugacity.
Also, note that the constructor and the various "set" methods of this class enforce molality > 0,
so there are no division-by-zero problems.  Also, note that we never compute derivatives of the
activity coefficients, or the activity of water with respect to the molalities, as they are
assumed to be quite small.  Finally, the surface_pot_expr will also always be positive.
  */
  // correctly size and zero
  jac.resize(_num_in_algebraic_system, _num_in_algebraic_system);
  const Real nw = _basis_molality[0];
 
  // jac(a, b) = d(R_a) / d(m_b), where a corresponds to a molality
  for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
  {
    const unsigned basis_of_a = _basis_index[a];
    for (unsigned b = 0; b < _num_basis_in_algebraic_system; ++b)
    {
      const unsigned basis_of_b = _basis_index[b];
 
      // contribution from mole_additions
      if (basis_of_a != _charge_balance_basis_index)
        jac(a, b) -= dmole_additions(basis_of_a, basis_of_b);
      else
      {
        for (unsigned i = 0; i < _num_basis; ++i)
        {
          if (i == _charge_balance_basis_index)
            continue;
          else if (_mgd.basis_species_charge[i] ==
                   0.0) // certainly includes water, minerals and gases
            continue;
          else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
            jac(a, b) += _mgd.basis_species_charge[i] * dmole_additions(i, basis_of_b) /
                         _mgd.basis_species_charge[basis_of_a];
        }
      }
 
      // contribution from explicit nw, in case where mass of solvent water is an unknown
      if (basis_of_b == 0)
      {
        if (basis_of_a != _charge_balance_basis_index)
        {
          // use a short-cut here: dR/dnw = m + sum_eqm[stoi*mol] = (R + bulk + addition - kin)/nw
          Real numerator = res(a) + _bulk_moles_old[basis_of_a] + mole_additions(basis_of_a);
          for (unsigned kin = 0; kin < _num_kin; ++kin)
            numerator -= _mgd.kin_stoichiometry(kin, basis_of_a) * _kin_moles[kin];
          jac(a, 0) += numerator / nw;
        }
        else
        {
          for (unsigned i = 0; i < _num_basis; ++i)
          {
            if (i == _charge_balance_basis_index)
              continue;
            else if (_mgd.basis_species_charge[i] ==
                     0.0) // certainly includes water, minerals and gases
              continue;
            else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
              continue;
            else
            {
              Real molal = _basis_molality[i];
              for (unsigned j = 0; j < _num_eqm; ++j)
                molal += _mgd.eqm_stoichiometry(j, i) * _eqm_molality[j];
              jac(a, 0) +=
                  _mgd.basis_species_charge[i] * molal / _mgd.basis_species_charge[basis_of_a];
            }
          }
          jac(a, 0) += _basis_molality[basis_of_a];
          for (unsigned j = 0; j < _num_eqm; ++j)
            jac(a, 0) += _mgd.eqm_stoichiometry(j, basis_of_a) * _eqm_molality[j];
        }
      }
 
      // contribution from molality of basis_of_b
      if (basis_of_b != 0)
      {
        if (a == b)
          jac(a, b) += nw; // remember b != 0
        for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
        {
          jac(a, b) += nw * _mgd.eqm_stoichiometry(eqm_j, basis_of_a) * _eqm_molality[eqm_j] *
                       _mgd.eqm_stoichiometry(eqm_j, basis_of_b) / _basis_molality[basis_of_b];
        }
        if (basis_of_a == _charge_balance_basis_index)
        {
          // additional terms from the charge-balance additions
          for (unsigned i = 0; i < _num_basis; ++i)
          {
            if (i == _charge_balance_basis_index)
              continue;
            else if (_mgd.basis_species_charge[i] ==
                     0.0) // certainly includes water, minerals and gases
              continue;
            else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
              continue;
            else
            {
              const Real prefactor =
                  _mgd.basis_species_charge[i] * nw / _mgd.basis_species_charge[basis_of_a];
              if (i == basis_of_b)
                jac(a, b) += prefactor;
              for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
              {
                jac(a, b) += prefactor * _mgd.eqm_stoichiometry(eqm_j, i) * _eqm_molality[eqm_j] *
                             _mgd.eqm_stoichiometry(eqm_j, basis_of_b) /
                             _basis_molality[basis_of_b];
              }
            }
          }
        }
      }
    }
  }
 
  // jac(a, b) = d(R_a) / d(surface_pot), where a corresponds to a molality
  for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
  {
    const unsigned basis_of_a = _basis_index[a];
    for (unsigned s = 0; s < _num_surface_pot; ++s)
    {
      const unsigned b = s + _num_basis_in_algebraic_system;
      // derivative of nw * _mgd.eqm_stoichiometry(eqm_j, basis_i) * _eqm_molality[eqm_j],
      // where _eqm_molality = (_surface_pot_expr)^(2 * charge) * stuff
      for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
        if (_mgd.surface_sorption_related[eqm_j] && _mgd.surface_sorption_number[eqm_j] == s)
          jac(a, b) += nw * _mgd.eqm_stoichiometry(eqm_j, basis_of_a) * 2.0 *
                       _mgd.eqm_species_charge[eqm_j] * _eqm_molality[eqm_j] /
                       _surface_pot_expr[_mgd.surface_sorption_number[eqm_j]];
    }
  }
 
  // jac(a, b) = d(R_a) / d(_kin_moles) where a corresponds to a molality
  for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
  {
    const unsigned basis_of_a = _basis_index[a];
 
    // contribution from mole_additions
    for (unsigned kin = 0; kin < _num_kin; ++kin)
    {
      const unsigned ind = _num_basis + kin;
      const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kin;
      if (basis_of_a != _charge_balance_basis_index)
        jac(a, b) -= dmole_additions(basis_of_a, ind);
      else
      {
        for (unsigned i = 0; i < _num_basis; ++i)
        {
          if (i == _charge_balance_basis_index)
            continue;
          else if (_mgd.basis_species_charge[i] ==
                   0.0) // certainly includes water, minerals and gases
            continue;
          else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
            jac(a, b) += _mgd.basis_species_charge[i] * dmole_additions(i, ind) /
                         _mgd.basis_species_charge[basis_of_a];
        }
      }
    }
 
    // contribution from sum_kin[stoi*kin_moles]
    for (unsigned kin = 0; kin < _num_kin; ++kin)
    {
      const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kin;
      jac(a, b) += _mgd.kin_stoichiometry(kin, basis_of_a);
    }
 
    // special additional contribution for charge-balance from kinetic stuff
    if (basis_of_a == _charge_balance_basis_index)
    {
      for (unsigned i = 0; i < _num_basis; ++i)
      {
        if (i == _charge_balance_basis_index)
          continue;
        else if (_mgd.basis_species_charge[i] ==
                 0.0) // certainly includes water, minerals and gases
          continue;
        else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
          continue;
        else
        {
          const Real prefactor =
              _mgd.basis_species_charge[i] / _mgd.basis_species_charge[basis_of_a];
          for (unsigned kin = 0; kin < _num_kin; ++kin)
          {
            const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kin;
            jac(a, b) += prefactor * _mgd.kin_stoichiometry(kin, i);
          }
        }
      }
    }
  }
 
  // jac(a, b) = d(R_a) / d(variable_b) where a corresponds to a surface potential
  for (unsigned s = 0; s < _num_surface_pot; ++s)
  {
    const unsigned a = s + _num_basis_in_algebraic_system;
 
    for (unsigned b = 0; b < _num_basis_in_algebraic_system; ++b)
    {
      // derivative of nw * _mgd.eqm_species_charge[j] * _eqm_molality[j];
      const unsigned basis_of_b = _basis_index[b];
      if (basis_of_b == 0) // special case: mass of solvent water is an unknown
      {
        for (unsigned j = 0; j < _num_eqm; ++j)
          if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == s)
            jac(a, b) += _mgd.eqm_species_charge[j] * _eqm_molality[j];
      }
      else
      {
        for (unsigned j = 0; j < _num_eqm; ++j)
          if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == s)
            jac(a, b) += nw * _mgd.eqm_species_charge[j] * _eqm_molality[j] *
                         _mgd.eqm_stoichiometry(j, basis_of_b) / _basis_molality[basis_of_b];
      }
    }
 
    const Real coef = surfacePotPrefactor(s);
    // derivative of coef * (x - 1/x) wrt x, where x = _surface_pot_expr
    jac(a, a) += coef * (1.0 + 1.0 / std::pow(_surface_pot_expr[s], 2.0));
    // derivative of nw * _mgd.eqm_species_charge[j] * _eqm_molality[j];
    // where _eqm_molality = (_surface_pot_expr)^(2 * charge) * stuff
    for (unsigned j = 0; j < _num_eqm; ++j)
      if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == s)
        jac(a, a) += nw * std::pow(_mgd.eqm_species_charge[j], 2.0) * 2.0 * _eqm_molality[j] /
                     _surface_pot_expr[s];
  }
 
  // jac(a, b) = d(R_a) / d(variable_b) where a corresponds to a kinetic species
  for (unsigned kin = 0; kin < _num_kin; ++kin)
  {
    const unsigned a = _num_basis_in_algebraic_system + _num_surface_pot + kin;
    jac(a, a) += 1.0; // deriv wrt _kin_moles[kin]
 
    // contribution from mole_additions
    const unsigned ind_of_addition = _num_basis + kin;
    for (unsigned b = 0; b < _num_basis_in_algebraic_system; ++b)
    {
      const unsigned basis_of_b = _basis_index[b];
      jac(a, b) -= dmole_additions(ind_of_addition, basis_of_b);
    }
    for (unsigned kinp = 0; kinp < _num_kin; ++kinp)
    {
      const unsigned b = _num_basis_in_algebraic_system + _num_surface_pot + kinp;
      const unsigned indp = _num_basis + kinp;
      jac(a, b) -= dmole_additions(ind_of_addition, indp);
    }
  }
}

Referenced by GeochemicalSolver::solveSystem().

◆ computeRemainingBasisActivities()

void GeochemicalSystem::computeRemainingBasisActivities ( std::vector< Real > & basis_activity ) const

private

Compute the activity of water and put in basis_activity[0], and use the _basis_activity_coef and _basis_molality to compute the remaining basis activities (for those that are not minerals, gases or aqueous species provided with an activity)

Definition at line 939 of file GeochemicalSystem.C.

{
  if (!_basis_activity_known[0])
    basis_activity[0] = _gac.waterActivity();
  for (unsigned basis_ind = 1; basis_ind < _num_basis; ++basis_ind) // don't loop over water
    if (!_basis_activity_known[basis_ind]) // basis_activity_known = true for minerals, gases and
                                           // species with activities provided by the user
      basis_activity[basis_ind] = _basis_activity_coef[basis_ind] * _basis_molality[basis_ind];
}

Referenced by computeConsistentConfiguration(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ computeSorbingSurfaceArea()

void GeochemicalSystem::computeSorbingSurfaceArea ( std::vector< Real > & sorbing_surface_area ) const

private

Compute sorbing_surface_area depending on the current molality of the sorbing minerals.

Definition at line 1842 of file GeochemicalSystem.C.

{
  for (unsigned sp = 0; sp < _num_surface_pot; ++sp)
  {
    sorbing_surface_area[sp] = _mgd.surface_sorption_area[sp];
    if (_mgd.basis_species_index.count(_mgd.surface_sorption_name[sp]) == 1)
    {
      const unsigned basis_ind = _mgd.basis_species_index.at(_mgd.surface_sorption_name[sp]);
      const Real grams =
          _mgd.basis_species_molecular_weight[basis_ind] * _basis_molality[basis_ind];
      sorbing_surface_area[sp] *= grams;
    }
  }
}

Referenced by computeConsistentConfiguration(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ computeTransportedBulkFromMolalities()

void GeochemicalSystem::computeTransportedBulkFromMolalities ( std::vector< Real > & transported_bulk ) const

Computes the value of transported bulk moles for all basis species using the existing molalities.

Note that this is probably gives rubbish results unless the system is consistent (eg, the solve has converged). Also note that this does not include contributions from kinetic species

Definition at line 1152 of file GeochemicalSystem.C.

{
  transported_bulk.resize(_num_basis);
  const Real nw = _basis_molality[0];
  for (unsigned i = 0; i < _num_basis; ++i)
  {
    transported_bulk[i] = 0.0;
    if (i == 0)
      transported_bulk[i] = GeochemistryConstants::MOLES_PER_KG_WATER;
    else if (_mgd.basis_species_transported[i])
    {
      if (_mgd.basis_species_mineral[i])
        transported_bulk[i] = _basis_molality[i] / nw; // because of multiplication by nw, below
      else
        transported_bulk[i] = _basis_molality[i];
    }
    else
      transported_bulk[i] = 0.0;
    for (unsigned j = 0; j < _num_eqm; ++j)
      if (_mgd.eqm_species_transported[j])
        transported_bulk[i] += _mgd.eqm_stoichiometry(j, i) * _eqm_molality[j];
    transported_bulk[i] *= nw;
  }
}

Referenced by getTransportedBulkInOriginalBasis().

◆ enforceChargeBalance() [1/2]

void GeochemicalSystem::enforceChargeBalance ( )

Enforces charge balance by altering the constraint_value and bulk_moles_old of the charge-balance species.

Use with caution, since this overwrites the constraint values provided in the constructor, and also changes bulk_moles_old which is usually the value of bulk moles from the previous time-step

Definition at line 1069 of file GeochemicalSystem.C.

{
  enforceChargeBalance(_constraint_value, _bulk_moles_old);
}

Referenced by enforceChargeBalance(), GeochemicalSolver::solveSystem(), and TEST_F().

◆ enforceChargeBalance() [2/2]

void GeochemicalSystem::enforceChargeBalance	(	std::vector< Real > &	constraint_value,
		std::vector< Real > &	bulk_moles_old
	)		const

private

Enforces charge balance by altering the constraint_value and bulk_moles_old of the charge-balance species.

Note that this changes the important constraint_value and the old bulk_moles_old (usually from the previous time-step) so should be used with caution.

Definition at line 1075 of file GeochemicalSystem.C.

{
  const Real tot_charge = getTotalChargeOld();
  constraint_value[_charge_balance_basis_index] -=
      tot_charge / _mgd.basis_species_charge[_charge_balance_basis_index];
  bulk_moles_old[_charge_balance_basis_index] = constraint_value[_charge_balance_basis_index];
}

◆ enforceChargeBalanceIfSimple()

void GeochemicalSystem::enforceChargeBalanceIfSimple	(	std::vector< Real > &	constraint_value,
		std::vector< Real > &	bulk_moles_old
	)		const

private

If all non-charged basis species are provided with a bulk number of moles, alter the bulk number of moles (old) of the charge_balance_species so that charges are balanced.

Parameters

constraint_value	the values of the constraints provided by the user
bulk_moles_old	the old values of bulk moles (from previous time-step)

Definition at line 1036 of file GeochemicalSystem.C.

{
  Real tot_charge = 0.0;
  for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
    if (_mgd.basis_species_charge[basis_i] != 0.0)
    {
      if (_constraint_meaning[basis_i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
        tot_charge += _mgd.basis_species_charge[basis_i] * constraint_value[basis_i];
      else
        return;
    }
  // kinetic species are counted in the bulk
  // all charged basis species must have been provided with a MOLES_BULK_SPECIES value, so we can
  // easily enforce charge neutrality
  tot_charge -= _mgd.basis_species_charge[_charge_balance_basis_index] *
                constraint_value[_charge_balance_basis_index];
  constraint_value[_charge_balance_basis_index] =
      -tot_charge / _mgd.basis_species_charge[_charge_balance_basis_index];
  bulk_moles_old[_charge_balance_basis_index] = constraint_value[_charge_balance_basis_index];
}

Referenced by alterSystemBecauseBulkChanged(), initialize(), performSwapNoCheck(), setChargeBalanceSpecies(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ getAlgebraicBasisValues()

std::vector< Real > GeochemicalSystem::getAlgebraicBasisValues ( ) const

Returns: v, a vector of length _num_basis_in_algebraic_system, the elements of which are the current values of the algebraic variables: molalities only (not surface potentials or kinetic species)

Definition at line 696 of file GeochemicalSystem.C.

{
  std::vector<Real> var(_num_basis_in_algebraic_system);
  for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
    var[a] = _basis_molality[_basis_index[a]];
  return var;
}

Referenced by GeochemicalSolver::reduceInitialResidual().

◆ getAlgebraicIndexOfBasisSystem()

const std::vector< unsigned > & GeochemicalSystem::getAlgebraicIndexOfBasisSystem ( ) const

Returns: v, a vector of length _num_basis, where v[i] = the algebraic index of the i^th basis species. Note that unless the i^th basis species is in the algebraic system, this is meaningless

Definition at line 677 of file GeochemicalSystem.C.

{
  return _algebraic_index;
}

◆ getAlgebraicVariableDenseValues()

DenseVector< Real > GeochemicalSystem::getAlgebraicVariableDenseValues ( ) const

Returns: v, a DenseVector of length getNumInAlgebraicSystem, the elements of which are the current values of the algebraic variables

Definition at line 705 of file GeochemicalSystem.C.

{
  DenseVector<Real> var(_num_in_algebraic_system);
  for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
    var(a) = _basis_molality[_basis_index[a]];
  for (unsigned s = 0; s < _num_surface_pot; ++s)
    var(s + _num_basis_in_algebraic_system) = _surface_pot_expr[s];
  for (unsigned k = 0; k < _num_kin; ++k)
    var(k + _num_basis_in_algebraic_system + _num_surface_pot) = _kin_moles[k];
  return var;
}

◆ getAlgebraicVariableValues()

std::vector< Real > GeochemicalSystem::getAlgebraicVariableValues ( ) const

Returns: v, a vector of length getNumInAlgebraicSystem, the elements of which are the current values of the algebraic variables

Definition at line 683 of file GeochemicalSystem.C.

{
  std::vector<Real> var(_num_basis_in_algebraic_system + _num_surface_pot + _num_kin);
  for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
    var[a] = _basis_molality[_basis_index[a]];
  for (unsigned s = 0; s < _num_surface_pot; ++s)
    var[s + _num_basis_in_algebraic_system] = _surface_pot_expr[s];
  for (unsigned k = 0; k < _num_kin; ++k)
    var[k + _num_basis_in_algebraic_system + _num_surface_pot] = _kin_moles[k];
  return var;
}

Referenced by GeochemicalSolver::reduceInitialResidual(), and GeochemicalSolver::solveAndUnderrelax().

◆ getBasisActivity() [1/2]

const std::vector< Real > & GeochemicalSystem::getBasisActivity ( ) const

Returns: the activity for the basis species

Definition at line 882 of file GeochemicalSystem.C.

{
  return _basis_activity;
}

◆ getBasisActivity() [2/2]

Real GeochemicalSystem::getBasisActivity ( unsigned i ) const

Returns: the activity for the i^th basis species

Definition at line 870 of file GeochemicalSystem.C.

{
  if (i >= _num_basis)
    mooseError("Cannot retrieve basis activity for species ",
               i,
               " since there are only ",
               _num_basis,
               " basis species");
  return _basis_activity[i];
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryConsoleOutput::output(), and TEST_F().

◆ getBasisActivityCoefficient() [1/2]

const std::vector< Real > & GeochemicalSystem::getBasisActivityCoefficient ( ) const

Returns: the activity coefficients for the basis species. These are not defined for water, minerals, gases, or aqueous species that have been provided an activity by the user

Definition at line 980 of file GeochemicalSystem.C.

{
  return _basis_activity_coef;
}

◆ getBasisActivityCoefficient() [2/2]

Real GeochemicalSystem::getBasisActivityCoefficient ( unsigned i ) const

Returns: the activity coefficient for the i^th basis species. This is not defined for water, minerals, gases, or aqueous species that have been provided an activity by the user

Definition at line 968 of file GeochemicalSystem.C.

{
  if (i >= _num_basis)
    mooseError("Cannot retrieve basis activity coefficient for species ",
               i,
               " since there are only ",
               _num_basis,
               " basis species");
  return _basis_activity_coef[i];
}

Referenced by GeochemistryConsoleOutput::output().

◆ getBasisActivityKnown()

const std::vector< bool > & GeochemicalSystem::getBasisActivityKnown ( ) const

Returns: vector v, where v[i] = true if the activity of basis species i is fixed, either by the user providing an activity value, or a fugacity value, or because the i^th basis species is a mineral

Definition at line 864 of file GeochemicalSystem.C.

{
  return _basis_activity_known;
}

Referenced by GeochemicalSolver::swapNeeded().

◆ getBasisIndexOfAlgebraicSystem()

const std::vector< unsigned > & GeochemicalSystem::getBasisIndexOfAlgebraicSystem ( ) const

Returns: v, a vector of length _num_basis, where v[i] = the basis index of the i^th species in the algebraic system. Note that for i > _num_basis_in_algebraic_system, the value of v is undefined.

Definition at line 671 of file GeochemicalSystem.C.

{
  return _basis_index;
}

◆ getBulkMolesOld()

const std::vector< Real > & GeochemicalSystem::getBulkMolesOld ( ) const

Returns: the number of bulk-composition moles. Note this will typically be the "old" number of bulk moles (from the previous time-step) unless addtoBulkMoles or updateOldWithCurrent or other similar methods have just been called. Note that this contains contributions from kinetic species, in contrast to the Bethke approach.

Definition at line 803 of file GeochemicalSystem.C.

{
  return _bulk_moles_old;
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryTimeDependentReactor::newTemperature(), GeochemistryConsoleOutput::output(), GeochemistryTimeDependentReactor::preSolveFlush(), GeochemistryTimeDependentReactor::removeCurrentSpecies(), TEST_F(), TEST_F(), TEST_F(), and TEST_F().

◆ getBulkOldInOriginalBasis()

DenseVector< Real > GeochemicalSystem::getBulkOldInOriginalBasis ( ) const

Returns: the number of bulk-composition moles in the original basis. Note this will typically be the "old" number of bulk moles (from the previous time-step) unless addtoBulkMoles or updateOldWithCurrent or other similar methods have just been called. Note that this contains contributions from kinetic species, in contrast to the Bethke approach.

Definition at line 809 of file GeochemicalSystem.C.

{
  DenseVector<Real> result(_bulk_moles_old);
  if (_mgd.swap_to_original_basis.n() == 0)
    return result; // no swaps have been performed
  _mgd.swap_to_original_basis.vector_mult_transpose(result, _bulk_moles_old);
  return result;
}

Referenced by GeochemistryConsoleOutput::output().

◆ getChargeBalanceBasisIndex()

unsigned GeochemicalSystem::getChargeBalanceBasisIndex ( ) const

return the index of the charge-balance species in the basis list

Definition at line 489 of file GeochemicalSystem.C.

{
  return _charge_balance_basis_index;
}

Referenced by GeochemistryConsoleOutput::output(), GeochemicalSolver::solveSystem(), and GeochemicalSolver::swapNeeded().

◆ getConstraintMeaning()

const std::vector< GeochemicalSystem::ConstraintMeaningEnum > & GeochemicalSystem::getConstraintMeaning ( ) const

Returns: a vector of the constraint meanings for the basis species

Definition at line 2094 of file GeochemicalSystem.C.

{
  return _constraint_meaning;
}

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

◆ getEquilibriumActivity()

Real GeochemicalSystem::getEquilibriumActivity ( unsigned eqm_ind ) const

Returns the value of activity for the equilibrium species with index eqm_index.

Parameters

eqm_index the index in mgd for the equilibrium species of interest

Definition at line 2319 of file GeochemicalSystem.C.

{
  if (eqm_ind >= _num_eqm)
    mooseError("Cannot retrieve activity for equilibrium species ",
               eqm_ind,
               " since there are only ",
               _num_eqm,
               " equilibrium species");
  if (_mgd.eqm_species_mineral[eqm_ind])
    return 1.0;
  else if (_mgd.eqm_species_gas[eqm_ind])
  {
    Real log10f = 0.0;
    for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
      log10f += _mgd.eqm_stoichiometry(eqm_ind, basis_i) * std::log10(_basis_activity[basis_i]);
    log10f -= getLog10K(eqm_ind);
    return std::pow(10.0, log10f);
  }
  else
    return _eqm_activity_coef[eqm_ind] * _eqm_molality[eqm_ind];
}

Referenced by addKineticRates(), computeAndGetEquilibriumActivity(), GeochemistryQuantityAux::computeValue(), and TEST_F().

◆ getEquilibriumActivityCoefficient() [1/2]

const std::vector< Real > & GeochemicalSystem::getEquilibriumActivityCoefficient ( ) const

Returns: the activity coefficients for the equilibrium species. These are not defined for minerals

Definition at line 962 of file GeochemicalSystem.C.

{
  return _eqm_activity_coef;
}

◆ getEquilibriumActivityCoefficient() [2/2]

Real GeochemicalSystem::getEquilibriumActivityCoefficient ( unsigned j ) const

Returns: the activity coefficient for the j^th equilibrium species. This is not defined for minerals

Definition at line 950 of file GeochemicalSystem.C.

{
  if (j >= _num_eqm)
    mooseError("Cannot retrieve activity coefficient for equilibrium species ",
               j,
               " since there are only ",
               _num_eqm,
               " equilibrium species");
  return _eqm_activity_coef[j];
}

Referenced by GeochemistryConsoleOutput::output(), and TEST_F().

◆ getEquilibriumMolality() [1/2]

const std::vector< Real > & GeochemicalSystem::getEquilibriumMolality ( ) const

Returns: the molalities of the equilibrium species. These are not defined for minerals or gases

Definition at line 900 of file GeochemicalSystem.C.

{
  return _eqm_molality;
}

◆ getEquilibriumMolality() [2/2]

Real GeochemicalSystem::getEquilibriumMolality ( unsigned j ) const

Returns: the molality of the j^th equilibrium species. This is not defined for minerals or gases

Definition at line 888 of file GeochemicalSystem.C.

{
  if (j >= _num_eqm)
    mooseError("Cannot retrieve molality for equilibrium species ",
               j,
               " since there are only ",
               _num_eqm,
               " equilibrium species");
  return _eqm_molality[j];
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryConsoleOutput::output(), GeochemistryTimeDependentReactor::preSolveDump(), GeochemicalSolver::swapNeeded(), and TEST_F().

◆ getInAlgebraicSystem()

const std::vector< bool > & GeochemicalSystem::getInAlgebraicSystem ( ) const

Returns: a vector of length _num_basis whose entries determine whether the basis species is in the algebraic system

Definition at line 665 of file GeochemicalSystem.C.

{
  return _in_algebraic_system;
}

Referenced by GeochemicalSolver::swapNeeded().

◆ getIonicStrength()

Real GeochemicalSystem::getIonicStrength ( ) const

Get the ionic strength.

Definition at line 1789 of file GeochemicalSystem.C.

{
  return _is.ionicStrength(_mgd, _basis_molality, _eqm_molality, _kin_moles);
}

Referenced by GeochemistryConsoleOutput::output().

◆ getKineticLog10K() [1/2]

const std::vector< Real > & GeochemicalSystem::getKineticLog10K ( ) const

Returns: the value of log10(equilibrium constant) for the each kinetic species

Definition at line 558 of file GeochemicalSystem.C.

{
  return _kin_log10K;
}

◆ getKineticLog10K() [2/2]

Real GeochemicalSystem::getKineticLog10K ( unsigned kin ) const

Parameters

kin	the index of the kinetic species (must be < _num_kin)

Returns: the value of log10(equilibrium constant) for the given kinetic species

Definition at line 546 of file GeochemicalSystem.C.

{
  if (kin >= _num_kin)
    mooseError("Cannot retrieve log10K for kinetic species ",
               kin,
               " since there are only ",
               _num_kin,
               " kinetic species");
  return _kin_log10K[kin];
}

Referenced by GeochemistryConsoleOutput::output().

◆ getKineticMoles() [1/2]

const std::vector< Real > & GeochemicalSystem::getKineticMoles ( ) const

Returns: the mole number for all kinetic species

Definition at line 933 of file GeochemicalSystem.C.

{
  return _kin_moles;
}

◆ getKineticMoles() [2/2]

Real GeochemicalSystem::getKineticMoles ( unsigned kin ) const

Parameters

kin	the index of the kinetic species (must be < _num_kin)

Returns: the mole number for the given kinetic species

Definition at line 906 of file GeochemicalSystem.C.

{
  if (k >= _num_kin)
    mooseError("Cannot retrieve moles for kinetic species ",
               k,
               " since there are only ",
               _num_kin,
               " kinetic species");
  return _kin_moles[k];
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryConsoleOutput::output(), GeochemistryTimeDependentReactor::postSolveFlowThrough(), GeochemistryTimeDependentReactor::preSolveFlush(), and GeochemistryTimeDependentReactor::removeCurrentSpecies().

◆ getLog10K()

Real GeochemicalSystem::getLog10K ( unsigned j ) const

Parameters

j	the index of the equilibrium species

Returns: the value of log10(equilibrium constant) for the j^th equilibrium species

Definition at line 495 of file GeochemicalSystem.C.

{
  if (j >= _num_eqm)
    mooseError("Cannot retrieve log10K for equilibrium species ",
               j,
               " since there are only ",
               _num_eqm,
               " equilibrium species");
  return _eqm_log10K[j];
}

Referenced by getEquilibriumActivity(), GeochemistryConsoleOutput::output(), and TEST_F().

◆ getModelGeochemicalDatabase()

const ModelGeochemicalDatabase & GeochemicalSystem::getModelGeochemicalDatabase ( ) const

Returns: reference to the underlying ModelGeochemicalDatabase

Definition at line 1571 of file GeochemicalSystem.C.

{
  return _mgd;
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryConsoleOutput::output(), GeochemistryConsoleOutput::outputNernstInfo(), GeochemicalSolver::solveSystem(), GeochemicalSolver::swapNeeded(), TEST_F(), and TEST_F().

◆ getNumBasisInAlgebraicSystem()

unsigned GeochemicalSystem::getNumBasisInAlgebraicSystem ( ) const

return the number of basis species in the algebraic system

Definition at line 653 of file GeochemicalSystem.C.

{
  return _num_basis_in_algebraic_system;
}

Referenced by GeochemicalSolver::solveSystem().

◆ getNumInAlgebraicSystem()

unsigned GeochemicalSystem::getNumInAlgebraicSystem ( ) const

Returns: the number in the algebraic system (number of basis species in algebraic system + number of surface potentials + number of kinetic species)

Definition at line 647 of file GeochemicalSystem.C.

{
  return _num_in_algebraic_system;
}

Referenced by GeochemicalSolver::solveSystem().

◆ getNumInBasis()

unsigned GeochemicalSystem::getNumInBasis ( ) const

returns the number of species in the basis

Definition at line 1693 of file GeochemicalSystem.C.

{
  return _num_basis;
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryConsoleOutput::output(), and TEST_F().

◆ getNumInEquilibrium()

unsigned GeochemicalSystem::getNumInEquilibrium ( ) const

returns the number of species in equilibrium with the basis components

Definition at line 1699 of file GeochemicalSystem.C.

{
  return _num_eqm;
}

Referenced by GeochemistryConsoleOutput::output(), and TEST_F().

◆ getNumKinetic()

unsigned GeochemicalSystem::getNumKinetic ( ) const

returns the number of kinetic species

Definition at line 540 of file GeochemicalSystem.C.

{
  return _num_kin;
}

Referenced by GeochemistryConsoleOutput::output().

◆ getNumRedox()

unsigned GeochemicalSystem::getNumRedox ( ) const

Returns: the number of redox species in disequibrium

Definition at line 507 of file GeochemicalSystem.C.

{
  return _num_redox;
}

◆ getNumSurfacePotentials()

unsigned GeochemicalSystem::getNumSurfacePotentials ( ) const

return the number of surface potentials

Definition at line 659 of file GeochemicalSystem.C.

{
  return _num_surface_pot;
}

Referenced by GeochemistryConsoleOutput::output().

◆ getOriginalRedoxLHS()

const std::string & GeochemicalSystem::getOriginalRedoxLHS ( ) const

Returns: the original redox left-hand-side after the initial basis swaps

Definition at line 2265 of file GeochemicalSystem.C.

{
  return _original_redox_lhs;
}

Referenced by GeochemistryConsoleOutput::outputNernstInfo().

◆ getRedoxLog10K()

Real GeochemicalSystem::getRedoxLog10K ( unsigned red ) const

Parameters

red	the index of the redox species in disequilibrium (red < _num_redox)

Returns: the value of log10(equilibrium constant) for the given disequilibrium-redox species

Definition at line 513 of file GeochemicalSystem.C.

{
  if (red >= _num_redox)
    mooseError("Cannot retrieve log10K for redox species ",
               red,
               " since there are only ",
               _num_redox,
               " redox species");
  return _redox_log10K[red];
}

Referenced by GeochemistryConsoleOutput::output(), and GeochemistryConsoleOutput::outputNernstInfo().

◆ getResidualComponent()

Real GeochemicalSystem::getResidualComponent	(	unsigned	algebraic_ind,
		const DenseVector< Real > &	mole_additions
	)		const

Return the residual of the algebraic system for the given algebraic index.

Parameters

algebraic_ind the algebraic index

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.

Definition at line 1178 of file GeochemicalSystem.C.

{
  if (algebraic_ind >= _num_in_algebraic_system)
    mooseError("Cannot retrieve residual for algebraic index ",
               algebraic_ind,
               " because there are only ",
               _num_basis_in_algebraic_system,
               " molalities in the algebraic system and ",
               _num_surface_pot,
               " surface potentials and ",
               _num_kin,
               " kinetic species");
  if (mole_additions.size() != _num_basis + _num_kin)
    mooseError("The increment in mole numbers (mole_additions) needs to be of size ",
               _num_basis,
               " + ",
               _num_kin,
               " but it is of size ",
               mole_additions.size());
 
  if (algebraic_ind < _num_basis_in_algebraic_system) // residual for basis molality
  {
    /*
     * Without the special things for water or the charge-balance species, we're trying to solve
     * bulk_new = bulk_old + mole_addition
     * where bulk_old is known (from previous time-step or, for a steady-state problem, the
     * constraint)
     * and mole_addition is given to this function
     * and bulk_new = nw * (m + sum_eqm[stoi * mol_eqm]) + sum_kin[stoi * mole_kin]
     * Hence, the residual is
     * R = -(bulk_old + mole_addition) + nw*(m + sum_eqm[stoi*mol_eqm]) + sum_kin[stoi*mole_kin]
     * This is seemingly different from Bethke Eqns(16.7)-(16.9), because Bethke bulk ("M") do not
     contain the sum_kin[stoi * mol_kin] terms.  Converting the above residual to Bethke's
     convention:
     * R = -((nw*(m + sum_eqm[stoi*mol_eqm]) + sum_kin[stoi*mole_kin])_old + mole_addition) + nw*(m
     + sum_eqm[stoi*mol_eqm]) + sum_kin[stoi*mole_kin]
     *   = -(M_old + sum_kin[stoi*mole_kin_old] + mole_addition) + M + sum_kin[stoi*mole_kin]
     *   = -(M_old + mole_addition) + M + sum_kin[stoi*kin_mole_addition]
     * which is exactly Eqns(16.7)-(16.9).
     * One problem with Bethke's approach is that if kin_mole_addition depends on the mole number
     of the kinetic species, there is a "hidden" variable (moles of kinetic species) which is kept
     constant during the solve.  Here, including the moles of kinetic species as additional
     variables allows greater accuracy.
     */
    const unsigned basis_i = _basis_index[algebraic_ind];
    Real res = 0.0;
    if (basis_i == 0)
      res += -(_bulk_moles_old[basis_i] + mole_additions(basis_i)) +
             _basis_molality[0] * GeochemistryConstants::MOLES_PER_KG_WATER;
    else if (basis_i == _charge_balance_basis_index)
    {
      res += _basis_molality[0] * _basis_molality[basis_i];
      for (unsigned i = 0; i < _num_basis; ++i)
      {
        if (i == _charge_balance_basis_index)
          continue;
        else if (_mgd.basis_species_charge[i] ==
                 0.0) // certainly includes water, minerals and gases
          continue;
        else if (_constraint_meaning[i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
        {
          // the molalities might not yet have converged so that the bulk moles of this species
          // currently equals _bulk_moles_old + mole_additions, but we know
          // that when the solve has converged it'll have to, so:
          res += _mgd.basis_species_charge[i] * (_bulk_moles_old[i] + mole_additions(i)) /
                 _mgd.basis_species_charge[basis_i];
        }
        else
        {
          // do not know the bulk moles of this species: use the current value for molality and
          // kin_moles.  Note, there is no mole_additions here since physically any mole_additions
          // get instantly get soaked up by the fixed activity (or fixed molality, etc), and
          // mathematically when we eventually come to add mole_additions to the bulk_moles_old
          // (via addtoBulkMoles, for instance) we immediately return without adding stuff.  End
          // result: charge-neutrality should not depend on mole_additions for fixed-activity
          // (molality, etc) species.
          res += _mgd.basis_species_charge[i] * computeBulkFromMolalities(i) /
                 _mgd.basis_species_charge[basis_i];
        }
      }
    }
    else
      res += -(_bulk_moles_old[basis_i] + mole_additions(basis_i)) +
             _basis_molality[0] * _basis_molality[basis_i];
    for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
      res += _basis_molality[0] * _mgd.eqm_stoichiometry(eqm_j, basis_i) * _eqm_molality[eqm_j];
    for (unsigned kin = 0; kin < _num_kin; ++kin)
      res += _mgd.kin_stoichiometry(kin, basis_i) * _kin_moles[kin];
    return res;
  }
  else if (algebraic_ind <
           _num_basis_in_algebraic_system + _num_surface_pot) // residual for surface potential
  {
    const unsigned sp = algebraic_ind - _num_basis_in_algebraic_system;
    Real res = surfacePotPrefactor(sp) * (_surface_pot_expr[sp] - 1.0 / _surface_pot_expr[sp]);
    for (unsigned j = 0; j < _num_eqm; ++j)
      if (_mgd.surface_sorption_related[j] && _mgd.surface_sorption_number[j] == sp)
        res += _basis_molality[0] * _mgd.eqm_species_charge[j] * _eqm_molality[j];
    return res;
  }
 
  // else: residual for kinetic
  const unsigned kin = algebraic_ind - _num_basis_in_algebraic_system - _num_surface_pot;
  Real res = _kin_moles[kin] - (_kin_moles_old[kin] + mole_additions(_num_basis + kin));
  return res;
}

Referenced by GeochemicalSolver::computeResidual().

◆ getSaturationIndices()

std::vector< Real > GeochemicalSystem::getSaturationIndices ( ) const

Returns: vector v, where v[j] = saturation index of equilibrium species j = log10(activity product) - log10K. If the equilibrium species is not a mineral then v[j] = 0

Definition at line 1598 of file GeochemicalSystem.C.

{
  std::vector<Real> si(_num_eqm);
  for (unsigned j = 0; j < _num_eqm; ++j)
    if (_mgd.eqm_species_mineral[j])
      si[j] = log10ActivityProduct(j) - _eqm_log10K[j];
    else
      si[j] = 0.0;
  return si;
}

Referenced by GeochemistryConsoleOutput::output(), GeochemicalSolver::swapNeeded(), and TEST_F().

◆ getSolventMassAndFreeMolalityAndMineralMoles()

const std::vector< Real > & GeochemicalSystem::getSolventMassAndFreeMolalityAndMineralMoles ( ) const

Returns: vector v, where v[i] = mass of solvent water (i=0), or v[i] = molality of the basis aqueous species i, or v[i] = free moles of basis mineral i, whichever is appropriate

Definition at line 831 of file GeochemicalSystem.C.

{
  return _basis_molality;
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryConsoleOutput::output(), GeochemistryTimeDependentReactor::postSolveFlowThrough(), GeochemistryTimeDependentReactor::preSolveDump(), GeochemistryTimeDependentReactor::preSolveFlush(), GeochemicalSolver::swapNeeded(), TEST_F(), and TEST_F().

◆ getSolventWaterMass()

Real GeochemicalSystem::getSolventWaterMass ( ) const

Returns the mass of solvent water.

Definition at line 797 of file GeochemicalSystem.C.

{
  return _basis_molality[0];
}

◆ getSorbingSurfaceArea()

const std::vector< Real > & GeochemicalSystem::getSorbingSurfaceArea ( ) const

Returns: the sorbing surface area (units: m^2) for each sorbing surface

Definition at line 1858 of file GeochemicalSystem.C.

{
  return _sorbing_surface_area;
}

Referenced by GeochemistryConsoleOutput::output().

◆ getStoichiometricIonicStrength()

Real GeochemicalSystem::getStoichiometricIonicStrength ( ) const

Get the stoichiometric ionic strength.

Definition at line 1795 of file GeochemicalSystem.C.

{
  return _is.stoichiometricIonicStrength(_mgd, _basis_molality, _eqm_molality, _kin_moles);
}

Referenced by GeochemistryConsoleOutput::output().

◆ getSurfaceCharge()

Real GeochemicalSystem::getSurfaceCharge ( unsigned sp ) const

Parameters

sp	surface number of interest. sp < _num_surface_pot

Returns: the specific charge of the surface (units: Coulombs/m^2) for the surface number

Definition at line 1827 of file GeochemicalSystem.C.

{
  if (sp >= _num_surface_pot)
    mooseError("Cannot retrieve the surface charge for surface ",
               sp,
               " since there are only ",
               _num_surface_pot,
               " surfaces involved in surface complexation");
  // pre = mol of charge per square metre.  To convert to Coulombs/m^2:
  const Real pre = surfacePotPrefactor(sp) / _sorbing_surface_area[sp] *
                   (-_surface_pot_expr[sp] + 1.0 / _surface_pot_expr[sp]);
  return pre * GeochemistryConstants::FARADAY;
}

Referenced by GeochemistryQuantityAux::computeValue(), and GeochemistryConsoleOutput::output().

◆ getSurfacePotential()

Real GeochemicalSystem::getSurfacePotential ( unsigned sp ) const

Parameters

sp	surface number of interest. sp < _num_surface_pot

Returns: the surface potential (units Volts) for the surface number

Definition at line 1813 of file GeochemicalSystem.C.

{
  if (sp >= _num_surface_pot)
    mooseError("Cannot retrieve the surface potential for surface ",
               sp,
               " since there are only ",
               _num_surface_pot,
               " surfaces involved in surface complexation");
  return -2.0 * GeochemistryConstants::GAS_CONSTANT *
         (_temperature + GeochemistryConstants::CELSIUS_TO_KELVIN) /
         GeochemistryConstants::FARADAY * std::log(_surface_pot_expr[sp]);
}

Referenced by GeochemistryQuantityAux::computeValue(), and GeochemistryConsoleOutput::output().

◆ getSwapper()

const GeochemistrySpeciesSwapper & GeochemicalSystem::getSwapper ( ) const

Returns: a reference to the swapper used by this object

Definition at line 2277 of file GeochemicalSystem.C.

{
  return _swapper;
}

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

◆ getTemperature()

Real GeochemicalSystem::getTemperature ( ) const

Returns: the temperature in degC

Definition at line 1864 of file GeochemicalSystem.C.

{
  return _temperature;
}

Referenced by GeochemistryQuantityAux::computeValue(), GeochemistryConsoleOutput::output(), GeochemistryConsoleOutput::outputNernstInfo(), and TEST_F().

◆ getTotalChargeOld()

Real GeochemicalSystem::getTotalChargeOld ( ) const

Returns: the total charge in the system. Note this is based on the old values of bulk moles and kinetic moles (ie, from the previous time-step) so to get the current values methods like addToBulkMoles should be called, or updateOldWithCurrent

Definition at line 1059 of file GeochemicalSystem.C.

{
  Real tot_charge = 0.0;
  for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
    tot_charge += _mgd.basis_species_charge[basis_i] * _bulk_moles_old[basis_i];
  // kinetic species already counted in bulk_moles
  return tot_charge;
}

Referenced by enforceChargeBalance(), GeochemistryConsoleOutput::output(), and TEST_F().

◆ getTransportedBulkInOriginalBasis()

DenseVector< Real > GeochemicalSystem::getTransportedBulkInOriginalBasis ( ) const

Returns: the number of bulk-composition moles that are transported in reactive-transport, expressed in the original basis. Note this is computed using the existing molalities, so the result might be junk if the system is inconsistent, but will be OK if, for instance, a solve has just converged. Also note that this does not include contributions from kinetic species

Definition at line 819 of file GeochemicalSystem.C.

{
  std::vector<Real> trans_bulk;
  computeTransportedBulkFromMolalities(trans_bulk);
  DenseVector<Real> result(trans_bulk);
  if (_mgd.swap_to_original_basis.n() == 0)
    return result; // no swaps have been performed
  _mgd.swap_to_original_basis.vector_mult_transpose(result, trans_bulk);
  return result;
}

Referenced by GeochemistryQuantityAux::computeValue(), and GeochemistryConsoleOutput::output().

◆ initBulkAndFree()

void GeochemicalSystem::initBulkAndFree	(	std::vector< Real > &	bulk_moles_old,
		std::vector< Real > &	basis_molality
	)		const

private

based on _constrained_value and _constraint_meaning, populate nw, bulk_moles_old and basis_molality with reasonable initial conditions that may be used during the Newton solve of the algebraic system

Parameters

bulk_moles_old	bulk composition number of moles of the basis species
basis_molality	zeroth component is mass of solvent water, other components are either molality of the basis aqueous species or number of moles of basis mineral, whichever is appropriate

Definition at line 718 of file GeochemicalSystem.C.

{
  for (unsigned i = 0; i < _num_basis; ++i)
  {
    // water is done first, so dividing by basis_molality[0] is OK
    const Real value = _constraint_value[i];
    const ConstraintMeaningEnum meaning = _constraint_meaning[i];
    switch (meaning)
    {
      case ConstraintMeaningEnum::MOLES_BULK_WATER:
      {
        bulk_moles_old[i] = value;
        basis_molality[i] = std::max(
            _min_initial_molality,
            0.999 * value /
                GeochemistryConstants::MOLES_PER_KG_WATER); // mass of solvent water (water is an
                                                            // algebraic variable). Guess used to
                                                            // initialize the Newton process
        break;
      }
      case ConstraintMeaningEnum::KG_SOLVENT_WATER:
      {
        bulk_moles_old[i] = value * GeochemistryConstants::MOLES_PER_KG_WATER /
                            0.999; // initial guess (water is not an algebraic variable).  Will be
                                   // determined exactly during the solve
        basis_molality[i] = value; // mass of solvent water
        break;
      }
      case ConstraintMeaningEnum::MOLES_BULK_SPECIES:
      {
        bulk_moles_old[i] = value;
        basis_molality[i] = std::max(
            _min_initial_molality,
            0.9 * value / basis_molality[0]); // initial guess (i is an algebraic variable).  This
                                              // is what we solve for in the Newton process
        break;
      }
      case ConstraintMeaningEnum::FREE_MOLALITY:
      {
        bulk_moles_old[i] =
            value * basis_molality[0] / 0.9; // initial guess (i is not an algebraic variable).
                                             // Will be determined exactly during the solve
        basis_molality[i] = value;
        break;
      }
      case ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES:
      {
        bulk_moles_old[i] = value / 0.9; // initial guess (i is not an algebraic variable).  Will
                                         // be determined exactly during the solve
        basis_molality[i] = value;       // note, this is *moles*, not molality
        break;
      }
      case ConstraintMeaningEnum::FUGACITY:
      {
        bulk_moles_old[i] = 0.0; // initial guess (i is not an algebraic variable).  will be
                                 // determined exactly after the solve
        basis_molality[i] =
            0.0; // never used in any solve process, but since this is a gas should be zero, and
                 // setting this explicitly elimiates if(species=gas) checks in various loops
        break;
      }
      case ConstraintMeaningEnum::ACTIVITY:
      {
        bulk_moles_old[i] = value / 0.9; // initial guess (i is not an algebraic variable).  Will
                                         // be determined exactly during the solve
        if (i == 0)
          basis_molality[i] = 1.0; // assumption
        else
          basis_molality[i] =
              value /
              0.9; // assume activity_coefficient = 0.9.  this will be updated during the solve
        break;
      }
    }
  }
}

Referenced by initialize().

◆ initialize()

void GeochemicalSystem::initialize ( )

initialize all variables, ready for a Newton solve of the algebraic system

Definition at line 445 of file GeochemicalSystem.C.

{
  buildTemperatureDependentQuantities(_temperature);
  enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
  buildAlgebraicInfo(_in_algebraic_system,
                     _num_basis_in_algebraic_system,
                     _num_in_algebraic_system,
                     _algebraic_index,
                     _basis_index);
  initBulkAndFree(_bulk_moles_old, _basis_molality);
  buildKnownBasisActivities(_basis_activity_known, _basis_activity);
 
  _eqm_molality.assign(_num_eqm, 0.0);
  _surface_pot_expr.assign(_num_surface_pot, 1.0);
 
  computeConsistentConfiguration();
}

Referenced by checkAndInitialize().

◆ log10ActivityProduct()

Real GeochemicalSystem::log10ActivityProduct ( unsigned eqm_j ) const

private

Returns: log10(activity product) for equilibrium species eqm_j

Definition at line 998 of file GeochemicalSystem.C.

{
  Real log10ap = 0.0;
  for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
    log10ap += _mgd.eqm_stoichiometry(eqm_j, basis_i) * std::log10(_basis_activity[basis_i]);
  return log10ap;
}

Referenced by computeEqmMolalities(), and getSaturationIndices().

◆ log10KineticActivityProduct()

Real GeochemicalSystem::log10KineticActivityProduct ( unsigned kin ) const

Parameters

kin	the index of the kinetic species (must be < _num_kin)

Returns: the value of log10(activity product) for the given kinetic species

Definition at line 564 of file GeochemicalSystem.C.

{
  if (kin >= _num_kin)
    mooseError("Cannot retrieve activity product for kinetic species ",
               kin,
               " since there are only ",
               _num_kin,
               " kinetic species");
  Real log10ap = 0.0;
  for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
    log10ap += _mgd.kin_stoichiometry(kin, basis_i) * std::log10(_basis_activity[basis_i]);
  return log10ap;
}

Referenced by addKineticRates(), and GeochemistryConsoleOutput::output().

◆ log10RedoxActivityProduct()

Real GeochemicalSystem::log10RedoxActivityProduct ( unsigned red ) const

Parameters

red	the index of the redox species in disequilibrium (red < _num_redox)

Returns: the value of log10(activity product) for the given disequilibrium-redox species

Definition at line 525 of file GeochemicalSystem.C.

{
  if (red >= _num_redox)
    mooseError("Cannot retrieve activity product for redox species ",
               red,
               " since there are only ",
               _num_redox,
               " redox species");
  Real log10ap = 0.0;
  for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
    log10ap += _mgd.redox_stoichiometry(red, basis_i) * std::log10(_basis_activity[basis_i]);
  return log10ap;
}

Referenced by GeochemistryConsoleOutput::output(), and GeochemistryConsoleOutput::outputNernstInfo().

◆ operator=()

GeochemicalSystem & GeochemicalSystem::operator= ( const GeochemicalSystem & src )

inline

Copy assignment operator.

Almost all of the following is trivial. The most important non-trivial feature is copying src._mgd into our _mgd. Note this method gets called when dest = src and dest is already constructed. (Code such as GeochemicalSystem dest = src uses the copy constructor which simply sets the _mgd reference in dest equal to the _mgd reference in src, and does not make a copy of the data within src._mgd)

Definition at line 157 of file GeochemicalSystem.h.

  {
    if (this == &src) // trivial a=a situation
      return *this;
    // check for bad assignment situations.  Other "const" things that
    // needn't be the same (but probably actually are) include: _swap_out and _swap_in
    if (_num_basis != src._num_basis || _num_eqm != src._num_eqm || _num_redox != src._num_redox ||
        _num_surface_pot != src._num_surface_pot || _num_kin != src._num_kin ||
        _original_charge_balance_species != src._original_charge_balance_species)
      mooseError("GeochemicalSystem: copy assigment operator called with inconsistent fundamental "
                 "properties");
    // actually do the copying
    _mgd = src._mgd;
    _swapper = src._swapper;
    _gac = src._gac;
    _is = src._is;
    _charge_balance_species = src._charge_balance_species;
    _charge_balance_basis_index = src._charge_balance_basis_index;
    _constrained_species = src._constrained_species;
    _constraint_value = src._constraint_value;
    _original_constraint_value = src._original_constraint_value;
    _constraint_unit = src._constraint_unit;
    _constraint_user_meaning = src._constraint_user_meaning;
    _constraint_meaning = src._constraint_meaning;
    _eqm_log10K = src._eqm_log10K;
    _redox_log10K = src._redox_log10K;
    _kin_log10K = src._kin_log10K;
    _num_basis_in_algebraic_system = src._num_basis_in_algebraic_system;
    _num_in_algebraic_system = src._num_in_algebraic_system;
    _in_algebraic_system = src._in_algebraic_system;
    _algebraic_index = src._algebraic_index;
    _basis_index = src._basis_index;
    _bulk_moles_old = src._bulk_moles_old;
    _basis_molality = src._basis_molality;
    _basis_activity_known = src._basis_activity_known;
    _basis_activity = src._basis_activity;
    _eqm_molality = src._eqm_molality;
    _basis_activity_coef = src._basis_activity_coef;
    _eqm_activity_coef = src._eqm_activity_coef;
    _eqm_activity = src._eqm_activity;
    _surface_pot_expr = src._surface_pot_expr;
    _sorbing_surface_area = src._sorbing_surface_area;
    _kin_moles = src._kin_moles;
    _kin_moles_old = src._kin_moles_old;
    _iters_to_make_consistent = src._iters_to_make_consistent;
    _temperature = src._temperature;
    _min_initial_molality = src._min_initial_molality;
    _original_redox_lhs = src._original_redox_lhs;
    return *this;
  };

◆ operator==()

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

inline

Definition at line 213 of file GeochemicalSystem.h.

  {
    return (_mgd == rhs._mgd) && (_num_basis == rhs._num_basis) && (_num_eqm == rhs._num_eqm) &&
           (_num_redox == rhs._num_redox) && (_num_surface_pot == rhs._num_surface_pot) &&
           (_num_kin == rhs._num_kin) && (_swapper == rhs._swapper) &&
           (_swap_out == rhs._swap_out) && (_swap_in == rhs._swap_in) && (_gac == rhs._gac) &&
           (_is == rhs._is) && (_charge_balance_species == rhs._charge_balance_species) &&
           (_original_charge_balance_species == rhs._original_charge_balance_species) &&
           (_charge_balance_basis_index == rhs._charge_balance_basis_index) &&
           (_constrained_species == rhs._constrained_species) &&
           (_constraint_value == rhs._constraint_value) &&
           (_original_constraint_value == rhs._original_constraint_value) &&
           (_constraint_unit == rhs._constraint_unit) &&
           (_constraint_user_meaning == rhs._constraint_user_meaning) &&
           (_constraint_meaning == rhs._constraint_meaning) && (_eqm_log10K == rhs._eqm_log10K) &&
           (_redox_log10K == rhs._redox_log10K) && (_kin_log10K == rhs._kin_log10K) &&
           (_num_basis_in_algebraic_system == rhs._num_basis_in_algebraic_system) &&
           (_num_in_algebraic_system == rhs._num_in_algebraic_system) &&
           (_in_algebraic_system == rhs._in_algebraic_system) &&
           (_algebraic_index == rhs._algebraic_index) && (_basis_index == rhs._basis_index) &&
           (_bulk_moles_old == rhs._bulk_moles_old) && (_basis_molality == rhs._basis_molality) &&
           (_basis_activity_known == rhs._basis_activity_known) &&
           (_basis_activity == rhs._basis_activity) && (_eqm_molality == rhs._eqm_molality) &&
           (_basis_activity_coef == rhs._basis_activity_coef) &&
           (_eqm_activity_coef == rhs._eqm_activity_coef) && (_eqm_activity == rhs._eqm_activity) &&
           (_surface_pot_expr == rhs._surface_pot_expr) &&
           (_sorbing_surface_area == rhs._sorbing_surface_area) && (_kin_moles == rhs._kin_moles) &&
           (_kin_moles_old == rhs._kin_moles_old) &&
           (_iters_to_make_consistent == rhs._iters_to_make_consistent) &&
           (_temperature == rhs._temperature) &&
           (_min_initial_molality == rhs._min_initial_molality) &&
           (_original_redox_lhs == rhs._original_redox_lhs);
  }

◆ performSwap()

void GeochemicalSystem::performSwap	(	unsigned	swap_out_of_basis,
		unsigned	swap_into_basis
	)

Perform the basis swap, and ensure that the resulting system is consistent.

Note that this alters the bulk_moles_old of species. If you are confident that your configuration is reasonably correct, you should therefore ensure bulk_moles_old is set to the actual bulk_moles in your system prior to calling performSwap. Alternatively, if you are unsure of the correctness, just let performSwap use bulk_moles_old to set the bulk moles of your new basis.

Parameters

swap_out_of_basis	index of basis species to remove from basis
swap_into_basis	index of equilibrium species to add to basis

Definition at line 1610 of file GeochemicalSystem.C.

{
  if (swap_out_of_basis == 0)
    mooseException("GeochemicalSystem: attempting to swap out water and replace it by ",
                   _mgd.eqm_species_name[swap_into_basis],
                   ".  This could be because the algorithm would like to "
                   "swap out the charge-balance species, in which case you should choose a "
                   "different charge-balance species");
  if (swap_out_of_basis == _charge_balance_basis_index)
    mooseException("GeochemicalSystem: attempting to swap the charge-balance species out of "
                   "the basis");
  if (_mgd.basis_species_gas[swap_out_of_basis])
    mooseException("GeochemicalSystem: attempting to swap a gas out of the basis");
  if (_mgd.eqm_species_gas[swap_into_basis])
    mooseException("GeochemicalSystem: attempting to swap a gas into the basis");
  // perform the swap
  performSwapNoCheck(swap_out_of_basis, swap_into_basis);
}

Referenced by GeochemistryTimeDependentReactor::preSolveDump(), GeochemicalSolver::solveSystem(), and TEST_F().

◆ performSwapNoCheck()

void GeochemicalSystem::performSwapNoCheck	(	unsigned	swap_out_of_basis,
		unsigned	swap_into_basis
	)

Perform the basis swap, and ensure that the resulting system is consistent.

Note: no sanity-checks regarding swapping out the charge-balance species, etc, are made. Note that this alters the bulk_moles_old of species. If you are confident that your configuration is reasonably correct, you should therefore ensure bulk_moles_old is set to the actual bulk_moles in your system prior to calling performSwap. Alternatively, if you are unsure of the correctness, just let performSwap use bulk_moles_old to set the bulk moles of your new basis.

Parameters

swap_out_of_basis	index of basis species to remove from basis
swap_into_basis	index of equilibrium species to add to basis

Definition at line 1630 of file GeochemicalSystem.C.

{
  DenseVector<Real> bm = _bulk_moles_old;
  _swapper.performSwap(_mgd, bm, swap_out_of_basis, swap_into_basis);
 
  // the swap_into_basis species now has fixed bulk moles irrespective of what the
  // swap_out_of_basis species had fixed
  _constraint_meaning[swap_out_of_basis] = ConstraintMeaningEnum::MOLES_BULK_SPECIES;
 
  // the bulk moles will have changed for many components.  The _bulk_moles_old values will be
  // set in computeConsistentConfiguration, below
  for (unsigned basis_i = 0; basis_i < _num_basis; ++basis_i)
    if (_constraint_meaning[basis_i] == ConstraintMeaningEnum::MOLES_BULK_SPECIES ||
        _constraint_meaning[basis_i] == ConstraintMeaningEnum::MOLES_BULK_WATER)
    {
      _constraint_value[basis_i] = bm(basis_i);
      _original_constraint_value[basis_i] = bm(basis_i);
    }
 
  // In the following, the molalities are carefully swapped so that the configuration stays
  // reasonably close to the configuration prior to the swap.  This may be advantageous as it
  // initializes the Newton process with a better starting guess than what might have been used
  Real molality_of_species_just_swapped_in =
      _eqm_molality[swap_into_basis]; // is positive, or zero iff (mineral or gas)
  if (_mgd.basis_species_mineral[swap_out_of_basis] || _mgd.basis_species_gas[swap_out_of_basis] ||
      _eqm_molality[swap_into_basis] == 0.0)
  {
    // the species just swapped in is a mineral or a gas, so its equilibium molality is
    // undefined: make a guess for its molality
    molality_of_species_just_swapped_in =
        std::max(_min_initial_molality, 0.9 * bm(swap_out_of_basis));
  }
  Real molality_of_species_just_swapped_out =
      _basis_molality[swap_out_of_basis]; // can be negative if a consumed mineral
  _basis_molality[swap_out_of_basis] = molality_of_species_just_swapped_in;
  _eqm_molality[swap_into_basis] =
      std::max(0.0, molality_of_species_just_swapped_out); // this gets recalculated in
                                                           // computeConsistentConfiguration, below
 
  // depending if minerals were swapped in or out of the basis, the known activity may have
  // changed
  buildKnownBasisActivities(_basis_activity_known, _basis_activity);
 
  // the equilibrium constants will have changed due to the swap
  buildTemperatureDependentQuantities(_temperature);
 
  // due to re-orderings in mgd, the charge-balance basis index might have changed
  _charge_balance_basis_index = _mgd.basis_species_index.at(_charge_balance_species);
  // charge balance might be able to be performed easily
  enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
 
  // the algebraic system has probably changed
  buildAlgebraicInfo(_in_algebraic_system,
                     _num_basis_in_algebraic_system,
                     _num_in_algebraic_system,
                     _algebraic_index,
                     _basis_index);
 
  // finally compute a consistent configuration, based on the basis molalities, etc, above
  computeConsistentConfiguration();
}

Referenced by changeConstraintToBulk(), GeochemistryConsoleOutput::outputNernstInfo(), and performSwap().

◆ revertToOriginalChargeBalanceSpecies()

bool GeochemicalSystem::revertToOriginalChargeBalanceSpecies ( )

Changes the charge-balance species to the original that is specified in the constructor (this might not be possible since the original charge-balance species might have been swapped out)

Returns: true if the charge-balance species is changed

Definition at line 1776 of file GeochemicalSystem.C.

{
  if (_mgd.basis_species_index.find(_original_charge_balance_species) ==
      _mgd.basis_species_index.end())
    return false; // original charge-balance species no longer in basis
  const unsigned original_index = _mgd.basis_species_index.at(_original_charge_balance_species);
  if (original_index == _charge_balance_basis_index)
    return false; // current charge-balance species is the original
  setChargeBalanceSpecies(original_index);
  return true;
}

◆ setAlgebraicVariables()

void GeochemicalSystem::setAlgebraicVariables ( const DenseVector< Real > & algebraic_var )

Set the variables in the algebraic system (molalities and potentially the mass of solvent water, and surface potentials, and kinetic mole numbers if any) to algebraic_var.

All elements of this vector must be postivie. This function also calls computeConsistentConfiguration.

Parameters

algebraic_var the values to set the algebraic variables to

Definition at line 1085 of file GeochemicalSystem.C.

{
  if (algebraic_var.size() != _num_in_algebraic_system)
    mooseError("Incorrect size in setAlgebraicVariables");
  for (unsigned a = 0; a < _num_in_algebraic_system; ++a)
    if (algebraic_var(a) <= 0.0)
      mooseError("Cannot set algebraic variables to non-positive values such as ",
                 algebraic_var(a));
 
  for (unsigned a = 0; a < _num_basis_in_algebraic_system; ++a)
    _basis_molality[_basis_index[a]] = algebraic_var(a);
  for (unsigned s = 0; s < _num_surface_pot; ++s)
    _surface_pot_expr[s] = algebraic_var(s + _num_basis_in_algebraic_system);
  for (unsigned k = 0; k < _num_kin; ++k)
    _kin_moles[k] = algebraic_var(k + _num_basis_in_algebraic_system + _num_surface_pot);
 
  computeConsistentConfiguration();
}

Referenced by GeochemicalSolver::reduceInitialResidual(), and GeochemicalSolver::solveSystem().

◆ setChargeBalanceSpecies()

void GeochemicalSystem::setChargeBalanceSpecies ( unsigned new_charge_balance_index )

private

Set the charge-balance species to the basis index provided.

No checks are made on the sanity of the desired change. Before setting, the _constraint_value mole number of the old charge-balance species is set to the value provided in the constructor. Then _charge_balance_basis_index and _charge_balance_species is set appropriately

Definition at line 1705 of file GeochemicalSystem.C.

{
  // because the original mole number of the charge balance species may have been changed due to
  // enforcing charge balance:
  _constraint_value[_charge_balance_basis_index] =
      _original_constraint_value[_charge_balance_basis_index];
  _bulk_moles_old[_charge_balance_basis_index] = _constraint_value[_charge_balance_basis_index];
  // now change the charge balance info
  _charge_balance_basis_index = new_charge_balance_index;
  _charge_balance_species = _mgd.basis_species_name[_charge_balance_basis_index];
  // enforce charge-balance if easily possible
  enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
}

Referenced by alterChargeBalanceSpecies(), and revertToOriginalChargeBalanceSpecies().

◆ setConstraintValue()

void GeochemicalSystem::setConstraintValue	(	unsigned	basis_ind,
		Real	value
	)

Set the constraint value for the basis species.

If molalities or activities are changed, this method uses computeConsistentConfiguration to result in a consistent configuration, while if bulk composition is altered it uses alterSystemBecauseBulkChanged to result in a consistent configuration (charge-balance is performed if all constraints are BULK-type)

Parameters

basis_ind	the index of the basis species
value	the value to set the constraint to

Definition at line 2199 of file GeochemicalSystem.C.

{
  if (basis_ind >= _num_basis)
    mooseError("Cannot setConstraintValue for species ",
               basis_ind,
               " because there are only ",
               _num_basis,
               " basis species");
  _constraint_value[basis_ind] = value;
  _original_constraint_value[basis_ind] = value;
  switch (_constraint_meaning[basis_ind])
  {
    case ConstraintMeaningEnum::MOLES_BULK_SPECIES:
    case ConstraintMeaningEnum::MOLES_BULK_WATER:
    {
      alterSystemBecauseBulkChanged();
      break;
    }
    case ConstraintMeaningEnum::KG_SOLVENT_WATER:
    case ConstraintMeaningEnum::FREE_MOLALITY:
    case ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES:
    {
      // the changes resulting from this change are very similar to setAlgebraicVariables
      _basis_molality[basis_ind] = value;
      computeConsistentConfiguration();
      break;
    }
    case ConstraintMeaningEnum::FUGACITY:
    {
      // the changes resulting from this change are very similar to setAlgebraicVariables
      _basis_activity[basis_ind] = value;
      _basis_molality[basis_ind] = 0.0;
      computeConsistentConfiguration();
      break;
    }
    case ConstraintMeaningEnum::ACTIVITY:
    {
      // the changes resulting from this change are very similar to setAlgebraicVariables
      _basis_activity[basis_ind] = value;
      _basis_molality[basis_ind] = _basis_activity[basis_ind] / _basis_activity_coef[basis_ind];
      computeConsistentConfiguration();
      break;
    }
  }
}

Referenced by addToBulkMoles(), changeConstraintToBulk(), GeochemistryTimeDependentReactor::execute(), and TEST_F().

◆ setKineticMoles()

void GeochemicalSystem::setKineticMoles	(	unsigned	kin,
		Real	moles
	)

Sets the current AND old mole number for a kinetic species.

Note: this does not compute a consistent configuration (viz, the bulk mole composition is not updated): use computeConsistentConfiguration if you need to.

Parameters

kin	the index of the kinetic species (must be < _num_kin)
moles	the number of moles

Definition at line 918 of file GeochemicalSystem.C.

{
  if (k >= _num_kin)
    mooseError("Cannot set moles for kinetic species ",
               k,
               " since there are only ",
               _num_kin,
               " kinetic species");
  if (moles <= 0.0)
    mooseError("Mole number for kinetic species must be positive, not ", moles);
  _kin_moles[k] = moles;
  _kin_moles_old[k] = moles;
}

Referenced by checkAndInitialize(), GeochemistryTimeDependentReactor::postSolveFlowThrough(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and TEST_F().

◆ setMineralRelatedFreeMoles()

void GeochemicalSystem::setMineralRelatedFreeMoles ( Real value )

Set the free mole number of mineral-related species to the value provided.

This sets the mole number of basis minerals, the molality of unoccupied sorption sites (whether in the basis or not), and the molality of sorbed species to the value provided. It does not set the molality of equilibrium minerals, which is always zero. NOTE: calling this method can easily result in a completely inconsistent GeochemicalSystem because no computeConsistentConfiguration() is called. If you need a consistent configuration, please call that method after calling this one

Parameters

value the mole number (or molality, for sorption-related species)

Definition at line 2283 of file GeochemicalSystem.C.

{
  // mole numbers of basis minerals
  for (unsigned i = 0; i < _num_basis; ++i)
    if (_mgd.basis_species_mineral[i])
      _basis_molality[i] = value;
 
  // mole numbers of sorption sites
  for (const auto & name_info :
       _mgd.surface_complexation_info) // all minerals involved in surface complexation
    for (const auto & name_frac :
         name_info.second.sorption_sites) // all sorption sites on the given mineral
    {
      if (_mgd.basis_species_index.count(name_frac.first))
        _basis_molality[_mgd.basis_species_index.at(name_frac.first)] = value;
      else
        _eqm_molality[_mgd.eqm_species_index.at(name_frac.first)] = value;
    }
 
  // mole numbers of sorbed equilibrium species
  for (unsigned j = 0; j < _num_eqm; ++j)
  {
    if (_mgd.eqm_species_mineral[j])
      _eqm_molality[j] = 0.0; // Note: explicitly setting to zero, irrespective of user input,
                              // to ensure consistency with the rest of the code
    if (_mgd.surface_sorption_related[j])
      _eqm_molality[j] = value;
  }
 
  // mole numbers of kinetic species
  for (unsigned k = 0; k < _num_kin; ++k)
    if (_mgd.kin_species_mineral[k])
      _kin_moles[k] = value;
}

Referenced by GeochemistryTimeDependentReactor::postSolveFlowThrough(), and GeochemistryTimeDependentReactor::preSolveDump().

◆ setModelGeochemicalDatabase()

void GeochemicalSystem::setModelGeochemicalDatabase ( const ModelGeochemicalDatabase & mgd )

Copies a ModelGeochemicalDatabase into our _mgd structure.

Parameters

mgd	reference to the ModelGeochemicalDatabase that will be copied into _mgd

Definition at line 2271 of file GeochemicalSystem.C.

{
  _mgd = mgd;
}

Referenced by GeochemistryConsoleOutput::outputNernstInfo(), and TEST_F().

◆ setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles()

void GeochemicalSystem::setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles	(	const std::vector< std::string > &	names,
		const std::vector< Real > &	values,
		const std::vector< bool > &	constraints_from_molalities
	)

Sets:

solvent water mass
free molality of all basis species and equilibrium species that are not (water or gas or mineral)
free number of moles of all minerals (if any)
surface potential expressions for all surface-sorption sites (if any)
mole number and old mole number of the kinetic species (if any) Then computes bulk mole composition, activities, activity coefficients and sorbing surface areas so that the GeochemicalSystem is kept self-consistent.

This method is designed to be run at the start of a simulation, or during a "recover" operator.

Parameters

names	the names of all the basis species, equilibrium species, surface-sorption minerals and kinetic species that are in the system. Note, there must be exactly num_basis + num_eqm + num_surface_pot + num_kin of these, ie, all species must be provided with a value
values	the values of the molalities (or solvent water mass for water, or free number of moles for minerals/kinetics, or surface_potential_expr for surface-sorption sites). There must be an equal number of these as "names".
contraints_from_molalities	whether the constraints initial provided to the GeochemicalSystem should be updated by the values provided. This must have size num_basis.

This method is reasonably nontrivial:

it assumes that temperature-dependent quantities (equilibrium constants, Debye-Huckel parameters, etc) have been set prior to calling this method, eg in the constructor or setTemperature()
it assumes the algebraic info has been built during instantiation, or during swaps, etc
it does not assume the value provided are "good", although since they most usually come from a previous solve, they are typically pretty "good". It does check for non-negativity: molality for basis gases must be zero; molality for basis minerals must be non-negative; molality for every other basis species must be positive; molality for equilibrium gases and equilibrium minerals are set to zero irrespective of values provided; molality of other equilibrium species must be non-negative; surface_potential_expressions must be positive; mole number of kinetic species must be positive. (In the previous sentence, "molality" is molality, kg solvent water, or free mineral moles, whichever is appropriate.)
if constraints_from_molalities is false then: if the original constraint was KG_SOLVENT_WATER, FREE_MOLALITY or FREE_MOLES_MINERAL_SPECIES then the constraint take precedence over the molality provided to this method, so the molality is ignored and the constraint value used instead.
if constraints_from_molalities is true then: if the original constraint was KG_SOLVENT_WATER, FREE_MOLALITY or FREE_MOLES_MINERAL_SPECIES then the constraint value is set to the value provided by to this method; if the original constraint was MOLES_BULK_WATER or MOLES_BULK_SPECIES then the constraint value is set to the value computed from the molalities provided to this method; if the original constraint was ACTIVITY, then the constraint value is set to activity_coefficient * molality_provided_to_this_method
Note the possibilities for ignoring the values provided to this method mentioned in the preceeding paragraphs (setting to zero for equilibrium minerals and gases, and constraints overriding the values provided)!

Definition at line 1878 of file GeochemicalSystem.C.

{
  // assume temperature-dependent quantities have been built during instantiation
  // assume algebraic info has been built during instantiation
 
  /*
   * STEP 0: Check sizes are correct
   */
  const unsigned num_names = names.size();
  if (num_names != values.size())
    mooseError("When setting all molalities, names and values must have same size");
  if (num_names != _num_basis + _num_eqm + _num_surface_pot + _num_kin)
    mooseError("When setting all molalities, values must be provided for every species and surface "
               "potentials");
  if (constraints_from_molalities.size() != _num_basis)
    mooseError("constraints_from_molalities has size ",
               constraints_from_molalities.size(),
               " while number of basis species is ",
               _num_basis);
 
  /*
   * STEP 1: Read values into _basis_molality, _eqm_molality and _surface_pot_expr and
   * _kin_moles
   */
  // The is no guarantee is made that the user-supplied molalities are "good", except they must
  // be non-negative. There are lots of "finds" in the loops below, which is slow, but it
  // ensures all species are given molalities.  This method is designed to be called only every
  // once-in-a-while (eg, at the start of a simulation)
  for (const auto & name_index : _mgd.basis_species_index)
  {
    const unsigned ind =
        std::distance(names.begin(), std::find(names.begin(), names.end(), name_index.first));
    if (ind < num_names)
    {
      if (_mgd.basis_species_gas[name_index.second])
      {
        if (values[ind] != 0.0)
          mooseError("Molality for gas ",
                     name_index.first,
                     " cannot be ",
                     values[ind],
                     ": it must be zero");
        else
          _basis_molality[name_index.second] = values[ind];
      }
      else if (_mgd.basis_species_mineral[name_index.second])
      {
        if (values[ind] < 0.0)
          mooseError("Molality for mineral ",
                     name_index.first,
                     " cannot be ",
                     values[ind],
                     ": it must be non-negative");
        else
          _basis_molality[name_index.second] = values[ind];
      }
      else if (values[ind] <= 0.0)
        mooseError("Molality for species ",
                   name_index.first,
                   " cannot be ",
                   values[ind],
                   ": it must be positive");
      else
        _basis_molality[name_index.second] = values[ind];
    }
    else
      mooseError("Molality (or free mineral moles, etc - whatever is appropriate) for species ",
                 name_index.first,
                 " needs to be provided when setting all molalities");
  }
  for (const auto & name_index : _mgd.eqm_species_index)
  {
    const unsigned ind =
        std::distance(names.begin(), std::find(names.begin(), names.end(), name_index.first));
    if (ind < num_names)
    {
      if (_mgd.eqm_species_gas[name_index.second] || _mgd.eqm_species_mineral[name_index.second])
        _eqm_molality[name_index.second] =
            0.0; // Note: explicitly setting to zero, irrespective of user input.  The reason
                 // for doing this is that during a restore, a basis species with positive
                 // molality could become a secondary species, which should have zero molality
      else if (values[ind] < 0.0)
        mooseError("Molality for species ",
                   name_index.first,
                   " cannot be ",
                   values[ind],
                   ": it must be non-negative");
      else
        _eqm_molality[name_index.second] = values[ind];
    }
    else
      mooseError("Molality for species ",
                 name_index.first,
                 " needs to be provided when setting all molalities");
  }
  for (unsigned sp = 0; sp < _num_surface_pot; ++sp)
  {
    const unsigned ind =
        std::distance(names.begin(),
                      std::find(names.begin(),
                                names.end(),
                                _mgd.surface_sorption_name[sp] + "_surface_potential_expr"));
    if (ind < num_names)
    {
      if (values[ind] <= 0.0)
        mooseError("Surface-potential expression for mineral ",
                   _mgd.surface_sorption_name[sp],
                   " cannot be ",
                   values[ind],
                   ": it must be positive");
      _surface_pot_expr[sp] = values[ind];
    }
    else
      mooseError("Surface potential for mineral ",
                 _mgd.surface_sorption_name[sp],
                 " needs to be provided when setting all molalities");
  }
  for (const auto & name_index : _mgd.kin_species_index)
  {
    const unsigned ind =
        std::distance(names.begin(), std::find(names.begin(), names.end(), name_index.first));
    if (ind < num_names)
      setKineticMoles(name_index.second, values[ind]);
    else
      mooseError("Moles for species ",
                 name_index.first,
                 " needs to be provided when setting all molalities");
  }
 
  /*
   * STEP 2: Alter _constraint_values if necessary
   */
  // If some of the constraints_from_molalities are false, then the molalities provided to this
  // method may have to be modified to satisfy the contraints: this alters _basis_molality so
  // must occur first
  for (unsigned i = 0; i < _num_basis; ++i)
  {
    const ConstraintMeaningEnum meaning = _constraint_meaning[i];
    if (meaning == ConstraintMeaningEnum::KG_SOLVENT_WATER ||
        meaning == ConstraintMeaningEnum::FREE_MOLALITY ||
        meaning == ConstraintMeaningEnum::FREE_MOLES_MINERAL_SPECIES)
    {
      if (constraints_from_molalities[i])
      {
        // molalities provided to this method dictate the contraints:
        _constraint_value[i] = _basis_molality[i];
        _original_constraint_value[i] = _constraint_value[i];
      }
      else
        // contraints take precedence over the molalities provided to this method:
        _basis_molality[i] = _constraint_value[i];
    }
  }
 
  // Potentially alter _constraint_value for the BULK contraints:
  for (unsigned i = 0; i < _num_basis; ++i)
  {
    const ConstraintMeaningEnum meaning = _constraint_meaning[i];
    if (meaning == ConstraintMeaningEnum::MOLES_BULK_WATER ||
        meaning == ConstraintMeaningEnum::MOLES_BULK_SPECIES)
      if (constraints_from_molalities[i])
      {
        // the constraint value should be overridden by the molality-computed bulk mole number
        _constraint_value[i] = computeBulkFromMolalities(i);
        _original_constraint_value[i] = _constraint_value[i];
      }
  }
 
  // Potentially alter _constraint_value for ACTIVITY contraints:
  for (unsigned i = 0; i < _num_basis; ++i)
    if (constraints_from_molalities[i])
    {
      const ConstraintMeaningEnum meaning = _constraint_meaning[i];
      if (meaning == ConstraintMeaningEnum::FUGACITY)
        mooseError("Gas fugacity cannot be determined from molality: constraints_from_molalities "
                   "must be set false for all gases");
      else if (meaning == ConstraintMeaningEnum::ACTIVITY)
      {
        if (i == 0)
          mooseError("Water activity cannot be determined from molalities: "
                     "constraints_from_molalities "
                     "must be set to false for water if activity of water is fixed");
        // the constraint value should be overidden by the molality provided to this method
        _constraint_value[i] = _basis_activity_coef[i] * _basis_molality[i];
        _original_constraint_value[i] = _constraint_value[i];
      }
    }
 
  /*
   * STEP 3: Follow the initialize() and computeConsistentConfiguration() methods
   */
  // assume done already: buildTemperatureDependentQuantities
  enforceChargeBalanceIfSimple(_constraint_value, _bulk_moles_old);
  // assume done already: buildAlgebraicInfo
  // should not be done, as basis_molality is set by this method instead: initBulkAndFree
  buildKnownBasisActivities(_basis_activity_known, _basis_activity);
  // should not be done, as these are set by this method: _eqm_molality.assign
  // should not be done, as these are set by this method: _surface_pot_expr.assign
  _gac.setInternalParameters(_temperature, _mgd, _basis_molality, _eqm_molality, _kin_moles);
  _gac.buildActivityCoefficients(_mgd, _basis_activity_coef, _eqm_activity_coef);
  // for constraints_from_molalities = false then the following: (1) overrides the basis
  // molality provided to this method; (2) produces a slightly inconsistent equilibrium
  // geochemical system because basis_activity_coef was computed on the basis of the basis
  // molalities provided to this method
  updateBasisMolalityForKnownActivity(_basis_molality);
  computeRemainingBasisActivities(_basis_activity);
  // should not be done, as these are set by this method: computeEqmMolalities
  computeBulk(_bulk_moles_old);
  // should not be done, as these are set by this method: computeFreeMineralMoles
  computeSorbingSurfaceArea(_sorbing_surface_area);
}

Referenced by TEST_F(), TEST_F(), TEST_F(), and TEST_F().

◆ setTemperature()

void GeochemicalSystem::setTemperature ( Real temperature )

Sets the temperature to the given quantity.

This also:

calculates new values of equilibrium constants for equilibrium, redox and kinetic species
instructs the activity-coefficient calculator to set its internal temperature and calculate new values of the Debye-Huckel parameters. However, it does NOT update activity coefficients, basis molalities for species of known activities, basis activities, equilibrium molalities, bulk mole number of species with fixed free molality, free mineral mole numbers, or surface sorbing area. Hence, the state of the equilibrium geochemical system will be left inconsistent after a call to setTemperature, and you should computeConsistentConfiguration() to rectify this, if needed.
Parameters

temperature the temperature in degC

Definition at line 1870 of file GeochemicalSystem.C.

{
  _temperature = temperature;
  buildTemperatureDependentQuantities(_temperature);
  _gac.setInternalParameters(_temperature, _mgd, _basis_molality, _eqm_molality, _kin_moles);
}

Referenced by GeochemistryTimeDependentReactor::execute(), GeochemistryTimeDependentReactor::postSolveExchanger(), and TEST_F().

◆ surfacePotPrefactor()

Real GeochemicalSystem::surfacePotPrefactor ( unsigned sp ) const

private

Parameters

sp	Surface potential number. sp < _num_surface_pot

Returns: (1/2) * A / F * sqrt(R * T_k * eps * eps_0 * rho * I). This appears in the surface-potential equation.

Definition at line 1801 of file GeochemicalSystem.C.

{
  return 0.5 * _sorbing_surface_area[sp] / GeochemistryConstants::FARADAY *
         std::sqrt(8.0 * GeochemistryConstants::GAS_CONSTANT *
                   (_temperature + GeochemistryConstants::CELSIUS_TO_KELVIN) *
                   GeochemistryConstants::PERMITTIVITY_FREE_SPACE *
                   GeochemistryConstants::DIELECTRIC_CONSTANT_WATER *
                   GeochemistryConstants::DENSITY_WATER *
                   _is.ionicStrength(_mgd, _basis_molality, _eqm_molality, _kin_moles));
}

Referenced by computeJacobian(), getResidualComponent(), and getSurfaceCharge().

◆ surfaceSorptionModifier()

Real GeochemicalSystem::surfaceSorptionModifier ( unsigned eqm_j ) const

private

Parameters

eqm_j the index of the equilibrium species

Returns: the modifier to the mass-balance equation for equilibrium species j, which is exp(-z * psi * F / R / TK), where z is the charge of the equilibrium species j, psi is the surface potential relevant to the equilibrium species j, F the Faraday constant, R the gas constant, and TK the temperature in Kelvin. If equilibrium species j has no relevant surface potential, unity is returned. Note that exp(-z * psi * F / R / TK) = (_surface_pot_expr)^(2z)

Definition at line 1025 of file GeochemicalSystem.C.

{
  if (eqm_j >= _num_eqm)
    return 1.0;
  if (!_mgd.surface_sorption_related[eqm_j])
    return 1.0;
  return std::pow(_surface_pot_expr[_mgd.surface_sorption_number[eqm_j]],
                  2.0 * _mgd.eqm_species_charge[eqm_j]);
}

Referenced by computeEqmMolalities().

◆ updateBasisMolalityForKnownActivity()

void GeochemicalSystem::updateBasisMolalityForKnownActivity ( std::vector< Real > & basis_molality ) const

private

For basis aqueous species (not water, minerals or gases) with activity set by the user, set basis_molality = activity / activity_coefficient.

Definition at line 986 of file GeochemicalSystem.C.

{
  for (unsigned i = 1; i < _num_basis; ++i) // don't loop over water
  {
    if (_basis_activity_known[i] && !_mgd.basis_species_mineral[i] && !_mgd.basis_species_gas[i])
      basis_molality[i] =
          _basis_activity[i] / _basis_activity_coef[i]; // just those for
                                                        // which activity is provided by the user
  }
}

Referenced by computeConsistentConfiguration(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ updateOldWithCurrent()

void GeochemicalSystem::updateOldWithCurrent ( const DenseVector< Real > & mole_additions )

Usually used at the end of a solve, to provide correct initial conditions to the next time-step, this method:

sets kin_moles_old = kin_moles
updates the constraint_values and bulk_moles_old with mole_additions, for BULK-type species
ensures charge balance holds if all species are BULK-type species
computes basis_molality for the BULK-type mineral species

computes bulk_moles_old from the molalities for non-BULK-type species Note that it is possible that you would like to enforceChargeBalance() immediately after calling this method: that is fine, there will be no negative consequences if the system has been solved

Parameters

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. Only the first _num_basis of these are used.

Definition at line 2350 of file GeochemicalSystem.C.

{
  if (mole_additions.size() != _num_basis + _num_kin)
    mooseError("The increment in mole numbers (mole_additions) needs to be of size ",
               _num_basis,
               " + ",
               _num_kin,
               " but it is of size ",
               mole_additions.size());
 
  // Copy the current kin_moles to the old values
  for (unsigned k = 0; k < _num_kin; ++k)
    _kin_moles_old[k] = _kin_moles[k];
 
  // The following:
  // - just returns if constraint is not BULK type.
  // - ensures charge balance if simple
  // - sets the constraint value and bulk_moles_old (for BULK-type species)
  // - computes basis_molality for the BULK-type mineral species
  for (unsigned i = 0; i < _num_basis; ++i)
    addToBulkMoles(i, mole_additions(i));
 
  // The following:
  // - sets bulk_moles_old to the Constraint value for BULK-type species (duplicating the above
  // function)
  // - for the other species, computes bulk_moles_old from the molalities
  computeBulk(_bulk_moles_old);
 
  // It is possible that the user would also like to enforceChargeBalance() now, and that is fine
  // - there will be no negative consequences
}

Referenced by GeochemicalSolver::solveSystem(), and TEST_F().

Member Data Documentation

◆ _algebraic_index

std::vector<unsigned> GeochemicalSystem::_algebraic_index

private

_algebraic_index[i] = index in the algebraic system of the basis species i. _basis_index[_algebraic_index[i]] = i

Definition at line 890 of file GeochemicalSystem.h.

Referenced by changeConstraintToBulk(), getAlgebraicIndexOfBasisSystem(), initialize(), operator=(), operator==(), and performSwapNoCheck().

◆ _basis_activity

std::vector<Real> GeochemicalSystem::_basis_activity

private

values of activity (for water, minerals and aqueous basis species) or fugacity (for gases)

Definition at line 905 of file GeochemicalSystem.h.

Referenced by addKineticRates(), changeConstraintToBulk(), computeConsistentConfiguration(), getBasisActivity(), getBasisActivity(), getEquilibriumActivity(), initialize(), log10ActivityProduct(), log10KineticActivityProduct(), log10RedoxActivityProduct(), operator=(), operator==(), performSwapNoCheck(), setConstraintValue(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and updateBasisMolalityForKnownActivity().

◆ _basis_activity_coef

std::vector<Real> GeochemicalSystem::_basis_activity_coef

private

basis activity coefficients

Definition at line 909 of file GeochemicalSystem.h.

Referenced by computeConsistentConfiguration(), computeRemainingBasisActivities(), getBasisActivityCoefficient(), getBasisActivityCoefficient(), operator=(), operator==(), setConstraintValue(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and updateBasisMolalityForKnownActivity().

◆ _basis_activity_known

std::vector<bool> GeochemicalSystem::_basis_activity_known

private

whether basis_activity[i] is fixed by the user

Definition at line 903 of file GeochemicalSystem.h.

Referenced by addKineticRates(), changeConstraintToBulk(), computeRemainingBasisActivities(), getBasisActivityKnown(), initialize(), operator=(), operator==(), performSwapNoCheck(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and updateBasisMolalityForKnownActivity().

◆ _basis_index

std::vector<unsigned> GeochemicalSystem::_basis_index

private

_basis_index[i] = index in the basis of the algebraic system species i, for i<num_basis_in_algebraic_system. _basis_index[_algebraic_index[i]] == i

Definition at line 892 of file GeochemicalSystem.h.

Referenced by changeConstraintToBulk(), computeJacobian(), getAlgebraicBasisValues(), getAlgebraicVariableDenseValues(), getAlgebraicVariableValues(), getBasisIndexOfAlgebraicSystem(), getResidualComponent(), initialize(), operator=(), operator==(), performSwapNoCheck(), and setAlgebraicVariables().

◆ _basis_molality

std::vector<Real> GeochemicalSystem::_basis_molality

private

IMPORTANT: this is.

number of kg of solvent water as the first element
molality of basis aqueous species, or free moles of basis mineral, whichever is appropriate
always zero for gases

Definition at line 901 of file GeochemicalSystem.h.

Referenced by addKineticRates(), alterChargeBalanceSpecies(), alterSystemBecauseBulkChanged(), computeBulkFromMolalities(), computeConsistentConfiguration(), computeFreeMineralMoles(), computeJacobian(), computeRemainingBasisActivities(), computeSorbingSurfaceArea(), computeTransportedBulkFromMolalities(), getAlgebraicBasisValues(), getAlgebraicVariableDenseValues(), getAlgebraicVariableValues(), getIonicStrength(), getResidualComponent(), getSolventMassAndFreeMolalityAndMineralMoles(), getSolventWaterMass(), getStoichiometricIonicStrength(), initialize(), operator=(), operator==(), performSwapNoCheck(), setAlgebraicVariables(), setConstraintValue(), setMineralRelatedFreeMoles(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), setTemperature(), and surfacePotPrefactor().

◆ _bulk_moles_old

std::vector<Real> GeochemicalSystem::_bulk_moles_old

private

Number of bulk moles of basis species (this includes contributions from kinetic species, in contrast to Bethke)

Definition at line 894 of file GeochemicalSystem.h.

Referenced by alterSystemBecauseBulkChanged(), computeConsistentConfiguration(), computeFreeMineralMoles(), computeJacobian(), enforceChargeBalance(), getBulkMolesOld(), getBulkOldInOriginalBasis(), getResidualComponent(), getTotalChargeOld(), initialize(), operator=(), operator==(), performSwapNoCheck(), setChargeBalanceSpecies(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and updateOldWithCurrent().

◆ _charge_balance_basis_index

unsigned GeochemicalSystem::_charge_balance_basis_index

private

The index in the list of basis species corresponding to the charge-balance species. This gets altered as the simulation progresses, due to swaps, and alterChargeBalanceSpecies.

Definition at line 864 of file GeochemicalSystem.h.

Referenced by alterChargeBalanceSpecies(), checkAndInitialize(), computeJacobian(), enforceChargeBalance(), enforceChargeBalanceIfSimple(), getChargeBalanceBasisIndex(), getResidualComponent(), operator=(), operator==(), performSwap(), performSwapNoCheck(), revertToOriginalChargeBalanceSpecies(), and setChargeBalanceSpecies().

◆ _charge_balance_species

std::string GeochemicalSystem::_charge_balance_species

private

The species used to enforce charge balance.

Definition at line 860 of file GeochemicalSystem.h.

Referenced by checkAndInitialize(), operator=(), operator==(), performSwapNoCheck(), and setChargeBalanceSpecies().

◆ _constrained_species

std::vector<std::string> GeochemicalSystem::_constrained_species

private

Names of the species in that have their values fixed to _constraint_value with meanings _constraint_meaning. In the constructor, this is ordered to have the same ordering as the basis species.

Definition at line 866 of file GeochemicalSystem.h.

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

◆ _constraint_meaning

std::vector<ConstraintMeaningEnum> GeochemicalSystem::_constraint_meaning

private

The meaning of the values in _constraint_value. In the constructor, this is ordered to have the same ordering as the basis species.

Definition at line 876 of file GeochemicalSystem.h.

Referenced by addToBulkMoles(), alterChargeBalanceSpecies(), alterSystemBecauseBulkChanged(), buildAlgebraicInfo(), buildKnownBasisActivities(), changeConstraintToBulk(), checkAndInitialize(), closeSystem(), computeBulk(), computeFreeMineralMoles(), computeJacobian(), enforceChargeBalanceIfSimple(), getConstraintMeaning(), getResidualComponent(), initBulkAndFree(), operator=(), operator==(), performSwapNoCheck(), setConstraintValue(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ _constraint_unit

std::vector<GeochemistryUnitConverter::GeochemistryUnit> GeochemicalSystem::_constraint_unit

private

Units of the _constraint_value when the GeochemicalSystem is constructed. This is used during the constructor to convert the constraint_values into mole units.

Definition at line 872 of file GeochemicalSystem.h.

Referenced by checkAndInitialize(), GeochemicalSystem(), operator=(), and operator==().

◆ _constraint_user_meaning

std::vector<ConstraintUserMeaningEnum> GeochemicalSystem::_constraint_user_meaning

private

The user-defined meaning of the values in _constraint_value. In the constructor, this is ordered to have the same ordering as the basis species. During the process of unit conversion, from the user-supplied _constraint_unit, to mole-based units, _constraint_user_meaning is used to populate _constraint_meaning, and henceforth usually only _constraint_meaning is used in the code.

Definition at line 874 of file GeochemicalSystem.h.

Referenced by checkAndInitialize(), GeochemicalSystem(), operator=(), and operator==().

◆ _constraint_value

std::vector<Real> GeochemicalSystem::_constraint_value

private

Numerical values of the constraints on _constraint_species. In the constructor, this is ordered to have the same ordering as the basis species.

Definition at line 868 of file GeochemicalSystem.h.

Referenced by addToBulkMoles(), alterSystemBecauseBulkChanged(), buildKnownBasisActivities(), checkAndInitialize(), computeBulk(), computeFreeMineralMoles(), enforceChargeBalance(), initBulkAndFree(), initialize(), operator=(), operator==(), performSwapNoCheck(), setChargeBalanceSpecies(), setConstraintValue(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ _eqm_activity

std::vector<Real> GeochemicalSystem::_eqm_activity

private

equilibrium activities.

NOTE: for computational efficiency, these are not computed until computeAndGetEqmActivity is called, and they are currently only ever used for output purposes and computing kinetic rates

Definition at line 917 of file GeochemicalSystem.h.

Referenced by addKineticRates(), computeAndGetEquilibriumActivity(), operator=(), and operator==().

◆ _eqm_activity_coef

std::vector<Real> GeochemicalSystem::_eqm_activity_coef

private

equilibrium activity coefficients

Definition at line 911 of file GeochemicalSystem.h.

Referenced by computeConsistentConfiguration(), computeEqmMolalities(), getEquilibriumActivity(), getEquilibriumActivityCoefficient(), getEquilibriumActivityCoefficient(), operator=(), operator==(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ _eqm_log10K

std::vector<Real> GeochemicalSystem::_eqm_log10K

private

equilibrium constant of the equilibrium species

Definition at line 878 of file GeochemicalSystem.h.

Referenced by buildTemperatureDependentQuantities(), computeEqmMolalities(), getLog10K(), getSaturationIndices(), operator=(), and operator==().

◆ _eqm_molality

std::vector<Real> GeochemicalSystem::_eqm_molality

private

◆ _gac

GeochemistryActivityCoefficients& GeochemicalSystem::_gac

private

Object to compute the activity coefficients and activity of water.

Definition at line 856 of file GeochemicalSystem.h.

Referenced by computeConsistentConfiguration(), computeRemainingBasisActivities(), operator=(), operator==(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and setTemperature().

◆ _in_algebraic_system

std::vector<bool> GeochemicalSystem::_in_algebraic_system

private

if _in_algebraic_system[i] == true then we need to solve for the i^th basis-species molality

Definition at line 888 of file GeochemicalSystem.h.

Referenced by changeConstraintToBulk(), getInAlgebraicSystem(), initialize(), operator=(), operator==(), and performSwapNoCheck().

◆ _is

GeochemistryIonicStrength& GeochemicalSystem::_is

private

Object that provides the ionic strengths.

Definition at line 858 of file GeochemicalSystem.h.

Referenced by getIonicStrength(), getStoichiometricIonicStrength(), operator=(), operator==(), and surfacePotPrefactor().

◆ _iters_to_make_consistent

unsigned GeochemicalSystem::_iters_to_make_consistent

private

Iterations to make the initial configuration consistent.

Note that because equilibrium molality depends on activity and activity coefficients, and water activity and activity coefficients depend on molality, it is a nontrivial task to compute all these so that the algorithm starts with a consistent initial condition from which to solve the algebraic system. Usually the algorithm doesn't even attempt to make a consistent initial condition (a suitable default is iters_to_make_consistent=0), because solving the algebraic system includes so many approximations anyway.

Definition at line 940 of file GeochemicalSystem.h.

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

◆ _kin_log10K

std::vector<Real> GeochemicalSystem::_kin_log10K

private

equilibrium constant of the kinetic species

Definition at line 882 of file GeochemicalSystem.h.

Referenced by addKineticRates(), buildTemperatureDependentQuantities(), getKineticLog10K(), getKineticLog10K(), operator=(), and operator==().

◆ _kin_moles

std::vector<Real> GeochemicalSystem::_kin_moles

private

◆ _kin_moles_old

std::vector<Real> GeochemicalSystem::_kin_moles_old

private

old mole number of kinetic species

Definition at line 930 of file GeochemicalSystem.h.

Referenced by getResidualComponent(), operator=(), operator==(), setKineticMoles(), and updateOldWithCurrent().

◆ _mgd

ModelGeochemicalDatabase& GeochemicalSystem::_mgd

private

The minimal geochemical database.

Definition at line 838 of file GeochemicalSystem.h.

Referenced by addKineticRates(), alterChargeBalanceSpecies(), buildAlgebraicInfo(), buildKnownBasisActivities(), buildTemperatureDependentQuantities(), changeConstraintToBulk(), changeConstraintToBulk(), checkAndInitialize(), computeBulkFromMolalities(), computeConsistentConfiguration(), computeEqmMolalities(), computeFreeMineralMoles(), computeJacobian(), computeSorbingSurfaceArea(), computeTransportedBulkFromMolalities(), enforceChargeBalance(), enforceChargeBalanceIfSimple(), getBulkOldInOriginalBasis(), getEquilibriumActivity(), getIonicStrength(), getModelGeochemicalDatabase(), getResidualComponent(), getSaturationIndices(), getStoichiometricIonicStrength(), getTotalChargeOld(), getTransportedBulkInOriginalBasis(), log10ActivityProduct(), log10KineticActivityProduct(), log10RedoxActivityProduct(), operator=(), operator==(), performSwap(), performSwapNoCheck(), revertToOriginalChargeBalanceSpecies(), setChargeBalanceSpecies(), setMineralRelatedFreeMoles(), setModelGeochemicalDatabase(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), setTemperature(), surfacePotPrefactor(), surfaceSorptionModifier(), and updateBasisMolalityForKnownActivity().

◆ _min_initial_molality

Real GeochemicalSystem::_min_initial_molality

private

Minimum molality ever used in an initial guess.

Definition at line 944 of file GeochemicalSystem.h.

Referenced by initBulkAndFree(), operator=(), operator==(), and performSwapNoCheck().

◆ _num_basis

const unsigned GeochemicalSystem::_num_basis

private

number of basis species

Definition at line 840 of file GeochemicalSystem.h.

Referenced by addKineticRates(), addToBulkMoles(), alterChargeBalanceSpecies(), alterSystemBecauseBulkChanged(), buildAlgebraicInfo(), buildKnownBasisActivities(), changeConstraintToBulk(), changeConstraintToBulk(), checkAndInitialize(), closeSystem(), computeBulk(), computeFreeMineralMoles(), computeJacobian(), computeRemainingBasisActivities(), computeTransportedBulkFromMolalities(), enforceChargeBalanceIfSimple(), getBasisActivity(), getBasisActivityCoefficient(), getEquilibriumActivity(), getNumInBasis(), getResidualComponent(), getTotalChargeOld(), initBulkAndFree(), log10ActivityProduct(), log10KineticActivityProduct(), log10RedoxActivityProduct(), operator=(), operator==(), performSwapNoCheck(), setConstraintValue(), setMineralRelatedFreeMoles(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), updateBasisMolalityForKnownActivity(), and updateOldWithCurrent().

◆ _num_basis_in_algebraic_system

unsigned GeochemicalSystem::_num_basis_in_algebraic_system

private

number of unknown molalities in the algebraic system. Note: surface potentials (if any) are extra quantities in the algebraic system

Definition at line 884 of file GeochemicalSystem.h.

Referenced by buildAlgebraicInfo(), changeConstraintToBulk(), computeJacobian(), getAlgebraicBasisValues(), getAlgebraicVariableDenseValues(), getAlgebraicVariableValues(), getNumBasisInAlgebraicSystem(), getResidualComponent(), initialize(), operator=(), operator==(), performSwapNoCheck(), and setAlgebraicVariables().

◆ _num_eqm

const unsigned GeochemicalSystem::_num_eqm

private

◆ _num_in_algebraic_system

unsigned GeochemicalSystem::_num_in_algebraic_system

private

number of unknowns in the algebraic system (includes molalities and surface potentials)

Definition at line 886 of file GeochemicalSystem.h.

Referenced by changeConstraintToBulk(), computeJacobian(), getAlgebraicVariableDenseValues(), getNumInAlgebraicSystem(), getResidualComponent(), initialize(), operator=(), operator==(), performSwapNoCheck(), and setAlgebraicVariables().

◆ _num_kin

const unsigned GeochemicalSystem::_num_kin

private

◆ _num_redox

const unsigned GeochemicalSystem::_num_redox

private

number of redox species

Definition at line 844 of file GeochemicalSystem.h.

Referenced by buildTemperatureDependentQuantities(), getNumRedox(), getRedoxLog10K(), log10RedoxActivityProduct(), operator=(), and operator==().

◆ _num_surface_pot

const unsigned GeochemicalSystem::_num_surface_pot

private

number of surface potentials

Definition at line 846 of file GeochemicalSystem.h.

Referenced by buildAlgebraicInfo(), computeJacobian(), computeSorbingSurfaceArea(), getAlgebraicVariableDenseValues(), getAlgebraicVariableValues(), getNumSurfacePotentials(), getResidualComponent(), getSurfaceCharge(), getSurfacePotential(), initialize(), operator=(), operator==(), setAlgebraicVariables(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ _original_charge_balance_species

const std::string GeochemicalSystem::_original_charge_balance_species

private

The species used to enforce charge balance, as provided in the constructor.

Definition at line 862 of file GeochemicalSystem.h.

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

◆ _original_constraint_value

std::vector<Real> GeochemicalSystem::_original_constraint_value

private

Numerical values of the constraints on _constraint_species. In the constructor, this is ordered to have the same ordering as the basis species. Since values can change due to charge-balance, this holds the original values set by the user.

Definition at line 870 of file GeochemicalSystem.h.

Referenced by checkAndInitialize(), operator=(), operator==(), performSwapNoCheck(), setChargeBalanceSpecies(), setConstraintValue(), and setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles().

◆ _original_redox_lhs

std::string GeochemicalSystem::_original_redox_lhs

private

The left-hand-side of the redox equations in _mgd after initial swaps.

Definition at line 946 of file GeochemicalSystem.h.

Referenced by checkAndInitialize(), getOriginalRedoxLHS(), operator=(), and operator==().

◆ _redox_log10K

std::vector<Real> GeochemicalSystem::_redox_log10K

private

equilibrium constant of the redox species

Definition at line 880 of file GeochemicalSystem.h.

Referenced by buildTemperatureDependentQuantities(), getRedoxLog10K(), operator=(), and operator==().

◆ _sorbing_surface_area

std::vector<Real> GeochemicalSystem::_sorbing_surface_area

private

surface areas of the sorbing minerals

Definition at line 926 of file GeochemicalSystem.h.

Referenced by computeConsistentConfiguration(), getSorbingSurfaceArea(), getSurfaceCharge(), operator=(), operator==(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and surfacePotPrefactor().

◆ _surface_pot_expr

std::vector<Real> GeochemicalSystem::_surface_pot_expr

private

surface potential expressions.

These are not the surface potentials themselves. Instead they are exp(-surface_potential * Faraday / 2 / R / T_k). Hence _surface_pot_expr >= 0 (like molalities) and the surface-potential residual is close to linear in _surface_pot_expr if equilibrium molalities are large

Definition at line 924 of file GeochemicalSystem.h.

Referenced by computeJacobian(), getAlgebraicVariableDenseValues(), getAlgebraicVariableValues(), getResidualComponent(), getSurfaceCharge(), getSurfacePotential(), initialize(), operator=(), operator==(), setAlgebraicVariables(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), and surfaceSorptionModifier().

◆ _swap_in

const std::vector<std::string> GeochemicalSystem::_swap_in

private

Species to immediately add to the basis in favor of _swap_out.

Definition at line 854 of file GeochemicalSystem.h.

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

◆ _swap_out

const std::vector<std::string> GeochemicalSystem::_swap_out

private

Species to immediately remove from the basis in favor of _swap_in.

Definition at line 852 of file GeochemicalSystem.h.

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

◆ _swapper

GeochemistrySpeciesSwapper& GeochemicalSystem::_swapper

private

swapper that can swap species into and out of the basis

Definition at line 850 of file GeochemicalSystem.h.

Referenced by checkAndInitialize(), getSwapper(), operator=(), operator==(), and performSwapNoCheck().

◆ _temperature

Real GeochemicalSystem::_temperature

private

The temperature in degC.

Definition at line 942 of file GeochemicalSystem.h.

Referenced by addKineticRates(), computeConsistentConfiguration(), getSurfacePotential(), getTemperature(), initialize(), operator=(), operator==(), performSwapNoCheck(), setSolventMassAndFreeMolalityAndMineralMolesAndSurfacePotsAndKineticMoles(), setTemperature(), and surfacePotPrefactor().

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

geochemistry/include/utils/GeochemicalSystem.h
geochemistry/src/utils/GeochemicalSystem.C

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ ConstraintMeaningEnum

◆ ConstraintUserMeaningEnum

Constructor & Destructor Documentation

◆ GeochemicalSystem() [1/3]

◆ GeochemicalSystem() [2/3]

◆ GeochemicalSystem() [3/3]

Member Function Documentation

◆ addKineticRates()

◆ addToBulkMoles()

◆ alterChargeBalanceSpecies()

◆ alterSystemBecauseBulkChanged()

◆ buildAlgebraicInfo()

◆ buildKnownBasisActivities()

◆ buildTemperatureDependentQuantities()

◆ changeConstraintToBulk() [1/2]

◆ changeConstraintToBulk() [2/2]

◆ checkAndInitialize()

◆ closeSystem()

◆ computeAndGetEquilibriumActivity()

◆ computeBulk()

◆ computeBulkFromMolalities()

◆ computeConsistentConfiguration()

◆ computeEqmMolalities()

◆ computeFreeMineralMoles()

◆ computeJacobian()

◆ computeRemainingBasisActivities()

◆ computeSorbingSurfaceArea()

◆ computeTransportedBulkFromMolalities()

◆ enforceChargeBalance() [1/2]

◆ enforceChargeBalance() [2/2]

◆ enforceChargeBalanceIfSimple()

◆ getAlgebraicBasisValues()

◆ getAlgebraicIndexOfBasisSystem()

◆ getAlgebraicVariableDenseValues()

◆ getAlgebraicVariableValues()

◆ getBasisActivity() [1/2]

◆ getBasisActivity() [2/2]

◆ getBasisActivityCoefficient() [1/2]

◆ getBasisActivityCoefficient() [2/2]

◆ getBasisActivityKnown()

◆ getBasisIndexOfAlgebraicSystem()

◆ getBulkMolesOld()

◆ getBulkOldInOriginalBasis()

◆ getChargeBalanceBasisIndex()

◆ getConstraintMeaning()

◆ getEquilibriumActivity()

◆ getEquilibriumActivityCoefficient() [1/2]

◆ getEquilibriumActivityCoefficient() [2/2]

◆ getEquilibriumMolality() [1/2]

◆ getEquilibriumMolality() [2/2]

◆ getInAlgebraicSystem()

◆ getIonicStrength()

◆ getKineticLog10K() [1/2]

◆ getKineticLog10K() [2/2]

◆ getKineticMoles() [1/2]

◆ getKineticMoles() [2/2]

◆ getLog10K()

◆ getModelGeochemicalDatabase()

◆ getNumBasisInAlgebraicSystem()

◆ getNumInAlgebraicSystem()

◆ getNumInBasis()

◆ getNumInEquilibrium()

◆ getNumKinetic()

◆ getNumRedox()

◆ getNumSurfacePotentials()

◆ getOriginalRedoxLHS()

◆ getRedoxLog10K()

◆ getResidualComponent()

◆ getSaturationIndices()

◆ getSolventMassAndFreeMolalityAndMineralMoles()

◆ getSolventWaterMass()

◆ getSorbingSurfaceArea()

◆ getStoichiometricIonicStrength()

◆ getSurfaceCharge()