Coverage for src/gwtransport/logremoval.py: 87%
31 statements
« prev ^ index » next coverage.py v7.8.2, created at 2025-06-06 15:45 +0000
« prev ^ index » next coverage.py v7.8.2, created at 2025-06-06 15:45 +0000
1"""
2Functions for calculating log removal efficiency in water treatment systems.
4This module provides utilities to calculate log removal values for different
5configurations of water treatment systems, particularly focusing on parallel flow
6arrangements where multiple treatment processes operate simultaneously on different
7fractions of the total flow.
9Log removal is a standard measure in water treatment that represents the
10reduction of pathogen concentration on a logarithmic scale. For example,
11a log removal of 3 represents a 99.9% reduction in pathogen concentration.
13Functions:
14calculate_parallel_log_removal : Calculate the weighted average log removal for
15 a system with parallel flows.
17Notes
18-----
19For systems in series, log removals are typically summed directly, while for
20parallel systems, a weighted average based on flow distribution is required.
21"""
23import numpy as np
24from scipy import stats
25from scipy.special import digamma, gamma
28def parallel_mean(log_removals, flow_fractions=None):
29 """
30 Calculate the weighted average log removal for a system with parallel flows.
32 This function computes the overall log removal efficiency of a parallel
33 filtration system. If flow_fractions is not provided, it assumes equal
34 distribution of flow across all paths.
36 The calculation uses the formula:
38 Total Log Removal = -log₁₀(sum(F_i * 10^(-LR_i)))
40 Where:
41 - F_i = fraction of flow through system i (decimal, sum to 1.0)
42 - LR_i = log removal of system i
44 Parameters
45 ----------
46 log_removals : array_like
47 Array of log removal values for each parallel flow.
48 Each value represents the log₁₀ reduction of pathogens.
50 flow_fractions : array_like, optional
51 Array of flow fractions for each parallel flow.
52 Must sum to 1.0 and be the same length as log_removals.
53 If None, equal flow distribution is assumed (default is None).
55 Returns
56 -------
57 float
58 The combined log removal value for the parallel system.
60 Raises
61 ------
62 ValueError
63 If log_removals is empty
64 If flow_fractions is provided and does not sum to 1.0 (within tolerance)
65 If flow_fractions is provided and has different length than log_removals
67 Notes
68 -----
69 Log removal is a logarithmic measure of pathogen reduction:
70 - Log 1 = 90% reduction
71 - Log 2 = 99% reduction
72 - Log 3 = 99.9% reduction
74 For parallel flows, the combined removal is typically less effective
75 than the best individual removal but better than the worst.
77 Examples
78 --------
79 >>> import numpy as np
80 >>> # Three parallel streams with equal flow and log removals of 3, 4, and 5
81 >>> log_removals = np.array([3, 4, 5])
82 >>> calculate_parallel_log_removal(log_removals)
83 3.431798275933005
85 >>> # Two parallel streams with weighted flow
86 >>> log_removals = np.array([3, 5])
87 >>> flow_fractions = np.array([0.7, 0.3])
88 >>> calculate_parallel_log_removal(log_removals, flow_fractions)
89 3.153044674980176
91 See Also
92 --------
93 For systems in series, log removals would be summed directly.
94 """
95 # Convert log_removals to numpy array if it isn't already
96 log_removals = np.asarray(log_removals, dtype=float)
98 # Check if log_removals is empty
99 if len(log_removals) == 0:
100 msg = "At least one log removal value must be provided"
101 raise ValueError(msg)
103 # If flow_fractions is not provided, assume equal distribution
104 if flow_fractions is None:
105 # Calculate the number of parallel flows
106 n = len(log_removals)
107 # Create equal flow fractions
108 flow_fractions = np.full(n, 1.0 / n)
109 else:
110 # Convert flow_fractions to numpy array
111 flow_fractions = np.asarray(flow_fractions, dtype=float)
113 # Validate inputs
114 if len(log_removals) != len(flow_fractions):
115 msg = "log_removals and flow_fractions must have the same length"
116 raise ValueError(msg)
118 if not np.isclose(np.sum(flow_fractions), 1.0):
119 msg = "flow_fractions must sum to 1.0"
120 raise ValueError(msg)
122 # Convert log removal to decimal reduction values
123 decimal_reductions = 10 ** (-log_removals)
125 # Calculate weighted average decimal reduction
126 weighted_decimal_reduction = np.sum(flow_fractions * decimal_reductions)
128 # Convert back to log scale
129 return -np.log10(weighted_decimal_reduction)
132def gamma_pdf(r, rt_alpha, rt_beta, log_removal_rate):
133 """
134 Compute the probability density function (PDF) of log removal given a gamma distribution for the residence time.
136 gamma(rt_alpha, rt_beta) = gamma(apv_alpha, apv_beta / flow)
138 Parameters
139 ----------
140 r : array_like
141 Log removal values at which to compute the PDF.
142 rt_alpha : float
143 Shape parameter of the gamma distribution for residence time.
144 rt_beta : float
145 Scale parameter of the gamma distribution for residence time.
146 log_removal_rate : float
147 Coefficient for log removal calculation (R = log_removal_rate * log10(T)).
149 Returns
150 -------
151 pdf_values : ndarray
152 PDF values corresponding to the input r values.
153 """
154 # Compute the transformed PDF
155 t_values = 10 ** (r / log_removal_rate)
157 return (
158 (np.log(10) / (log_removal_rate * gamma(rt_alpha) * (rt_beta**rt_alpha)))
159 * (t_values**rt_alpha)
160 * np.exp(-t_values / rt_beta)
161 )
164def gamma_cdf(r, rt_alpha, rt_beta, log_removal_rate):
165 """
166 Compute the cumulative distribution function (CDF) of log removal given a gamma distribution for the residence time.
168 gamma(rt_alpha, rt_beta) = gamma(apv_alpha, apv_beta / flow)
170 Parameters
171 ----------
172 r : array_like
173 Log removal values at which to compute the CDF.
174 alpha : float
175 Shape parameter of the gamma distribution for residence time.
176 beta : float
177 Scale parameter of the gamma distribution for residence time.
178 log_removal_rate : float
179 Coefficient for log removal calculation (R = log_removal_rate * log10(T)).
181 Returns
182 -------
183 cdf_values : ndarray
184 CDF values corresponding to the input r values.
185 """
186 # Compute t values corresponding to r values
187 t_values = 10 ** (r / log_removal_rate)
189 # Use the gamma CDF directly
190 return stats.gamma.cdf(t_values, a=rt_alpha, scale=rt_beta)
193def gamma_mean(rt_alpha, rt_beta, log_removal_rate):
194 """
195 Compute the mean of the log removal distribution given a gamma distribution for the residence time.
197 gamma(rt_alpha, rt_beta) = gamma(apv_alpha, apv_beta / flow)
199 Parameters
200 ----------
201 rt_alpha : float
202 Shape parameter of the gamma distribution for residence time.
203 rt_beta : float
204 Scale parameter of the gamma distribution for residence time.
205 log_removal_rate : float
206 Coefficient for log removal calculation (R = log_removal_rate * log10(T)).
208 Returns
209 -------
210 mean : float
211 Mean value of the log removal distribution.
212 """
213 # Calculate E[R] = log_removal_rate * E[log10(T)]
214 # For gamma distribution: E[ln(T)] = digamma(alpha) + ln(beta_adjusted)
215 # Convert to log10: E[log10(T)] = E[ln(T)] / ln(10)
217 return (log_removal_rate / np.log(10)) * (digamma(rt_alpha) + np.log(rt_beta))
220def gamma_find_flow_for_target_mean(target_mean, apv_alpha, apv_beta, log_removal_rate):
221 """
222 Find the flow rate flow that produces a specified target mean log removal given a gamma distribution for the residence time.
224 gamma(rt_alpha, rt_beta) = gamma(apv_alpha, apv_beta / flow)
226 Parameters
227 ----------
228 target_mean : float
229 Target mean log removal value.
230 apv_alpha : float
231 Shape parameter of the gamma distribution for residence time.
232 apv_beta : float
233 Scale parameter of the gamma distribution for pore volume.
234 log_removal_rate : float
235 Coefficient for log removal calculation (R = log_removal_rate * log10(T)).
237 Returns
238 -------
239 flow : float
240 Flow rate that produces the target mean log removal.
242 Notes
243 -----
244 This function uses the analytical solution derived from the mean formula.
245 From E[R] = (log_removal_rate/ln(10)) * (digamma(alpha) + ln(beta) - ln(Q)),
246 we can solve for Q to get:
247 flow = beta * exp(ln(10)*target_mean/log_removal_rate - digamma(alpha))
248 """
249 # Rearranging the mean formula to solve for Q:
250 # target_mean = (log_removal_rate/ln(10)) * (digamma(alpha) + ln(beta) - ln(Q))
251 # ln(Q) = digamma(alpha) + ln(beta) - (ln(10)*target_mean/log_removal_rate)
252 # Q = beta * exp(-(ln(10)*target_mean/log_removal_rate - digamma(alpha)))
253 return apv_beta * np.exp(digamma(apv_alpha) - (np.log(10) * target_mean) / log_removal_rate)