Skip to content

wrappers

outer_fixed_point 

outer_fixed_point(
    *,
    ionic_strength_dependence: bool = False,
    charges: NDArray | None = None,
    reference_ionic_str_species: NDArray | None = None,
    reference_ionic_str_solids: NDArray | None = None,
    dbh_values: dict[str, dict[str, NDArray]] | None = None,
    independent_component_activity: NDArray | None = None
)

Decorator for implementing the outer fixed point iteration method as developed described by Carrayrou, et al. 1.

Given a function capable of solving the equilibrium equations for a given set of species concentrations, the outer fixed point method fixes the activity coefficients and solve the equilibrium equations. After convergence is achieved for the species concentrations, the activity coefficients are updated and the process is repeated. If the activity coefficient do not change significantly between iterations, the process is considered converged and both species concentrations and activity coefficients are returned.

  1. Carrayrou, et al.: AIChE Journal 2022, 68, e17391 ↩

Parameters:

Name Type Description Default
ionic_strength_dependence bool

Flag indicating whether to consider ionic strength dependence, by default False.

 False
charges ndarray or None

Array of charges for each species, by default None.

 None
reference_ionic_str_species ndarray or None

Array of reference ionic strength for each species, by default None.

 None
reference_ionic_str_solids ndarray or None

Array of reference ionic strength for each solid, by default None.

 None
dbh_values dict[str, dict[str, ndarray]] or None

Dictionary containing the debye huckel parameters used for the species and solids, by default None.

 None
independent_component_activity ndarray or None

Array of independent component activities, by default None.

 None

Returns:

Name Type Description
decorator function

The outer fixed point decorator.

Examples:

Apply the outer_fixed_point decorator to a function:

@outer_fixed_point(
    ionic_strength_dependence=True,
    charges=[1,0,-1],
    ref_ionic_strength_species=[0.1,0.1,0.1],
    ref_ionic_strength_solids=[0.1,0.1,0.1],
    dbh_values={"species": {...}, "solids": {...}
    independent_component_activity=None
    )
def equilibrium_equations(concentration, log_beta, log_ks):
    # implementation of equilibrium equations
    pass

The resulting function can then be called with to compute the species concentrations at equilibrium:

result, final_log_beta, final_log_ks = equilibrium_equations(concentration, log_beta, log_ks)

If ionic_strength_dependence is set to True, the decorator will perform the outer fixed point iteration method, otherwise it will return the result of the decorated function.

Source code in src/libeq/outer_fixed_point/wrappers.py 
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def outer_fixed_point(
    *,
    ionic_strength_dependence: bool = False,
    charges: NDArray | None = None,
    reference_ionic_str_species: NDArray | None = None,
    reference_ionic_str_solids: NDArray | None = None,
    dbh_values: dict[str, dict[str, NDArray]] | None = None,
    independent_component_activity: NDArray | None = None,
):
    """
    Decorator for implementing the outer fixed point iteration method as developed described by Carrayrou, *et al.* [^1].

    Given a function capable of solving the equilibrium equations for a given set of species concentrations,
    the outer fixed point method fixes the activity coefficients and solve the equilibrium equations.
    After convergence is achieved for the species concentrations, the activity coefficients are updated and the process is repeated.
    If the activity coefficient do not change significantly between iterations, the process is considered converged and both species
    concentrations and activity coefficients are returned.

    [^1]: Carrayrou, *et al.*: *AIChE Journal* 2022, **68**, e17391

    Parameters
    ----------
    ionic_strength_dependence : bool, optional
        Flag indicating whether to consider ionic strength dependence, by default False.
    charges : numpy.ndarray or None, optional
        Array of charges for each species, by default None.
    reference_ionic_str_species : numpy.ndarray or None, optional
        Array of reference ionic strength for each species, by default None.
    reference_ionic_str_solids : numpy.ndarray or None, optional
        Array of reference ionic strength for each solid, by default None.
    dbh_values : dict[str, dict[str, numpy.ndarray]] or None, optional
        Dictionary containing the debye huckel parameters used for the species and solids, by default None.
    independent_component_activity : numpy.ndarray or None, optional
        Array of independent component activities, by default None.

    Returns
    -------
    decorator : function
        The outer fixed point decorator.

    Examples
    --------
    Apply the outer_fixed_point decorator to a function:

    ```python
    @outer_fixed_point(
        ionic_strength_dependence=True,
        charges=[1,0,-1],
        ref_ionic_strength_species=[0.1,0.1,0.1],
        ref_ionic_strength_solids=[0.1,0.1,0.1],
        dbh_values={"species": {...}, "solids": {...}
        independent_component_activity=None
        )
    def equilibrium_equations(concentration, log_beta, log_ks):
        # implementation of equilibrium equations
        pass
    ```

    The resulting function can then be called with to compute the species concentrations at equilibrium:

    ```python
    result, final_log_beta, final_log_ks = equilibrium_equations(concentration, log_beta, log_ks)
    ```

    If ionic_strength_dependence is set to True, the decorator will perform the outer fixed point iteration method,
    otherwise it will return the result of the decorated function.

    """

    def decorator(func):
        def unwrapped(concentration, **kwargs):
            log_beta = kwargs.pop("log_beta")
            log_ks = kwargs.pop("log_ks")
            return (
                func(concentration, log_beta=log_beta, log_ks=log_ks, **kwargs),
                log_beta,
                log_ks,
            )

        def wrapper(
            concentration,
            **kwargs,
        ):
            if independent_component_activity is not None:
                transposed_activity = independent_component_activity[:, np.newaxis]
            else:
                transposed_activity = None

            log_beta = kwargs.pop("log_beta")
            log_ks = kwargs.pop("log_ks")
            og_log_beta = log_beta.copy()
            og_log_ks = log_ks.copy()
            stoichiometry = kwargs["stoichiometry"]

            n_components = stoichiometry.shape[0]
            n_species = stoichiometry.shape[1]

            concentrations = species_concentration(
                concentration, log_beta, stoichiometry, full=True
            )
            ionic = _ionic_fn(
                _select_species_concentration(concentrations, n_components, n_species),
                charges,
                independent_component_activity=transposed_activity,
            )
            log_beta = _update_formation_constants(
                og_log_beta,
                ionic,
                reference_ionic_str_species,
                dbh_values["species"],
            )
            log_ks = _update_solubility_products(
                og_log_ks, ionic, reference_ionic_str_solids, dbh_values["solids"]
            )

            result = concentration
            iterations = 0
            while True:
                iterations += 1
                # Call the decorated function
                result = func(
                    result,
                    log_beta=log_beta,
                    log_ks=log_ks,
                    **kwargs,
                )

                # Code to be executed after the decorated function
                concentrations = species_concentration(
                    result, log_beta, stoichiometry, full=True
                )
                ionic = _ionic_fn(
                    _select_species_concentration(
                        concentrations, n_components, n_species
                    ),
                    charges,
                    independent_component_activity=transposed_activity,
                )

                old_log_beta = log_beta
                old_log_ks = log_ks

                log_beta = _update_formation_constants(
                    og_log_beta,
                    ionic,
                    reference_ionic_str_species,
                    dbh_values["species"],
                )
                log_ks = _update_solubility_products(
                    og_log_ks, ionic, reference_ionic_str_solids, dbh_values["solids"]
                )

                if _check_outer_point_convergence(
                    log_beta, old_log_beta, log_ks, old_log_ks
                ):
                    # print(func.__name__)
                    # print(f"Outer converged in {iterations} iterations")
                    # print("------------------------")
                    break
            return result, old_log_beta, old_log_ks

        def _distribution_ionic(
            concentration: NDArray,
            charges: NDArray,
            *,
            independent_component_activity: NDArray,
        ) -> NDArray:
            return _ionic(concentration, charges) + independent_component_activity

        def _ionic(concentration: NDArray, charges: NDArray, **kwargs) -> NDArray:
            return 0.5 * (concentration * (charges**2)).sum(axis=1, keepdims=True)

        if ionic_strength_dependence:
            if independent_component_activity is None:
                _ionic_fn = _ionic
            else:
                _ionic_fn = _distribution_ionic

            return wrapper
        else:
            return unwrapped

    return decorator