GeochemicalSystem.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "MultiMooseEnum.h"
#include "PertinentGeochemicalSystem.h"
#include "GeochemistrySpeciesSwapper.h"
#include "GeochemistryActivityCoefficientsDebyeHuckel.h"
#include "GeochemistryConstants.h"
#include "GeochemistryUnitConverter.h"

class GeochemicalSystem
{
public:
  enum class ConstraintMeaningEnum
  {
    MOLES_BULK_WATER,
    KG_SOLVENT_WATER,
    MOLES_BULK_SPECIES,
    FREE_MOLALITY,
    FREE_MOLES_MINERAL_SPECIES,
    FUGACITY,
    ACTIVITY
  };

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

  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);

  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)
  {
    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;
  };

  GeochemicalSystem(const GeochemicalSystem & src) = default;

  bool operator==(const GeochemicalSystem & rhs) const
  {
    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);
  }

  unsigned getNumInBasis() const;

  unsigned getNumInEquilibrium() const;

  unsigned getNumKinetic() const;

  void initialize();

  unsigned getChargeBalanceBasisIndex() const;

  bool alterChargeBalanceSpecies(Real threshold_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);

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

  void updateOldWithCurrent(const DenseVector<Real> & mole_additions);

  unsigned getNumInAlgebraicSystem() const;

  unsigned getNumBasisInAlgebraicSystem() const;

  unsigned getNumSurfacePotentials() const;

  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;

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

  DenseVector<Real> getBulkOldInOriginalBasis() const;

  DenseVector<Real> getTransportedBulkInOriginalBasis() const;

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

  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;

  const ModelGeochemicalDatabase & getModelGeochemicalDatabase() const;

  void setModelGeochemicalDatabase(const ModelGeochemicalDatabase & mgd);

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

  void setAlgebraicVariables(const DenseVector<Real> & algebraic_var);

  void enforceChargeBalance();

  bool revertToOriginalChargeBalanceSpecies();

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

  void performSwap(unsigned swap_out_of_basis, unsigned swap_into_basis);

  Real getIonicStrength() const;

  Real getStoichiometricIonicStrength() const;

  Real getSurfacePotential(unsigned sp) const;

  Real getSurfaceCharge(unsigned sp) const;

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

  Real getTemperature() const;

  void setTemperature(Real temperature);

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

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

  void closeSystem();

  void changeConstraintToBulk(unsigned basis_ind);

  void changeConstraintToBulk(unsigned basis_ind, Real value);

  void addToBulkMoles(unsigned basis_ind, Real value);

  void setConstraintValue(unsigned basis_ind, Real value);

  void computeConsistentConfiguration();

  const std::string & getOriginalRedoxLHS() const;

  void performSwapNoCheck(unsigned swap_out_of_basis, unsigned swap_into_basis);

  const GeochemistrySpeciesSwapper & getSwapper() const;

  void setMineralRelatedFreeMoles(Real value);

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

  Real getEquilibriumActivity(unsigned eqm_ind) const;

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

private:
  ModelGeochemicalDatabase & _mgd;
  const unsigned _num_basis;
  const unsigned _num_eqm;
  const unsigned _num_redox;
  const unsigned _num_surface_pot;
  const unsigned _num_kin;
  GeochemistrySpeciesSwapper & _swapper;
  const std::vector<std::string> _swap_out;
  const std::vector<std::string> _swap_in;
  GeochemistryActivityCoefficients & _gac;
  GeochemistryIonicStrength & _is;
  std::string _charge_balance_species;
  const std::string _original_charge_balance_species;
  unsigned _charge_balance_basis_index;
  std::vector<std::string> _constrained_species;
  std::vector<Real> _constraint_value;
  std::vector<Real> _original_constraint_value;
  std::vector<GeochemistryUnitConverter::GeochemistryUnit> _constraint_unit;
  std::vector<ConstraintUserMeaningEnum> _constraint_user_meaning;
  std::vector<ConstraintMeaningEnum> _constraint_meaning;
  std::vector<Real> _eqm_log10K;
  std::vector<Real> _redox_log10K;
  std::vector<Real> _kin_log10K;
  unsigned _num_basis_in_algebraic_system;
  unsigned _num_in_algebraic_system;
  std::vector<bool> _in_algebraic_system;
  std::vector<unsigned> _algebraic_index;
  std::vector<unsigned> _basis_index;
  std::vector<Real> _bulk_moles_old;
  std::vector<Real> _basis_molality;
  std::vector<bool> _basis_activity_known;
  std::vector<Real> _basis_activity;
  std::vector<Real> _eqm_molality;
  std::vector<Real> _basis_activity_coef;
  std::vector<Real> _eqm_activity_coef;
  std::vector<Real> _eqm_activity;
  std::vector<Real> _surface_pot_expr;
  std::vector<Real> _sorbing_surface_area;
  std::vector<Real> _kin_moles;
  std::vector<Real> _kin_moles_old;
  unsigned _iters_to_make_consistent;
  Real _temperature;
  Real _min_initial_molality;
  std::string _original_redox_lhs;

  void buildTemperatureDependentQuantities(Real temperature);

  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;

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

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

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

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

  Real log10ActivityProduct(unsigned eqm_j) const;

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

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

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

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

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

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

  void setChargeBalanceSpecies(unsigned new_charge_balance_index);

  Real surfaceSorptionModifier(unsigned eqm_j) const;

  Real surfacePotPrefactor(unsigned sp) const;

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

  Real computeBulkFromMolalities(unsigned basis_ind) const;

  void alterSystemBecauseBulkChanged();
};