Material designed to form a std::vector<std::vector> of mass fractions from primary-species concentrations and secondary-species concentrations for an equilibrium aqueous chemistry reaction system. More...

#include <PorousFlowMassFractionAqueousEquilibriumChemistry.h>

Inheritance diagram for PorousFlowMassFractionAqueousEquilibriumChemistry:

[legend]

Public Types
typedef DerivativeMaterialPropertyNameInterface::SymbolName	SymbolName

Public Member Functions
	PorousFlowMassFractionAqueousEquilibriumChemistry (const InputParameters &parameters)

const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialProperty (const std::string &name)

const GenericMaterialProperty< U, is_ad > &	getDefaultMaterialPropertyByName (const std::string &name)

void	validateDerivativeMaterialPropertyBase (const std::string &base)

const MaterialPropertyName	derivativePropertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	derivativePropertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	derivativePropertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	derivativePropertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

GenericMaterialProperty< U, is_ad > &	declarePropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, const SymbolName &c1, unsigned int v2, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivative (const std::string &base, unsigned int v1, unsigned int v2=libMesh::invalid_uint, unsigned int v3=libMesh::invalid_uint)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< VariableName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const std::vector< SymbolName > &c)

const GenericMaterialProperty< U, is_ad > &	getMaterialPropertyDerivativeByName (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2="", const SymbolName &c3="")

void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)

void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateCoupling (const MaterialPropertyName &base, const std::vector< VariableName > &c, bool validate_aux=true)

void	validateCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

void	validateNonlinearCoupling (const MaterialPropertyName &base, const VariableName &c1="", const VariableName &c2="", const VariableName &c3="")

const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	propertyName (const MaterialPropertyName &base, const std::vector< SymbolName > &c) const

const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	propertyNameFirst (const MaterialPropertyName &base, const SymbolName &c1) const

const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	propertyNameSecond (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2) const

const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

const MaterialPropertyName	propertyNameThird (const MaterialPropertyName &base, const SymbolName &c1, const SymbolName &c2, const SymbolName &c3) const

Static Public Member Functions
static InputParameters	validParams ()

Protected Member Functions
virtual void	initQpStatefulProperties () override

virtual void	computeQpProperties () override

Real	stoichiometry (unsigned reaction_num, unsigned primary_num) const
	The stoichiometric coefficient.

virtual void	computeQpSecondaryConcentrations ()
	Compute the secondary-species concentration as defined by the chemistry Must be overridden by derived classes.

virtual void	initQpSecondaryConcentrations ()
	Initialises (at _t_step = 0) the secondary concentrations.

void	findZeroConcentration (unsigned &zero_conc_index, unsigned &zero_count) const
	Checks gamp[i] = _primary_activity_coefficients[i] * (*_primary[i])[qp].

virtual void	dQpSecondaryConcentration_dprimary (unsigned reaction_num, std::vector< Real > &dsc) const
	Computes derivative of the secondary concentration with respect to the primary concentrations Must be overridden by derived classes.

virtual Real	dQpSecondaryConcentration_dT (unsigned reaction_num) const
	Computes derivative of the secondary concentration with respect to the temperature Must be overridden by derived classes.

void	build_mass_frac (unsigned int qp)
	Builds the mass-fraction variable matrix at the quad point.

Protected Attributes
MaterialProperty< std::vector< Real > > &	_sec_conc
	Secondary concentrations at quadpoint or nodes.

MaterialProperty< std::vector< std::vector< Real > > > &	_dsec_conc_dvar
	Derivative of the secondary concentrations with respect to the porous flow variables.

const MaterialProperty< Real > &	_temperature
	Temperature.

const MaterialProperty< std::vector< Real > > &	_dtemperature_dvar
	d(temperature)/(d porflow variable)

const unsigned int	_num_primary
	Number of primary species.

const unsigned int	_aq_ph
	Aqueous phase number.

const unsigned int	_aq_i
	Index (into _mf_vars) of the first of the primary species.

const unsigned int	_num_reactions
	Number of equations in the aqueous geochemistry system.

const bool	_equilibrium_constants_as_log10
	Whether the equilibium constants are written in their log10 form, or in absolute terms.

const unsigned	_num_equilibrium_constants
	Number of equilibrium_constants provided.

std::vector< const VariableValue * >	_equilibrium_constants
	Equilibrium constants (dimensionless)

const std::vector< Real >	_primary_activity_coefficients
	Activity coefficients for the primary species (dimensionless)

const std::vector< Real >	_reactions
	Stoichiometry defining the aqeuous geochemistry equilibrium reactions.

const std::vector< Real >	_secondary_activity_coefficients
	Activity coefficients for the secondary species.

GenericMaterialProperty< std::vector< std::vector< Real > >, is_ad > &	_mass_frac
	Mass fraction matrix at quadpoint or nodes.

GenericMaterialProperty< std::vector< std::vector< RealGradient > >, is_ad > *const	_grad_mass_frac
	Gradient of the mass fraction matrix at the quad points.

MaterialProperty< std::vector< std::vector< std::vector< Real > > > > *const	_dmass_frac_dvar
	Derivative of the mass fraction matrix with respect to the porous flow variables.

const unsigned int	_num_passed_mf_vars
	Number of mass-fraction variables provided by the user This needs to be num_phases(_num_components - 1), since the mass fraction of the final component in each phase is determined as 1 - sum*{components}(mass fraction of all other components in the phase)

std::vector< unsigned int >	_mf_vars_num
	The variable number of the mass-fraction variables.

std::vector< const GenericVariableValue< is_ad > * >	_mf_vars
	The mass-fraction variables.

std::vector< const GenericVariableGradient< is_ad > * >	_grad_mf_vars
	The gradient of the mass-fraction variables.

const unsigned int	_num_phases
	Number of phases.

const unsigned int	_num_components
	Number of fluid components.

const unsigned int	_num_var
	Number of PorousFlow variables.

Private Member Functions
bool	haveMaterialProperty (const std::string &prop_name)

std::vector< VariableName >	buildVariableVector (const VariableName &c1, const VariableName &c2, const VariableName &c3)

void	validateCouplingHelper (const MaterialPropertyName &base, const std::vector< VariableName > &c, const System &system, std::vector< VariableName > &missing)

bool	isNotObjectVariable (const VariableName &name)

Private Attributes
FEProblemBase &	_dmi_fe_problem

Detailed Description

Material designed to form a std::vector<std::vector> of mass fractions from primary-species concentrations and secondary-species concentrations for an equilibrium aqueous chemistry reaction system.

Definition at line 20 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Constructor & Destructor Documentation

◆ PorousFlowMassFractionAqueousEquilibriumChemistry()

PorousFlowMassFractionAqueousEquilibriumChemistry::PorousFlowMassFractionAqueousEquilibriumChemistry ( const InputParameters & parameters )

Definition at line 66 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

  : PorousFlowMassFraction(parameters),
    _sec_conc(_nodal_material
                  ? declareProperty<std::vector<Real>>("PorousFlow_secondary_concentration_nodal")
                  : declareProperty<std::vector<Real>>("PorousFlow_secondary_concentration_qp")),
    _dsec_conc_dvar(_nodal_material ? declareProperty<std::vector<std::vector<Real>>>(
                                          "dPorousFlow_secondary_concentration_nodal_dvar")
                                    : declareProperty<std::vector<std::vector<Real>>>(
                                          "dPorousFlow_secondary_concentration_qp_dvar")),
 
    _temperature(_nodal_material ? getMaterialProperty<Real>("PorousFlow_temperature_nodal")
                                 : getMaterialProperty<Real>("PorousFlow_temperature_qp")),
    _dtemperature_dvar(
        _nodal_material
            ? getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_nodal_dvar")
            : getMaterialProperty<std::vector<Real>>("dPorousFlow_temperature_qp_dvar")),
 
    _num_primary(_num_components - 1),
    _aq_ph(_dictator.aqueousPhaseNumber()),
    _aq_i(_aq_ph * _num_primary),
    _num_reactions(getParam<unsigned>("num_reactions")),
    _equilibrium_constants_as_log10(getParam<bool>("equilibrium_constants_as_log10")),
    _num_equilibrium_constants(coupledComponents("equilibrium_constants")),
    _equilibrium_constants(_num_equilibrium_constants),
    _primary_activity_coefficients(getParam<std::vector<Real>>("primary_activity_coefficients")),
    _reactions(getParam<std::vector<Real>>("reactions")),
    _secondary_activity_coefficients(getParam<std::vector<Real>>("secondary_activity_coefficients"))
{
  if (_dictator.numPhases() < 1)
    mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of fluid phases must "
               "not be zero");
 
  // correct number of equilibrium constants
  if (_num_equilibrium_constants != _num_reactions)
    mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of "
               "equilibrium constants is ",
               _num_equilibrium_constants,
               " which must be equal to the number of reactions (",
               _num_reactions,
               ")");
 
  // correct number of activity coefficients
  if (_primary_activity_coefficients.size() != _num_primary)
    mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of primary activity "
               "coefficients is ",
               _primary_activity_coefficients.size(),
               " which must be equal to the number of primary species (",
               _num_primary,
               ")");
 
  // correct number of stoichiometry coefficients
  if (_reactions.size() != _num_reactions * _num_primary)
    mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: The number of stoichiometric "
               "coefficients specified in 'reactions' (",
               _reactions.size(),
               ") must be equal to the number of reactions (",
               _num_reactions,
               ") multiplied by the number of primary species (",
               _num_primary,
               ")");
 
  // correct number of secondary activity coefficients
  if (_secondary_activity_coefficients.size() != _num_reactions)
    mooseError(
        "PorousFlowMassFractionAqueousEquilibriumChemistry: The number of secondary activity "
        "coefficients is ",
        _secondary_activity_coefficients.size(),
        " which must be equal to the number of secondary species (",
        _num_reactions,
        ")");
 
  // correct number of reactions
  if (_num_reactions != _dictator.numAqueousEquilibrium())
    mooseError("PorousFlowMassFractionAqueousEquilibriumChemistry: You have specified the number "
               "of reactions to be ",
               _num_reactions,
               " but the Dictator knows that the number of aqueous equilibrium reactions is ",
               _dictator.numAqueousEquilibrium());
 
  for (unsigned i = 0; i < _num_equilibrium_constants; ++i)
  {
    // If equilibrium_constants are elemental AuxVariables (or constants), we want to use
    // coupledGenericValue() rather than coupledGenericDofValue()
    const bool is_nodal = isCoupled("equilibrium_constants")
                              ? getFieldVar("equilibrium_constants", i)->isNodal()
                              : false;
 
    _equilibrium_constants[i] =
        (_nodal_material && is_nodal ? &coupledDofValues("equilibrium_constants", i)
                                     : &coupledValue("equilibrium_constants", i));
  }
}

Member Function Documentation

◆ build_mass_frac()

template<bool is_ad>

void PorousFlowMassFractionTempl< is_ad >::build_mass_frac ( unsigned int qp )

protectedinherited

Builds the mass-fraction variable matrix at the quad point.

Parameters

qp	the quad point

◆ computeQpProperties()

void PorousFlowMassFractionAqueousEquilibriumChemistry::computeQpProperties ( )

overrideprotectedvirtual

Reimplemented from PorousFlowMassFractionTempl< is_ad >.

Definition at line 167 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  // size all properties correctly and populate the non-aqueous phase info
  PorousFlowMassFraction::computeQpProperties();
 
  // size the secondary concentrations
  _sec_conc[_qp].resize(_num_reactions);
  _dsec_conc_dvar[_qp].resize(_num_reactions);
  for (unsigned r = 0; r < _num_reactions; ++r)
    _dsec_conc_dvar[_qp][r].assign(_num_var, 0.0);
 
  // Compute the secondary concentrations
  if (_t_step == 0 && !_app.isRestarting())
    initQpSecondaryConcentrations();
  else
    computeQpSecondaryConcentrations();
 
  // compute _mass_frac[_qp][_aq_ph]
  _mass_frac[_qp][_aq_ph][_num_components - 1] = 1.0; // the final component is H20
  for (unsigned i = 0; i < _num_primary; ++i)
  {
    _mass_frac[_qp][_aq_ph][i] = (*_mf_vars[_aq_i + i])[_qp];
    for (unsigned r = 0; r < _num_reactions; ++r)
      _mass_frac[_qp][_aq_ph][i] += stoichiometry(r, i) * _sec_conc[_qp][r];
 
    // remove mass-fraction from the H20 component
    _mass_frac[_qp][_aq_ph][_num_components - 1] -= _mass_frac[_qp][_aq_ph][i];
  }
 
  // Compute the derivatives of the secondary concentrations
  std::vector<std::vector<Real>> dsec(_num_reactions);
  std::vector<Real> dsec_dT(_num_reactions);
  for (unsigned r = 0; r < _num_reactions; ++r)
  {
    dQpSecondaryConcentration_dprimary(r, dsec[r]);
    dsec_dT[r] = dQpSecondaryConcentration_dT(r);
  }
 
  // Compute the derivatives of the mass_frac wrt the primary concentrations
  // This is used in _dmass_frac_dvar as well as _grad_mass_frac
  std::vector<std::vector<Real>> dmf(_num_components);
  for (unsigned i = 0; i < _num_components; ++i)
    dmf[i].assign(_num_primary, 0.0);
  for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
  {
    // run through the mass fractions (except the last one) adding to their derivatives
    // The special case is:
    dmf[wrt][wrt] = 1.0;
    // The secondary-species contributions are:
    for (unsigned i = 0; i < _num_primary; ++i)
      for (unsigned r = 0; r < _num_reactions; ++r)
        dmf[i][wrt] += stoichiometry(r, i) * dsec[r][wrt];
 
    // compute dmf[_num_components - 1]
    for (unsigned i = 0; i < _num_primary; ++i)
      dmf[_num_components - 1][wrt] -= dmf[i][wrt];
  }
 
  // Compute the derivatives of the mass_frac wrt the temperature
  // This is used in _dmass_frac_dvar
  std::vector<Real> dmf_dT(_num_components, 0.0);
  for (unsigned i = 0; i < _num_components - 1; ++i)
  {
    for (unsigned r = 0; r < _num_reactions; ++r)
      dmf_dT[i] += stoichiometry(r, i) * dsec_dT[r];
    dmf_dT[_num_components - 1] -= dmf_dT[i];
  }
 
  // compute _dmass_frac_dvar[_qp][_aq_ph] and _dsec_conc_dvar[_qp]
  for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
  {
    // derivative with respect to the "wrt"^th primary species concentration
    if (!_dictator.isPorousFlowVariable(_mf_vars_num[_aq_i + wrt]))
      continue;
    const unsigned pf_wrt = _dictator.porousFlowVariableNum(_mf_vars_num[_aq_i + wrt]);
 
    // run through the mass fractions, building the derivative using dmf
    for (unsigned i = 0; i < _num_components; ++i)
      (*_dmass_frac_dvar)[_qp][_aq_ph][i][pf_wrt] = dmf[i][wrt];
 
    // run through the secondary concentrations, using dsec in the appropriate places
    for (unsigned r = 0; r < _num_reactions; ++r)
      _dsec_conc_dvar[_qp][r][pf_wrt] = dsec[r][wrt];
  }
 
  // use the derivative wrt temperature
  for (unsigned i = 0; i < _num_components; ++i)
    for (unsigned v = 0; v < _num_var; ++v)
      (*_dmass_frac_dvar)[_qp][_aq_ph][i][v] += dmf_dT[i] * _dtemperature_dvar[_qp][v];
  for (unsigned r = 0; r < _num_reactions; ++r)
    for (unsigned v = 0; v < _num_var; ++v)
      _dsec_conc_dvar[_qp][r][v] += dsec_dT[r] * _dtemperature_dvar[_qp][v];
 
  // compute the gradient, if needed
  // NOTE: The derivative d(grad_mass_frac)/d(var) != d(mass_frac)/d(var) * grad_phi
  //       because mass fraction is a nonlinear function of the primary variables
  //       This means that the Jacobian in PorousFlowDispersiveFlux will be wrong
  if (!_nodal_material)
  {
    (*_grad_mass_frac)[_qp][_aq_ph][_num_components - 1] = 0.0;
    for (unsigned comp = 0; comp < _num_components - 1; ++comp)
    {
      (*_grad_mass_frac)[_qp][_aq_ph][comp] = 0.0;
      for (unsigned wrt = 0; wrt < _num_primary; ++wrt)
        (*_grad_mass_frac)[_qp][_aq_ph][comp] +=
            dmf[comp][wrt] * (*_grad_mf_vars[_aq_i + wrt])[_qp];
      (*_grad_mass_frac)[_qp][_aq_ph][_num_components - 1] -= (*_grad_mass_frac)[_qp][_aq_ph][comp];
    }
  }
}

Referenced by initQpStatefulProperties().

◆ computeQpSecondaryConcentrations()

void PorousFlowMassFractionAqueousEquilibriumChemistry::computeQpSecondaryConcentrations ( )

protectedvirtual

Compute the secondary-species concentration as defined by the chemistry Must be overridden by derived classes.

Definition at line 311 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  for (unsigned r = 0; r < _num_reactions; ++r)
  {
    _sec_conc[_qp][r] = 1.0;
    for (unsigned i = 0; i < _num_primary; ++i)
    {
      const Real gamp = _primary_activity_coefficients[i] * (*_mf_vars[_aq_i + i])[_qp];
      if (gamp <= 0.0)
      {
        if (stoichiometry(r, i) < 0.0)
          _sec_conc[_qp][r] = std::numeric_limits<Real>::max();
        else if (stoichiometry(r, i) == 0.0)
          _sec_conc[_qp][r] *= 1.0;
        else
        {
          _sec_conc[_qp][r] = 0.0;
          break;
        }
      }
      else
        _sec_conc[_qp][r] *= std::pow(gamp, stoichiometry(r, i));
    }
    _sec_conc[_qp][r] *=
        (_equilibrium_constants_as_log10 ? std::pow(10.0, (*_equilibrium_constants[r])[_qp])
                                         : (*_equilibrium_constants[r])[_qp]);
    _sec_conc[_qp][r] /= _secondary_activity_coefficients[r];
  }
}

Referenced by computeQpProperties().

◆ dQpSecondaryConcentration_dprimary()

void PorousFlowMassFractionAqueousEquilibriumChemistry::dQpSecondaryConcentration_dprimary	(	unsigned	reaction_num,
		std::vector< Real > &	dsc
	)		const

protectedvirtual

Computes derivative of the secondary concentration with respect to the primary concentrations Must be overridden by derived classes.

Parameters

reaction_num	The reaction number corresponding to the secondary-species concentration
dsc	dsc[i] = d(secondaryConcentration[reaction_num])/d(primary_species[i])

Definition at line 342 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  dsc.assign(_num_primary, 0.0);
 
  /*
   * the derivatives are straightforward if all primary > 0.
   *
   * If more than one primary = 0 then I set the derivatives to zero, even though it could be argued
   * that with certain stoichiometric coefficients you might have derivative = 0/0 and it might be
   * appropriate to set this to a non-zero finite value.
   *
   * If exactly one primary = 0 and its stoichiometry = 1 then the derivative wrt this one is
   * nonzero.
   * If exactly one primary = 0 and its stoichiometry > 1 then all derivatives are zero.
   * If exactly one primary = 0 and its stoichiometry < 1 then the derivative wrt this one is
   * infinity
   */
 
  unsigned zero_count = 0;
  unsigned zero_conc_index = 0;
  findZeroConcentration(zero_conc_index, zero_count);
 
  if (zero_count == 0)
  {
    for (unsigned i = 0; i < _num_primary; ++i)
      dsc[i] = stoichiometry(reaction_num, i) * _sec_conc[_qp][reaction_num] /
               (*_mf_vars[_aq_i + i])[_qp];
  }
  else
  {
    // count the number of primary <= 0, and record the one that's zero
    if (zero_count == 1 and stoichiometry(reaction_num, zero_conc_index) == 1.0)
    {
      Real conc_without_zero = 1.0;
      for (unsigned i = 0; i < _num_primary; ++i)
      {
        if (i == zero_conc_index)
          conc_without_zero *= _primary_activity_coefficients[i];
        else
          conc_without_zero *=
              std::pow(_primary_activity_coefficients[i] * (*_mf_vars[_aq_i + i])[_qp],
                       stoichiometry(reaction_num, i));
      }
      conc_without_zero *= (_equilibrium_constants_as_log10
                                ? std::pow(10.0, (*_equilibrium_constants[reaction_num])[_qp])
                                : (*_equilibrium_constants[reaction_num])[_qp]);
      conc_without_zero /= _secondary_activity_coefficients[reaction_num];
      dsc[zero_conc_index] = conc_without_zero;
    }
    else if (zero_count == 1 && stoichiometry(reaction_num, zero_conc_index) < 1.0)
      dsc[zero_conc_index] = std::numeric_limits<Real>::max();
 
    // all other cases have dsc = 0
  }
}

Referenced by computeQpProperties().

◆ dQpSecondaryConcentration_dT()

Real PorousFlowMassFractionAqueousEquilibriumChemistry::dQpSecondaryConcentration_dT ( unsigned reaction_num ) const

protectedvirtual

Computes derivative of the secondary concentration with respect to the temperature Must be overridden by derived classes.

Parameters

reaction_num	The reaction number corresponding to the secondary-species concentration
dsc	dsc[i] = d(secondaryConcentration[reaction_num])/d(primary_species[i])

Definition at line 400 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  return 0.0;
}

Referenced by computeQpProperties().

◆ findZeroConcentration()

void PorousFlowMassFractionAqueousEquilibriumChemistry::findZeroConcentration	(	unsigned &	zero_conc_index,
		unsigned &	zero_count
	)		const

protected

Checks gamp[i] = _primary_activity_coefficients[i] * (*_primary[i])[qp].

Returns: if all of these are positive, then zero_count = 0, zero_conc_index = 0 if one of these is zero, then zero_count = 1, zero_conc_index = the index of the zero gamp if more than one is zero, then zero_count = 2, and zero_conc_index is the index of the 2nd zero

Definition at line 287 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  zero_count = 0;
  for (unsigned i = 0; i < _num_primary; ++i)
  {
    if (_primary_activity_coefficients[i] * (*_mf_vars[_aq_i + i])[_qp] <= 0.0)
    {
      zero_count += 1;
      zero_conc_index = i;
      if (zero_count > 1)
        return;
    }
  }
  return;
}

Referenced by dQpSecondaryConcentration_dprimary().

◆ initQpSecondaryConcentrations()

void PorousFlowMassFractionAqueousEquilibriumChemistry::initQpSecondaryConcentrations ( )

protectedvirtual

Initialises (at _t_step = 0) the secondary concentrations.

Definition at line 305 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  _sec_conc[_qp].assign(_num_reactions, 0.0);
}

Referenced by computeQpProperties().

◆ initQpStatefulProperties()

void PorousFlowMassFractionAqueousEquilibriumChemistry::initQpStatefulProperties ( )

overrideprotectedvirtual

Reimplemented from PorousFlowMassFractionTempl< is_ad >.

Definition at line 161 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  computeQpProperties();
}

◆ stoichiometry()

Real PorousFlowMassFractionAqueousEquilibriumChemistry::stoichiometry	(	unsigned	reaction_num,
		unsigned	primary_num
	)		const

protected

The stoichiometric coefficient.

Parameters

reaction_num	Reaction number (0, ..., _num_reactions - 1)
primary_num	The number of the primary species (0, ..., _num_primary - 1)

Definition at line 279 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  const unsigned index = reaction_num * _num_primary + primary_num;
  return _reactions[index];
}

Referenced by computeQpProperties(), computeQpSecondaryConcentrations(), and dQpSecondaryConcentration_dprimary().

◆ validParams()

InputParameters PorousFlowMassFractionAqueousEquilibriumChemistry::validParams ( )

static

Definition at line 15 of file PorousFlowMassFractionAqueousEquilibriumChemistry.C.

{
  InputParameters params = PorousFlowMaterialVectorBase::validParams();
  params.addRequiredCoupledVar(
      "mass_fraction_vars",
      "List of variables that represent the mass fractions.  For the aqueous phase these are "
      "concentrations of the primary species with units m^{3}(chemical)/m^{3}(fluid phase).  For "
      "the other phases (if any) these will typically be initialised to zero and will not change "
      "throughout the simulation.  Format is 'f_ph0^c0 f_ph0^c1 f_ph0^c2 ... f_ph0^c(N-2) f_ph1^c0 "
      "f_ph1^c1 fph1^c2 ... fph1^c(N-2) ... fphP^c0 f_phP^c1 fphP^c2 ... fphP^c(N-2)' where "
      "N=number of primary species and P=num_phases, and it is assumed that "
      "f_ph^c(N-1)=1-sum(f_ph^c,{c,0,N-2}) so that f_ph^c(N-1) need not be given.");
  params.addRequiredParam<unsigned>("num_reactions",
                                    "Number of equations in the system of chemical reactions");
  params.addParam<bool>("equilibrium_constants_as_log10",
                        false,
                        "If true, the equilibrium constants are written in their log10 form, eg, "
                        "-2.  If false, the equilibrium constants are written in absolute terms, "
                        "eg, 0.01");
  params.addRequiredCoupledVar("equilibrium_constants",
                               "Equilibrium constant for each equation (dimensionless).  If these "
                               "are temperature dependent AuxVariables, the Jacobian will not be "
                               "exact");
  params.addRequiredParam<std::vector<Real>>(
      "primary_activity_coefficients",
      "Activity coefficients for the primary species (dimensionless) (one for each)");
  params.addRequiredParam<std::vector<Real>>(
      "reactions",
      "A matrix defining the aqueous reactions.  The matrix is entered as a long vector: the first "
      "row is "
      "entered first, followed by the second row, etc.  There should be num_reactions rows.  All "
      "primary species should appear only on the LHS of each reaction (and there should be just "
      "one secondary species on the RHS, by definition) so they may have negative coefficients.  "
      "Each row should have number of primary_concentrations entries, which are the stoichiometric "
      "coefficients.  The first coefficient must always correspond to the first primary species, "
      "etc");
  params.addRequiredParam<std::vector<Real>>(
      "secondary_activity_coefficients",
      "Activity coefficients for the secondary species (dimensionless) (one for each reaction)");
  params.addPrivateParam<std::string>("pf_material_type", "mass_fraction");
  params.addClassDescription(
      "This Material forms a std::vector<std::vector ...> of mass-fractions "
      "(total concentrations of primary species (m^{3}(primary species)/m^{3}(solution)) and since "
      "this is for an aqueous system only, mass-fraction equals volume-fraction) corresponding to "
      "an "
      "aqueous equilibrium chemistry system.  The first mass fraction is the "
      "concentration of the first primary species, etc, and the last mass "
      "fraction is the concentration of H2O.");
  return params;
}

Member Data Documentation

◆ _aq_i

const unsigned int PorousFlowMassFractionAqueousEquilibriumChemistry::_aq_i

protected

Index (into _mf_vars) of the first of the primary species.

Definition at line 95 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Referenced by computeQpProperties(), computeQpSecondaryConcentrations(), dQpSecondaryConcentration_dprimary(), and findZeroConcentration().

◆ _aq_ph

const unsigned int PorousFlowMassFractionAqueousEquilibriumChemistry::_aq_ph

protected

Aqueous phase number.

Definition at line 92 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Referenced by computeQpProperties().

◆ _dmass_frac_dvar

template<bool is_ad>

MaterialProperty<std::vector<std::vector<std::vector<Real> > > >* const PorousFlowMassFractionTempl< is_ad >::_dmass_frac_dvar

protectedinherited

Derivative of the mass fraction matrix with respect to the porous flow variables.

Definition at line 34 of file PorousFlowMassFraction.h.

Referenced by computeQpProperties().

◆ _dsec_conc_dvar

MaterialProperty<std::vector<std::vector<Real> > >& PorousFlowMassFractionAqueousEquilibriumChemistry::_dsec_conc_dvar

protected

Derivative of the secondary concentrations with respect to the porous flow variables.

Definition at line 80 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Referenced by computeQpProperties().

◆ _dtemperature_dvar

const MaterialProperty<std::vector<Real> >& PorousFlowMassFractionAqueousEquilibriumChemistry::_dtemperature_dvar

protected

d(temperature)/(d porflow variable)

Definition at line 86 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Referenced by computeQpProperties().

◆ _equilibrium_constants

std::vector<const VariableValue *> PorousFlowMassFractionAqueousEquilibriumChemistry::_equilibrium_constants

protected

Equilibrium constants (dimensionless)

Definition at line 107 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Referenced by computeQpSecondaryConcentrations(), dQpSecondaryConcentration_dprimary(), and PorousFlowMassFractionAqueousEquilibriumChemistry().

◆ _equilibrium_constants_as_log10

const bool PorousFlowMassFractionAqueousEquilibriumChemistry::_equilibrium_constants_as_log10

protected

Whether the equilibium constants are written in their log10 form, or in absolute terms.

Definition at line 101 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Referenced by computeQpSecondaryConcentrations(), and dQpSecondaryConcentration_dprimary().

◆ _grad_mass_frac

template<bool is_ad>

GenericMaterialProperty<std::vector<std::vector<RealGradient> >, is_ad>* const PorousFlowMassFractionTempl< is_ad >::_grad_mass_frac

protectedinherited

Gradient of the mass fraction matrix at the quad points.

Definition at line 31 of file PorousFlowMassFraction.h.

Referenced by computeQpProperties().

◆ _grad_mf_vars

template<bool is_ad>

std::vector<const GenericVariableGradient<is_ad> *> PorousFlowMassFractionTempl< is_ad >::_grad_mf_vars

protectedinherited

The gradient of the mass-fraction variables.

Definition at line 60 of file PorousFlowMassFraction.h.

Referenced by computeQpProperties(), and PorousFlowMassFractionTempl< is_ad >::PorousFlowMassFractionTempl().

◆ _mass_frac

template<bool is_ad>

GenericMaterialProperty<std::vector<std::vector<Real> >, is_ad>& PorousFlowMassFractionTempl< is_ad >::_mass_frac

protectedinherited

Mass fraction matrix at quadpoint or nodes.

Definition at line 28 of file PorousFlowMassFraction.h.

Referenced by computeQpProperties().

◆ _mf_vars

template<bool is_ad>

std::vector<const GenericVariableValue<is_ad> *> PorousFlowMassFractionTempl< is_ad >::_mf_vars

protectedinherited

The mass-fraction variables.

Definition at line 57 of file PorousFlowMassFraction.h.

Referenced by computeQpProperties(), computeQpSecondaryConcentrations(), dQpSecondaryConcentration_dprimary(), findZeroConcentration(), and PorousFlowMassFractionTempl< is_ad >::PorousFlowMassFractionTempl().

◆ _mf_vars_num

template<bool is_ad>

std::vector<unsigned int> PorousFlowMassFractionTempl< is_ad >::_mf_vars_num

protectedinherited

The variable number of the mass-fraction variables.

Definition at line 54 of file PorousFlowMassFraction.h.

Referenced by computeQpProperties(), and PorousFlowMassFractionTempl< is_ad >::PorousFlowMassFractionTempl().

◆ _num_components

const unsigned int PorousFlowMaterialVectorBase::_num_components

protectedinherited

Number of fluid components.

Definition at line 30 of file PorousFlowMaterialVectorBase.h.

Referenced by computeQpProperties(), PorousFlowAqueousPreDisChemistry::PorousFlowAqueousPreDisChemistry(), PorousFlowDiffusivityBaseTempl< is_ad >::PorousFlowDiffusivityBaseTempl(), and PorousFlowMassFractionTempl< is_ad >::PorousFlowMassFractionTempl().

◆ _num_equilibrium_constants

const unsigned PorousFlowMassFractionAqueousEquilibriumChemistry::_num_equilibrium_constants

protected

Number of equilibrium_constants provided.

Definition at line 104 of file PorousFlowMassFractionAqueousEquilibriumChemistry.h.

Referenced by PorousFlowMassFractionAqueousEquilibriumChemistry().

◆ _num_passed_mf_vars

template<bool is_ad>

const unsigned int PorousFlowMassFractionTempl< is_ad >::_num_passed_mf_vars

protectedinherited

Number of mass-fraction variables provided by the user This needs to be num_phases*(_num_components - 1), since the mass fraction of the final component in each phase is determined as 1 - sum{components}(mass fraction of all other components in the phase)

Definition at line 51 of file PorousFlowMassFraction.h.

Referenced by PorousFlowMassFractionTempl< is_ad >::PorousFlowMassFractionTempl().