GeochemistryTimeIndependentReactor

This UserObject is usually added via a TimeIndependentReactionSolver Action: please see that page for many examples. This UserObject's purpose is to solve a time-independent geochemical system.

Advanced users may wish to add GeochemistryTimeIndependentReactor objects manually to their input files. Here is an example of that:

# This is an example of an input file that does not utilize an action.  Its functionality is the same as HCl.i
# This solves for molalities in a system just containing HCl
[GlobalParams]
  point = '0 0 0'
[]

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx= 1
[]

[Variables]
  [u]
  []
[]

[Kernels]
  [u]
    type = Diffusion
    variable = u
  []
[]

[AuxVariables]
  [solution_temperature]
  []
  [kg_solvent_H2O]
  []
  [activity_H2O]
  []
  [bulk_moles_H2O]
  []
  [pH]
  []
  [molal_H+]
  []
  [molal_Cl-]
  []
  [molal_HCl]
  []
  [molal_OH-]
  []
  [mg_per_kg_H+]
  []
  [mg_per_kg_Cl-]
  []
  [mg_per_kg_HCl]
  []
  [mg_per_kg_OH-]
  []
  [activity_H+]
  []
  [activity_Cl-]
  []
  [activity_HCl]
  []
  [activity_OH-]
  []
  [bulk_moles_H+]
  []
  [bulk_moles_Cl-]
  []
  [bulk_moles_HCl]
  []
  [bulk_moles_OH-]
  []
[]

[AuxKernels]
  [solution_temperature]
    type = GeochemistryQuantityAux
    species = 'H+'
    reactor = reactor
    variable = solution_temperature
    quantity = temperature
  []
  [kg_solvent_H2O]
    type = GeochemistryQuantityAux
    species = 'H2O'
    reactor = reactor
    variable = kg_solvent_H2O
    quantity = molal
  []
  [activity_H2O]
    type = GeochemistryQuantityAux
    species = 'H2O'
    reactor = reactor
    variable = activity_H2O
    quantity = activity
  []
  [bulk_moles_H2O]
    type = GeochemistryQuantityAux
    species = 'H2O'
    reactor = reactor
    variable = bulk_moles_H2O
    quantity = bulk_moles
  []
  [pH]
    type = GeochemistryQuantityAux
    species = 'H+'
    reactor = reactor
    variable = pH
    quantity = neglog10a
  []
  [molal_H+]
    type = GeochemistryQuantityAux
    species = 'H+'
    reactor = reactor
    variable = 'molal_H+'
    quantity = molal
  []
  [molal_Cl-]
    type = GeochemistryQuantityAux
    species = 'Cl-'
    reactor = reactor
    variable = 'molal_Cl-'
    quantity = molal
  []
  [molal_HCl]
    type = GeochemistryQuantityAux
    species = 'HCl'
    reactor = reactor
    variable = 'molal_HCl'
    quantity = molal
  []
  [molal_OH-]
    type = GeochemistryQuantityAux
    species = 'OH-'
    reactor = reactor
    variable = 'molal_OH-'
    quantity = molal
  []
  [mg_per_kg_H+]
    type = GeochemistryQuantityAux
    species = 'H+'
    reactor = reactor
    variable = 'mg_per_kg_H+'
    quantity = mg_per_kg
  []
  [mg_per_kg_Cl-]
    type = GeochemistryQuantityAux
    species = 'Cl-'
    reactor = reactor
    variable = 'mg_per_kg_Cl-'
    quantity = mg_per_kg
  []
  [mg_per_kg_HCl]
    type = GeochemistryQuantityAux
    species = 'HCl'
    reactor = reactor
    variable = 'mg_per_kg_HCl'
    quantity = mg_per_kg
  []
  [mg_per_kg_OH-]
    type = GeochemistryQuantityAux
    species = 'OH-'
    reactor = reactor
    variable = 'mg_per_kg_OH-'
    quantity = mg_per_kg
  []
  [activity_H+]
    type = GeochemistryQuantityAux
    species = 'H+'
    reactor = reactor
    variable = 'activity_H+'
    quantity = activity
  []
  [activity_Cl-]
    type = GeochemistryQuantityAux
    species = 'Cl-'
    reactor = reactor
    variable = 'activity_Cl-'
    quantity = activity
  []
  [activity_HCl]
    type = GeochemistryQuantityAux
    species = 'HCl'
    reactor = reactor
    variable = 'activity_HCl'
    quantity = activity
  []
  [activity_OH-]
    type = GeochemistryQuantityAux
    species = 'OH-'
    reactor = reactor
    variable = 'activity_OH-'
    quantity = activity
  []
  [bulk_moles_H+]
    type = GeochemistryQuantityAux
    species = 'H+'
    reactor = reactor
    variable = 'bulk_moles_H+'
    quantity = bulk_moles
  []
  [bulk_moles_Cl-]
    type = GeochemistryQuantityAux
    species = 'Cl-'
    reactor = reactor
    variable = 'bulk_moles_Cl-'
    quantity = bulk_moles
  []
  [bulk_moles_HCl]
    type = GeochemistryQuantityAux
    species = 'HCl'
    reactor = reactor
    variable = 'bulk_moles_HCl'
    quantity = bulk_moles
  []
  [bulk_moles_OH-]
    type = GeochemistryQuantityAux
    species = 'OH-'
    reactor = reactor
    variable = 'bulk_moles_OH-'
    quantity = bulk_moles
  []
[]

[Postprocessors]
  [pH]
    type = PointValue
    variable = 'pH'
  []
  [solvent_mass]
    type = PointValue
    variable = 'kg_solvent_H2O'
  []
  [molal_Cl-]
    type = PointValue
    variable = 'molal_Cl-'
  []
  [mg_per_kg_HCl]
    type = PointValue
    variable = 'mg_per_kg_HCl'
  []
  [activity_OH-]
    type = PointValue
    variable = 'activity_OH-'
  []
  [bulk_H+]
    type = PointValue
    variable = 'bulk_moles_H+'
  []
  [temperature]
    type = PointValue
    variable = 'solution_temperature'
  []
[]

[Executioner]
  type = Steady
[]

[UserObjects]
  [definition]
    type = GeochemicalModelDefinition
    database_file = "../../../database/moose_geochemdb.json"
    basis_species = "H2O H+ Cl-"
    piecewise_linear_interpolation = true # to reproduce the GWB result
  []
  [reactor]
    type = GeochemistryTimeDependentReactor
    model_definition = definition
    charge_balance_species = "Cl-"
    constraint_species = "H2O              H+            Cl-"
    constraint_value = "  1.0              -2            1E-2"
    constraint_meaning = "kg_solvent_water log10activity bulk_composition"
    constraint_unit = "   kg               dimensionless moles"
    ramp_max_ionic_strength_initial = 0 # max_ionic_strength in such a simple problem does not need ramping
    abs_tol = 1E-15
  []
  [nnn]
    type = NearestNodeNumberUO
  []
[]

[Outputs]
  csv = true
  [console_output]
    type = GeochemistryConsoleOutput
    geochemistry_reactor = reactor
    nearest_node_number_UO = nnn
    solver_info = true
    execute_on = initial
  []
[]
(modules/geochemistry/test/tests/equilibrium_models/HCl_no_action.i)

Input Parameters

  • charge_balance_speciesCharge balance will be enforced on this basis species. This means that its bulk mole number may be changed from the initial value you provide in order to ensure charge neutrality. After the initial swaps have been performed, this must be in the basis, and it must be provided with a bulk_composition constraint_meaning.

    C++ Type:std::string

    Controllable:No

    Description:Charge balance will be enforced on this basis species. This means that its bulk mole number may be changed from the initial value you provide in order to ensure charge neutrality. After the initial swaps have been performed, this must be in the basis, and it must be provided with a bulk_composition constraint_meaning.

  • constraint_meaningMeanings of the numerical values given in constraint_value. kg_solvent_water: can only be applied to H2O and units must be kg. bulk_composition: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species but not in kinetic species, and units must be moles or mass (kg, g, etc). bulk_composition_with_kinetic: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species and in kinetic species, and units must be moles or mass (kg, g, etc). free_concentration: can be applied to all basis species that are not gas and not H2O and not mineral, and represents the total amount of the basis species existing freely (not as secondary species) within the solution, and units must be molal or mass_per_kg_solvent. free_mineral: can be applied to all mineral basis species, and represents the total amount of the mineral existing freely (precipitated) within the solution, and units must be moles, mass or cm3. activity and log10activity: can be applied to basis species that are not gas and not mineral and not sorbing sites, and represents the activity of the basis species (recall pH = -log10activity), and units must be dimensionless. fugacity and log10fugacity: can be applied to gases, and units must be dimensionless

    C++ Type:MultiMooseEnum

    Options:kg_solvent_water, bulk_composition, bulk_composition_with_kinetic, free_concentration, free_mineral, activity, log10activity, fugacity, log10fugacity

    Controllable:No

    Description:Meanings of the numerical values given in constraint_value. kg_solvent_water: can only be applied to H2O and units must be kg. bulk_composition: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species but not in kinetic species, and units must be moles or mass (kg, g, etc). bulk_composition_with_kinetic: can be applied to all non-gas species, and represents the total amount of the basis species contained as free species as well as the amount found in secondary species and in kinetic species, and units must be moles or mass (kg, g, etc). free_concentration: can be applied to all basis species that are not gas and not H2O and not mineral, and represents the total amount of the basis species existing freely (not as secondary species) within the solution, and units must be molal or mass_per_kg_solvent. free_mineral: can be applied to all mineral basis species, and represents the total amount of the mineral existing freely (precipitated) within the solution, and units must be moles, mass or cm3. activity and log10activity: can be applied to basis species that are not gas and not mineral and not sorbing sites, and represents the activity of the basis species (recall pH = -log10activity), and units must be dimensionless. fugacity and log10fugacity: can be applied to gases, and units must be dimensionless

  • constraint_speciesNames of the species that have their values fixed to constraint_value with meaning constraint_meaning. All basis species (after swap_into_basis and swap_out_of_basis) must be provided with exactly one constraint. These constraints are used to compute the configuration during the initial problem setup, and in time-dependent simulations they may be modified as time progresses.

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Names of the species that have their values fixed to constraint_value with meaning constraint_meaning. All basis species (after swap_into_basis and swap_out_of_basis) must be provided with exactly one constraint. These constraints are used to compute the configuration during the initial problem setup, and in time-dependent simulations they may be modified as time progresses.

  • constraint_unitUnits of the numerical values given in constraint_value. Dimensionless: should only be used for activity or fugacity constraints. Moles: mole number. Molal: moles per kg solvent water. kg: kilograms. g: grams. mg: milligrams. ug: micrograms. kg_per_kg_solvent: kilograms per kg solvent water. g_per_kg_solvent: grams per kg solvent water. mg_per_kg_solvent: milligrams per kg solvent water. ug_per_kg_solvent: micrograms per kg solvent water. cm3: cubic centimeters

    C++ Type:MultiMooseEnum

    Options:dimensionless, moles, molal, kg, g, mg, ug, kg_per_kg_solvent, g_per_kg_solvent, mg_per_kg_solvent, ug_per_kg_solvent, cm3

    Controllable:No

    Description:Units of the numerical values given in constraint_value. Dimensionless: should only be used for activity or fugacity constraints. Moles: mole number. Molal: moles per kg solvent water. kg: kilograms. g: grams. mg: milligrams. ug: micrograms. kg_per_kg_solvent: kilograms per kg solvent water. g_per_kg_solvent: grams per kg solvent water. mg_per_kg_solvent: milligrams per kg solvent water. ug_per_kg_solvent: micrograms per kg solvent water. cm3: cubic centimeters

  • constraint_valueNumerical value of the containts on constraint_species

    C++ Type:std::vector<double>

    Controllable:No

    Description:Numerical value of the containts on constraint_species

  • model_definitionThe name of the GeochemicalModelDefinition user object.

    C++ Type:UserObjectName

    Controllable:No

    Description:The name of the GeochemicalModelDefinition user object.

Required Parameters

  • abs_tol1e-10If the residual of the algebraic system (measured in mol) is lower than this value, the Newton process (that finds the aqueous configuration) is deemed to have converged

    Default:1e-10

    C++ Type:double

    Controllable:No

    Description:If the residual of the algebraic system (measured in mol) is lower than this value, the Newton process (that finds the aqueous configuration) is deemed to have converged

  • blockThe list of blocks (ids or names) that this object will be applied

    C++ Type:std::vector<SubdomainName>

    Controllable:No

    Description:The list of blocks (ids or names) that this object will be applied

  • boundaryThe list of boundaries (ids or names) from the mesh where this boundary condition applies

    C++ Type:std::vector<BoundaryName>

    Controllable:No

    Description:The list of boundaries (ids or names) from the mesh where this boundary condition applies

  • execute_onFINALThe list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM, ALWAYS.

    Default:FINAL

    C++ Type:ExecFlagEnum

    Options:NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM, ALWAYS

    Controllable:No

    Description:The list of flag(s) indicating when this object should be executed, the available options include NONE, INITIAL, LINEAR, NONLINEAR, TIMESTEP_END, TIMESTEP_BEGIN, FINAL, CUSTOM, ALWAYS.

  • extra_iterations_to_make_consistent0Extra iterations to make the molalities, activities, etc, consistent before commencing the Newton process to find the aqueous configuration

    Default:0

    C++ Type:unsigned int

    Controllable:No

    Description:Extra iterations to make the molalities, activities, etc, consistent before commencing the Newton process to find the aqueous configuration

  • ionic_str_using_basis_onlyFalseIf set to true, ionic strength and stoichiometric ionic strength will be computed using only the basis molalities, ignoring molalities of equilibrium species. Since basis molality is usually greater than equilibrium molality, and the whole Debye-Huckel concept of activity coefficients depending on ionic strength is only approximate in practice, setting this parameter true often results in a reasonable approximation. It can aid in convergence since it eliminates problems associated with unphysical huge equilibrium molalities that can occur during Newton iteration to the solution

    Default:False

    C++ Type:bool

    Controllable:No

    Description:If set to true, ionic strength and stoichiometric ionic strength will be computed using only the basis molalities, ignoring molalities of equilibrium species. Since basis molality is usually greater than equilibrium molality, and the whole Debye-Huckel concept of activity coefficients depending on ionic strength is only approximate in practice, setting this parameter true often results in a reasonable approximation. It can aid in convergence since it eliminates problems associated with unphysical huge equilibrium molalities that can occur during Newton iteration to the solution

  • max_initial_residual1000Attempt to alter the initial-guess molalities so that the initial residual for the Newton process (that finds the aqueous configuration) is less than this number of moles

    Default:1000

    C++ Type:double

    Controllable:No

    Description:Attempt to alter the initial-guess molalities so that the initial residual for the Newton process (that finds the aqueous configuration) is less than this number of moles

  • max_ionic_strength3Maximum value of ionic strength

    Default:3

    C++ Type:double

    Controllable:No

    Description:Maximum value of ionic strength

  • max_iter100Maximum number of Newton iterations allowed when finding the aqueous configuration

    Default:100

    C++ Type:unsigned int

    Controllable:No

    Description:Maximum number of Newton iterations allowed when finding the aqueous configuration

  • max_swaps_allowed20Maximum number of swaps allowed during a single attempt at finding the aqueous configuration. Usually only a handful of swaps are used: this parameter prevents endless cyclic swapping that prevents the algorithm from progressing

    Default:20

    C++ Type:unsigned int

    Controllable:No

    Description:Maximum number of swaps allowed during a single attempt at finding the aqueous configuration. Usually only a handful of swaps are used: this parameter prevents endless cyclic swapping that prevents the algorithm from progressing

  • min_initial_molality1e-20Minimum value of the initial-guess molality used in the Newton process to find the aqueous configuration

    Default:1e-20

    C++ Type:double

    Controllable:No

    Description:Minimum value of the initial-guess molality used in the Newton process to find the aqueous configuration

  • prevent_precipitationMineral species in this list will be prevented from precipitating, irrespective of their saturation index (unless they are initially in the basis)

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Mineral species in this list will be prevented from precipitating, irrespective of their saturation index (unless they are initially in the basis)

  • ramp_max_ionic_strength_initial20The number of iterations over which to progressively increase the maximum ionic strength (from zero to max_ionic_strength) during the initial equilibration. Increasing this can help in convergence of the Newton process, at the cost of spending more time finding the aqueous configuration.

    Default:20

    C++ Type:unsigned int

    Controllable:No

    Description:The number of iterations over which to progressively increase the maximum ionic strength (from zero to max_ionic_strength) during the initial equilibration. Increasing this can help in convergence of the Newton process, at the cost of spending more time finding the aqueous configuration.

  • rel_tol1e-200If the residual of the algebraic system (measured in mol) is lower than this value times the initial residual, the Newton process (that finds the aqueous configuration) is deemed to have converged

    Default:1e-200

    C++ Type:double

    Controllable:No

    Description:If the residual of the algebraic system (measured in mol) is lower than this value times the initial residual, the Newton process (that finds the aqueous configuration) is deemed to have converged

  • stoichiometric_ionic_str_using_Cl_onlyFalseIf set to true, the stoichiometric ionic strength will be set equal to Cl- molality (or max_ionic_strength if the Cl- molality is too big). This flag overrides ionic_str_using_basis_molality_only

    Default:False

    C++ Type:bool

    Controllable:No

    Description:If set to true, the stoichiometric ionic strength will be set equal to Cl- molality (or max_ionic_strength if the Cl- molality is too big). This flag overrides ionic_str_using_basis_molality_only

  • stoichiometry_tolerance1e-06Swapping involves inverting matrices via a singular value decomposition. During this process: (1) if abs(singular value) < stoi_tol * L1norm(singular values), then the matrix is deemed singular (so the basis swap is deemed invalid); (2) if abs(any stoichiometric coefficient) < stoi_tol then it is set to zero.

    Default:1e-06

    C++ Type:double

    Controllable:No

    Description:Swapping involves inverting matrices via a singular value decomposition. During this process: (1) if abs(singular value) < stoi_tol * L1norm(singular values), then the matrix is deemed singular (so the basis swap is deemed invalid); (2) if abs(any stoichiometric coefficient) < stoi_tol then it is set to zero.

  • swap_into_basisSpecies that should be removed from the model_definition's equilibrium species list and added to the basis. There must be the same number of species in swap_out_of_basis and swap_into_basis. These swaps are performed before any other computations during the initial problem setup. If this list contains more than one species, the swapping is performed one-by-one, starting with the first pair (swap_out_of_basis[0] and swap_into_basis[0]), then the next pair, etc

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Species that should be removed from the model_definition's equilibrium species list and added to the basis. There must be the same number of species in swap_out_of_basis and swap_into_basis. These swaps are performed before any other computations during the initial problem setup. If this list contains more than one species, the swapping is performed one-by-one, starting with the first pair (swap_out_of_basis[0] and swap_into_basis[0]), then the next pair, etc

  • swap_out_of_basisSpecies that should be removed from the model_definition's basis and be replaced with the swap_into_basis species

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Species that should be removed from the model_definition's basis and be replaced with the swap_into_basis species

  • swap_threshold0.1If the molality of a basis species in the algebraic system falls below swap_threshold * abs_tol then it is swapped out of the basis. The dimensions of swap_threshold are 1/kg(solvent water)

    Default:0.1

    C++ Type:double

    Controllable:No

    Description:If the molality of a basis species in the algebraic system falls below swap_threshold * abs_tol then it is swapped out of the basis. The dimensions of swap_threshold are 1/kg(solvent water)

  • temperature25The temperature (degC) of the aqueous solution

    Default:25

    C++ Type:double

    Controllable:No

    Description:The temperature (degC) of the aqueous solution

  • unique_node_executeFalseWhen false (default), block restricted objects will have the execute method called multiple times on a single node if the node lies on a interface between two subdomains.

    Default:False

    C++ Type:bool

    Controllable:No

    Description:When false (default), block restricted objects will have the execute method called multiple times on a single node if the node lies on a interface between two subdomains.

Optional Parameters

  • allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).

    Default:False

    C++ Type:bool

    Controllable:No

    Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).

  • control_tagsAdds user-defined labels for accessing object parameters via control logic.

    C++ Type:std::vector<std::string>

    Controllable:No

    Description:Adds user-defined labels for accessing object parameters via control logic.

  • enableTrueSet the enabled status of the MooseObject.

    Default:True

    C++ Type:bool

    Controllable:Yes

    Description:Set the enabled status of the MooseObject.

  • force_postauxFalseForces the UserObject to be executed in POSTAUX

    Default:False

    C++ Type:bool

    Controllable:No

    Description:Forces the UserObject to be executed in POSTAUX

  • force_preauxFalseForces the UserObject to be executed in PREAUX

    Default:False

    C++ Type:bool

    Controllable:No

    Description:Forces the UserObject to be executed in PREAUX

  • force_preicFalseForces the UserObject to be executed in PREIC during initial setup

    Default:False

    C++ Type:bool

    Controllable:No

    Description:Forces the UserObject to be executed in PREIC during initial setup

  • seed0The seed for the master random number generator

    Default:0

    C++ Type:unsigned int

    Controllable:No

    Description:The seed for the master random number generator

  • use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

    Default:False

    C++ Type:bool

    Controllable:No

    Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters