diff --git a/.gitignore b/.gitignore index 061c6832..80d9a117 100755 --- a/.gitignore +++ b/.gitignore @@ -9,7 +9,7 @@ *.pyc *.egg *.egg-info -*.egg-info/* +pydfnworks/*.egg-info/* *.ipynb_checkpoints/* # Build artifacts diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py index 31ccba13..db8f786b 100644 --- a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py +++ b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py @@ -1,540 +1,183 @@ import numpy as np from scipy import special import mpmath as mp +import os from pydfnworks.general.logging import local_print_log +import pydfnworks.dfnGraph.tdrw_finite_roubinet as roubinet +import pydfnworks.dfnGraph.tdrw_finite_dentz as dentz +import pydfnworks.dfnGraph.tdrw_finite_annulus as annulus +import pydfnworks.dfnGraph.tdrw_finite_from_file as from_file -def get_fracture_segments(transfer_time, - fracture_length, - b, - velocity, - matrix_diffusivity, - matrix_porosity, - plim=0.01): - """ - This function calculates the maximum segment length that can be used in the TDRW model, following the protocol proposed by Roubinet et al., 2010. - - - Parameters - --------------------- - transfer_time : float - Time is takes for a particle to diffusive through the matrix across the fracture spacing. transfer_time = fracture_spacing**2 / (2 * matrix_diffusivity) - - fracture_length : float - Length of the current edge segment in the graph [m] - - b : float - Aperture of the current edge segment in the graph [m] - - velocity : float - particle velocity along the edge segment in the graph [m/s] - - matrix_diffusivity: float - Matrix Diffusivity [m^2/s] - - matrix_porosity: float - Matrix Porosity - - plim : float - Parameter used to break up the fracture - - Returns - --------------------- - segment_length : float - Length of the segment of the edge to compute the matrix diffusion time - - num_segments : int - Number of segments of length segment_length on the original edge/fracture. We ake the ceiling of the value. - - Notes - --------------------- - Roubinet use a plim value of 0.1. However, we are better able to match the results of Sudicky & Frind with a parameter values when plim = 0.01. - - """ - - # Compute the maximum segment length - segment_length = ((b * np.sqrt(transfer_time)) / - (matrix_porosity * np.sqrt(matrix_diffusivity)) - ) * special.erfcinv(1.0 - plim) * velocity - # if that's bigger than the edge size, just return the edge length as a single segment - if segment_length >= fracture_length: - return fracture_length, 1 - else: - # otherwise, get the number of segments of length segment_length on the edge - num_segments = int(np.ceil(fracture_length / segment_length)) - return segment_length, num_segments - - -def t_diff_unlimited(a, tf, xi): - """ - This function returns one diffusion time sample for a and tf assuming an unlimited matrix block size. This is used throughout the limited block size sampling technique. - - Parameters - ------------------- - a : double - Constant parameter describing retention in the matrix. a = (matrix_porosity*matrix_diffusion)/aperture_length - - tf : double - Advective travel time - - xi : float - value between [0,1) - - Returns - -------------- - t_diff : float - Total time diffusing the matrix - - Notescx - ------------- - For a random sample, sample xi from U[0,1). - - """ - return ((a * tf) / special.erfcinv(xi))**2 +# Valid finite matrix diffusion models that require fracture_spacing +_FINITE_MODELS_REQUIRING_SPACING = ("roubinet", "dentz", "annulus") -def transition_probability_cdf(t_min, t_max, frac_spacing, matrix_diffusivity, - num_pts): - """ - This function calculates the cummulative probability density that a particle changes fractures, given the time needed to reach a penetration depth and an associated time with each probabilty. - - Parameters - --------------- - t_min : float - Minumum value of diffusion time [s] - - t_max : float - Maximum value of diffusion time [s] - - frac_spacing : float - Spacing between fractures [m] - - matrix_diffusivity : float - Matrix Diffusivity value [m^2/s] - - num_pts : int - Number of points in the logspace array between t_min and t_max - - Returns - -------------- - times : np.array - Array of diffusion times - - prob_cdf : np.array - Array of cummulative probabilities. They only go to 0.5 - - - Notes - ------------ - Negative probabilities can sometimes appear due to numerical instabilities in our laplace inverse transform. These are removed in the section below marked CLEAN UP. We also ensure that our final distribution is monotonically increasing - - """ - local_print_log("--> Building transition probability cdf") - - times = np.logspace(np.log10(t_min), np.log10(t_max), num=num_pts) - prob_cdf = np.zeros(num_pts) - - ## Parameters for Roubinet function for the transfer probability function described - ## in Roubinet et al. WRR 2010. Equation number 5. - l1 = frac_spacing / 2 - l2 = -l1 - delta_l = l1 - l2 - - ## Roubinet et al. WRR 2010. Equation number 5. - ## s is the Laplace variable - laplace_trans_prob_function = lambda s: (mp.exp(l1*mp.sqrt(s/matrix_diffusivity))/s) \ - *((1.0 - mp.exp(-2*l2*mp.sqrt(s/matrix_diffusivity)))/(1.0 - mp.exp(2.0 *(delta_l)*mp.sqrt(s/matrix_diffusivity)))) - - for i, t in enumerate(times): - prob_cdf[i] = mp.invertlaplace(laplace_trans_prob_function, - t, - method='talbot', - dps=16, - degree=36) - - if prob_cdf[i] >= 0.5: - break - - # # CLEAN UP the solution. - # # Sometime small negative values will appear due to numerical instabilities. Remove these. - ind = np.asarray(np.where(prob_cdf >= 0)).flatten() - prob_cdf = prob_cdf[ind] - times = times[ind] - - # # Monotonically increasing - prob_cdf = np.maximum.accumulate(prob_cdf) - # unique values - prob_cdf, ind = np.unique(prob_cdf, return_index=True) - times = times[ind] - - return times, prob_cdf - - -def transfer_probabilities(b_min, - b_max, - tf_min, - tf_max, - matrix_porosity, - matrix_diffusivity, - frac_spacing, - eps=1e-16, - num_pts=100): - """ Returns the CDF of transfer probabilities and assocaited times. Lower and upper bounds are first estimated using the physical parameters of the system. The bounds are then tighted based on the returned probabilities wherein we find a range with probabilities greater than eps and 0.5 +def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model, tdrw_filename): + """ Check that the provided tdrw values are physical Parameters - ------------ - b_min : float - Minimum aperture in the network - - b_max : float - Maximum aperture in the network - - tf_min : float - Minimum advective travel time in the network - - tf_max : float - Maximum advective travel time in the network - - matrix_porosity: float - Matrix Porosity - - matrix_diffusivity : float - Matrix Diffusivity value [m^2/s] - - frac_spacing : float - Spacing between fractures [m] - - eps : float - Default - 1e-16 - - num_pts : int - Number of points in the logspace array between t_min and t_max - - Returns - ------------ - trans_prob : dictionary - Dictionary elements - times : np.array - Array of diffusion times - prob_cdf : np.array - Array of cummulative probabilities. They only go to 0.5 - - Notes - ------------ - The initial bounds are huge, don't worry. They get tighted up just fine. - - """ - # estimate lower bound - # get the minimum factor - a_min = (matrix_porosity * np.sqrt(matrix_diffusivity)) / b_max - # sample at the lowest value of tf, with sample ~ 0 - t_diff_lb = t_diff_unlimited(a_min, tf_min, eps) - # Below this value the inversse laplace transform does not convergence - if t_diff_lb < 1e-12: - local_print_log(f"lb too low {t_diff_lb:0.2e} changing to 1e-12", 'warning') - t_diff_lb = 1e-12 - - # estimate upper bound - # get the maximum factor - a_max = (matrix_porosity * np.sqrt(matrix_diffusivity)) / b_min - # sample at the largest value of tf, with sample ~ 1 - t_diff_ub = t_diff_unlimited(a_max, tf_max, 1 - eps) - # if t_diff_ub > 1e30: - # print(f"ub too high {t_diff_ub} changing to 1e30") - # t_diff_ub = 1e30 - - local_print_log( - f"--> Initial bounds for diffusion times. Min: {t_diff_lb:0.2e}, Max: {t_diff_ub:0.2e}" - ) - - # Compute the transition probabilities across the whole range of times. - times, trans_cdf = transition_probability_cdf(t_diff_lb, t_diff_ub, - frac_spacing, - matrix_diffusivity, num_pts) - - # Restricts the - # We end up with a lot of values we don't need, close to 0 probability. - # walk through the CDF and find the time corresponding to the - # first probability greater than epsilon - for i, val in enumerate(trans_cdf): - if val > eps: - t_diff_lb = times[i] - break - - # Likewise we end up with a lot of values we don't need, close to 0.5 probability. - # walk through the CDF and find the time corresponding to the - # first probability greater than 0.5 - epsilon - - for i, val in reversed(list(enumerate(trans_cdf))): - if val <= 0.5: - t_diff_ub = times[i] - break - - local_print_log( - f"--> Final bounds for diffusion times. Min: {t_diff_lb:0.2e}, Max: {t_diff_ub:0.2e}" - ) - - # Recompute the transition probabilities across the restricted range of times. - times, trans_cdf = transition_probability_cdf(t_diff_lb, t_diff_ub, - frac_spacing, - matrix_diffusivity, num_pts) - - #convert to as dictionary - trans_prob = {"times": times, "cdf": trans_cdf} - # print(times) - # print(trans_cdf) - return trans_prob - - -def segment_matrix_diffusion(trans_prob, matrix_porosity, matrix_diffusivity, - b, velocity, segment_length, num_segments): - """ Computes the time delay for a particle due to matrix diffusion on a given edge in the graph. The edge might already be broken into multiple segments depending on the parameters of the simulation. - - Parameters - -------------- - trans_prob : dictionary - Dictionary elements - times : np.array - Array of diffusion times - prob_cdf : np.array - Array of cummulative probabilities. They only go to 0.5 - - matrix_porosity: float - Matrix Porosity + ---------- + matrix_porosity : float + Porosity of the rock matrix [-] matrix_diffusivity : float - Matrix Diffusivity value [m^2/s] - - b : float - Aperture of the current edge segment in the graph [m] + Effective diffusivity of the rock matrix [m^2/s] - velocity : float - particle velocity along the edge segment in the graph [m/s] + fracture_spacing : float + Characteristic matrix block half-width / fracture spacing [m] + Required for finite matrix models. - segment_length : float - Length of the segment of the edge to compute the matrix diffusion time - - num_segments : int - Number of segments of length segment_length on the original edge/fracture. We ake the ceiling of the value. + tdrw_model : str + Name of the TDRW matrix diffusion model. Valid options: + 'infinite', 'roubinet', 'dentz', 'annulus', 'from_file'. + tdrw_filename : str or None + Path to CDF file, required when tdrw_model == 'from_file'. Returns - -------------- - t_diff : float - Total time delay due to retention via matrix diffusion - - Notes - ------------- + ------- None + Notes + ----- + Logs errors and warnings via local_print_log. Does not raise + exceptions; callers should treat logged errors as fatal. """ + local_print_log("\n* TDRW for Matrix Diffusion is ON.") + + # --- Matrix Porosity --- + if matrix_porosity is None: + local_print_log("Error: Requested TDRW but no value for matrix_porosity was provided.", "error") + elif not (0 <= matrix_porosity <= 1): + local_print_log( + f"Error: Requested TDRW but matrix_porosity={matrix_porosity} is outside the valid range [0, 1].", + "error" + ) + else: + local_print_log(f"--> Matrix porosity value: {matrix_porosity:.2f} [-]") - a = (matrix_porosity * np.sqrt(matrix_diffusivity)) / b - - segment_time = segment_length / velocity - prob_min = min(trans_prob['cdf']) - prob_max = max(trans_prob['cdf']) - t_diff = 0 - for i in range(num_segments): - # sample a diffusion time from the unlimited scenario - xi = np.random.uniform(low=0, high=1) - t_diff_seg = t_diff_unlimited(a, segment_time, xi) - - # compare that time with limited times - # if it's less than the minimum time, there is 0 probability for - # transfer and we accept that time - if t_diff_seg < trans_prob['times'][0]: - limited_probability = 0 - - # if it's greater than the maximum time, there is probability 1 for - # transfer and we resample - elif t_diff_seg > trans_prob['times'][-1]: - limited_probability = 1 - - # otherwise, we get a new probability using the CDF based on the sampled time. + # --- Matrix Diffusivity --- + if matrix_diffusivity is None: + local_print_log("Error: Requested TDRW but no value for matrix_diffusivity was provided.", "error") + elif matrix_diffusivity <= 0: + local_print_log(f"Error: Non-positive matrix_diffusivity={matrix_diffusivity}. Exiting program.", "error") + else: + local_print_log(f"--> Matrix diffusivity value: {matrix_diffusivity:.2e} [m^2/s]") + + # --- TDRW Model Selection --- + if tdrw_model is None: + local_print_log("Warning: No TDRW model specified; default behavior may apply.", "warning") + + elif tdrw_model == "infinite": + local_print_log("--> Using infinite matrix diffusion model (no fracture spacing required).") + + elif tdrw_model in _FINITE_MODELS_REQUIRING_SPACING: + local_print_log(f"--> Using {tdrw_model.capitalize()} model for finite matrix diffusion.") + # --- Fracture Spacing Check --- + if fracture_spacing is None: + local_print_log( + f"Error: TDRW model '{tdrw_model}' requires a fracture_spacing value, but none was provided.", + "error" + ) + elif fracture_spacing <= 0: + local_print_log(f"Error: Non-positive fracture_spacing={fracture_spacing}. Exiting program.", "error") else: - limited_probability = 2 * np.interp( - t_diff_seg / 2, trans_prob['times'], trans_prob['cdf']) - - if limited_probability > 0: - # Draw a random number from U[0,1] for each particle. - # Particles transfer to a new fracture if this random number - # is less than the transfer probability. - # this is the inverse CDF method part of the algorithm. - xi = np.random.uniform(low=0, high=1) - if xi < limited_probability: - xi = np.random.uniform(low=prob_min, high=prob_max) - t_diff_seg = np.interp(xi, trans_prob['cdf'], - trans_prob['times']) + local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e} [m]") + if tdrw_model == "annulus": + local_print_log("--> Annulus geometry: absorbing boundary at fracture-matrix interface, reflecting boundary at block interior.") - t_diff += t_diff_seg - return t_diff + elif tdrw_model == "from_file": + if tdrw_filename is None: + local_print_log("Error. TDRW Filename not provided.", "error") + elif not os.path.isfile(tdrw_filename): + local_print_log(f"Error. File not found: {tdrw_filename}", "error") + else: + local_print_log(f"--> Sampling matrix diffusion time from: {tdrw_filename}.") + else: + local_print_log( + f"Error: Unknown TDRW model '{tdrw_model}'. Valid options are 'infinite', 'roubinet', 'dentz', 'annulus', or 'from_file' (case-sensitive).", + "error" + ) + + local_print_log("") # clean trailing newline + + +def _set_up_limited_matrix_diffusion(G, + tdrw_model, + fracture_spacing, + matrix_porosity, + matrix_diffusivity, + eps=1e-16, + num_pts=100, + tdrw_filename=None): + """ Sets up transition probabilities for limited block size matrix diffusion -def get_aperture_and_time_limits(G): - """ Walks through edges on the graph and returns min and max of aperture and advective travel times - Parameters - --------------------- - G : networkX graph + ---------- + G : networkX graph Graph provided by graph_flow modules - Returns - ------------- - b_min : float - Minimum value of aperture on the network - - b_max : float - Maximmum value of aperture on the network - - t_min : float - Minimum value of advective travel time on the network - - t_max : float - Maximmum value of advective travel time on the network + tdrw_model : str + Name of the TDRW matrix diffusion model. Valid options: + 'roubinet', 'dentz', 'annulus', 'from_file'. - Notes - -------------------- - Could be improved with nx.get_edge_attributes, maybe. JDH + fracture_spacing : float + Characteristic matrix block half-width [m]. + Required for 'roubinet', 'dentz', and 'annulus' models. - - """ - # get b_min, b_max, t_min, t_max - local_print_log("--> Getting b and t limits") - b_min = None - b_max = None - t_min = None - t_max = None - for u, v, d in G.edges(data=True): - if b_min is None: - b_min = d['b'] - elif d['b'] < b_min: - b_min = d['b'] - - if b_max is None: - b_max = d['b'] - elif d['b'] > b_max: - b_max = d['b'] - - if t_min is None: - t_min = d['time'] - elif d['time'] < t_min: - t_min = d['time'] - - if t_max is None: - t_max = d['time'] - elif d['time'] > t_max: - t_max = d['time'] - - local_print_log(f"--> b-min: {b_min:0.2e}, b-max: {b_max:0.2e}") - local_print_log(f"--> t-min: {t_min:0.2e}, t-max: {t_max:0.2e}") - return b_min, b_max, t_min, t_max - - -def set_up_limited_matrix_diffusion(G, - frac_spacing, - matrix_porosity, - matrix_diffusivity, - eps=1e-16, - num_pts=100): - """ Sets up transition probabilities for limited block size matrix diffusion - - Parameters - --------------------- - G : networkX graph - Graph provided by graph_flow modules - - fracture_length : float - Length of the current edge segment in the graph [m] - - matrix_porosity: float - Matrix Porosity + matrix_porosity : float + Matrix porosity [-] matrix_diffusivity : float - Matrix Diffusivity value [m^2/s] + Matrix diffusivity [m^2/s] - eps : float - Default - 1e-16 + eps : float + Small parameter for Roubinet model. Default 1e-16. - num_pts : int - Number of points in the logspace array between t_min and t_max - + num_pts : int + Number of points in the logspace CDF table. Default 100. - Returns - ------------- - trans_prob : dictionary - Dictionary elements - times : np.array - Array of diffusion times - prob_cdf : np.array - Array of cummulative probabilities. They only go to 0.5 - - Notes - -------------------- - None - - """ - - b_min, b_max, tf_min, tf_max = get_aperture_and_time_limits(G) - trans_prob = transfer_probabilities(b_min, b_max, tf_min, tf_max, - matrix_porosity, matrix_diffusivity, - frac_spacing, eps, num_pts) - return trans_prob - - -def limited_matrix_diffusion(self, G): - """ Matrix diffusion with limited block size - - Parameters - ---------- - G : NetworkX graph - graph obtained from graph_flow + tdrw_filename : str or None + Path to CDF file, required when tdrw_model == 'from_file'. Returns ------- - None + transfer_time : np.ndarray or None + Array of dimensionless diffusion return times for inverse CDF lookup. + + trans_prob : np.ndarray or None + Array of cumulative probabilities corresponding to transfer_time. Notes - ----------- - All parameters are attached to the particle class + ----- + The returned arrays are used with np.interp to map uniform random + samples to return times during particle transport. See + limited_matrix_diffusion_* functions in each model module. """ + if tdrw_model == 'roubinet': + b_min, b_max, tf_min, tf_max = roubinet.get_aperture_and_time_limits(G) + trans_prob = roubinet.transfer_probabilities(b_min, b_max, tf_min, tf_max, + matrix_porosity, matrix_diffusivity, + fracture_spacing, eps, num_pts) + transfer_time = fracture_spacing**2 / (2 * matrix_diffusivity) - frac_length = G.edges[self.curr_node, self.next_node]['length'] - b = G.edges[self.curr_node, self.next_node]['b'] - velocity = G.edges[self.curr_node, self.next_node]['velocity'] - - segment_length, num_segments = get_fracture_segments( - self.transfer_time, frac_length, b, velocity, self.matrix_diffusivity, - self.matrix_porosity) - - self.delta_t_md = segment_matrix_diffusion(self.trans_prob, - self.matrix_porosity, - self.matrix_diffusivity, b, - velocity, segment_length, - num_segments) + elif tdrw_model == "dentz": + transfer_time, trans_prob = dentz._make_inverse_cdf_spline_for_times() + elif tdrw_model == "annulus": + # annulus model: absorbing wall at fracture-matrix interface, + # reflecting wall at block interior; eps fixed inside the module + transfer_time, trans_prob = annulus._make_inverse_cdf_annulus(num_samples=num_pts) -def unlimited_matrix_diffusion(self, G): - """ Matrix diffusion with unlimited block size + elif tdrw_model == "from_file": + transfer_time, trans_prob = from_file._load_finite_time_cdf(tdrw_filename) - Parameters - ---------- - G : NetworkX graph - graph obtained from graph_flow - - Returns - ---------- - None - - Notes - ----------- - All parameters are attached to the particle class - - """ + else: + trans_prob = None + transfer_time = None - b = G.edges[self.curr_node, self.next_node]['b'] - a_nondim = self.matrix_porosity * np.sqrt(self.matrix_diffusivity) / b - xi = np.random.uniform(size=1, low=0, high=1)[0] - self.delta_t_md = ((a_nondim * self.delta_t / special.erfcinv(xi))**2) + return transfer_time, trans_prob diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py index b0660795..3cd1a3f6 100644 --- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py +++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py @@ -13,32 +13,34 @@ # pydfnworks graph modules modules import pydfnworks.dfnGraph.particle_io as io -from pydfnworks.dfnGraph.graph_tdrw import set_up_limited_matrix_diffusion from pydfnworks.dfnGraph.particle_class import Particle from pydfnworks.general.logging import local_print_log +from pydfnworks.dfnGraph.graph_tdrw import _check_tdrw_params, _set_up_limited_matrix_diffusion +from pydfnworks.dfnGraph.graph_transport_setup_functions import _create_neighbor_list, _get_initial_posititions, _check_control_planes -def track_particle(data, verbose=False): + +def track_particle(args): """ Tracks a single particle through the graph - all input parameters are in the dictionary named data + all input parameters are packed into a single tuple (data, verbose) + for compatibility with multiprocessing.Pool.imap_unordered. Parameters ---------- - data : dict - Dictionary of parameters the includes particle_number, initial_position, - tdrw_flag, matrix_porosity, matrix_diffusivity, cp_flag, control_planes, - direction, G, and nbrs_dict. - - verbose : bool - Toggles verbosity + args : tuple + (data, verbose) where data is a dict containing particle_number, + initial_position, tdrw_flag, matrix_porosity, matrix_diffusivity, + cp_flag, control_planes, direction, G, and nbrs_dict; + and verbose is a bool toggling verbosity. Returns ------- particle : object - Particle will full trajectory + Particle with full trajectory """ + data, verbose = args if verbose: p = mp.current_process() _, cpu_id = p.name.split("-") @@ -47,12 +49,21 @@ def track_particle(data, verbose=False): f"--> Particle {data['particle_number']} is starting on worker {cpu_id}" ) - particle = Particle(data["particle_number"], data["initial_position"], - data["tdrw_flag"], data["matrix_porosity"], - data["matrix_diffusivity"], data["fracture_spacing"], - data["trans_prob"], data["transfer_time"], - data["cp_flag"], data["control_planes"], - data["direction"]) + try: + particle = Particle(particle_number = data["particle_number"], + ip = data["initial_position"], + tdrw_flag = data["tdrw_flag"], + tdrw_model = data["tdrw_model"], + matrix_porosity = data["matrix_porosity"], + matrix_diffusivity = data["matrix_diffusivity"], + fracture_spacing = data["fracture_spacing"], + transfer_time = data["transfer_time"], + trans_prob = data["trans_prob"], + cp_flag = data["cp_flag"], + control_planes = data["control_planes"], + direction = data["direction"]) + except Exception as e: + local_print_log(f"Issue initializing particle {e}", "warning") # # get current process information global nbrs_dict @@ -65,181 +76,6 @@ def track_particle(data, verbose=False): return particle -def get_initial_posititions(G, initial_positions, nparticles): - """ Distributes initial particle positions - - Parameters - ---------- - - G : NetworkX graph - obtained from graph_flow - - initial_positions : str - distribution of initial conditions. options are uniform and flux (flux-weighted) - - nparticles : int - requested number of particles - - Returns - ------- - ip : numpy array - array nparticles long. Each element is the initial position for each particle - - """ - - inlet_nodes = [v for v in nx.nodes(G) if G.nodes[v]['inletflag']] - cnt = len(inlet_nodes) - local_print_log(f"--> There are {cnt} inlet nodes") - if cnt == 0: - error = "Error. There are no nodes in the inlet.\nExiting" - local_print_log(error, 'error') - - # Uniform Distribution for particles - if initial_positions == "uniform": - local_print_log("--> Using uniform initial positions.") - ip = np.zeros(nparticles).astype(int) - n = int(np.ceil(nparticles / cnt)) - local_print_log(f"--> {n} particles will be placed at every inflow node.\n") - ## this could be cleaned up using clever indexing - inflow_idx = 0 - inflow_cnt = 0 - for i in range(nparticles): - ip[i] = inlet_nodes[inflow_idx] - inflow_cnt += 1 - if inflow_cnt >= n: - inflow_idx += 1 - inflow_cnt = 0 - - ## flux weighted initial positions for particles - elif initial_positions == "flux": - local_print_log("--> Using flux-weighted initial positions.\n") - flux = np.zeros(cnt) - for i, u in enumerate(inlet_nodes): - for v in G.successors(u): - flux[i] += G.edges[u, v]['flux'] - flux /= flux.sum() - flux_cnts = [np.ceil(nparticles * i) for i in flux] - nparticles = int(sum(flux_cnts)) - ip = np.zeros(nparticles).astype(int) - ## Populate ip with Flux Cnts - ## this could be cleaned up using clever indexing - inflow_idx = 0 - inflow_cnt = 0 - for i in range(nparticles): - ip[i] = inlet_nodes[inflow_idx] - inflow_cnt += 1 - if inflow_cnt >= flux_cnts[inflow_idx]: - inflow_idx += 1 - inflow_cnt = 0 - - # Throw error if unknown initial position is provided - else: - error = f"Error. Unknown initial_positions input {initial_positions}. Options are uniform or flux \n" - local_print_log(error, 'error') - - return ip, nparticles - - -def create_neighbor_list(G): - """ Create a list of downstream neighbor vertices for every vertex on NetworkX graph obtained after running graph_flow - - Parameters - ---------- - G: NetworkX graph - Directed Graph obtained from output of graph_flow - - Returns - ------- - dict : nested dictionary. - - Notes - ----- - dict[n]['child'] is a list of vertices downstream to vertex n - dict[n]['prob'] is a list of probabilities for choosing a downstream node for vertex n - """ - - nbrs_dict = {} - - for u in nx.nodes(G): - - if G.nodes[u]['outletflag']: - continue - - node_list = [] - prob_list = [] - nbrs_dict[u] = {} - - for v in G.successors(u): - node_list.append(v) - prob_list.append(G.edges[u, v]['vol_flow_rate']) - - if node_list: - nbrs_dict[u]['child'] = node_list - nbrs_dict[u]['prob'] = np.asarray(prob_list) / sum(prob_list) - else: - nbrs_dict[u]['child'] = None - nbrs_dict[u]['prob'] = None - - return nbrs_dict - - -def check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing): - """ Check that the provided tdrw values are physiscal - - - Parameters - ---------- - G: NetworkX graph - Directed Graph obtained from output of graph_flow - - Returns - ------- - dict : nested dictionary. - - Notes - ----- - dict[n]['child'] is a list of vertices downstream to vertex n - dict[n]['prob'] is a list of probabilities for choosing a downstream node for vertex n - - """ - - if matrix_porosity is None: - error = f"Error. Requested TDRW but no value for matrix_porosity was provided\n" - local_print_log(error, 'error') - elif matrix_porosity < 0 or matrix_porosity > 1: - error = f"Error. Requested TDRW but value for matrix_porosity provided is outside of [0,1]. Value provided {matrix_porosity}\n" - local_print_log(error, 'error') - if matrix_diffusivity is None: - error = f"Error. Requested TDRW but no value for matrix_diffusivity was provided\n" - local_print_log(error, 'error') - - if fracture_spacing is not None: - if fracture_spacing <= 0: - error = f"Error. Non-positive value for fracture_spacing was provided.\nValue {fracture_spacing}\nExiting program" - local_print_log(error, 'error') - - -def check_control_planes(control_planes, direction): - control_plane_flag = False - if not type(control_planes) is list: - error = f"Error. provided controls planes are not a list\n" - local_print_log(error, 'error') - else: - # add None to indicate the end of the control plane list - control_plane_flag = True - - if direction is None: - error = f"Error. Primary direction not provided. Required for control planes\n" - local_print_log(error, 'error') - elif direction not in ['x', 'y', 'z']: - error = f"Error. Primary direction is not known. Acceptable values are x,y, and z\n" - local_print_log(error, 'error') - - local_print_log(f"--> Control Planes: {control_planes}") - local_print_log(f"--> Direction: {direction}") - return control_plane_flag - - def run_graph_transport(self, G, nparticles, @@ -249,12 +85,15 @@ def run_graph_transport(self, initial_positions="uniform", dump_traj=False, tdrw_flag=False, + tdrw_model = 'infinite', matrix_porosity=None, matrix_diffusivity=None, fracture_spacing=None, + tdrw_filename = None, control_planes=None, direction=None, - cp_filename='control_planes'): + cp_filename='control_planes', + verbose = False): """ Run particle tracking on the given NetworkX graph Parameters @@ -320,57 +159,67 @@ def run_graph_transport(self, self.print_log(error, 'error') self.print_log("--> Running Graph Particle Tracking") - - # Check parameters for TDRW - if tdrw_flag: - check_tdrw_params(matrix_porosity, matrix_diffusivity, - fracture_spacing) - self.print_log( - f"--> Running particle transport with TDRW.\n--> Matrix porosity {matrix_porosity}.\n--> Matrix Diffusivity {matrix_diffusivity} m^2/s" - ) if control_planes is None: control_plane_flag = False else: - control_plane_flag = check_control_planes( + control_plane_flag = _check_control_planes( control_planes=control_planes, direction=direction) self.print_log(f"--> Control Plane Flag {control_plane_flag}") self.print_log("--> Creating downstream neighbor list") global nbrs_dict - nbrs_dict = create_neighbor_list(G) + nbrs_dict = _create_neighbor_list(G) self.print_log("--> Getting initial Conditions") - ip, nparticles = get_initial_posititions(G, initial_positions, nparticles) - + ip, nparticles = _get_initial_posititions(G, initial_positions, nparticles) self.print_log(f"--> Starting particle tracking for {nparticles} particles") if dump_traj: self.print_log(f"--> Writing trajectory information to file") - if fracture_spacing is not None: - self.print_log(f"--> Using limited matrix block size for TDRW") - self.print_log(f"--> Fracture spacing {fracture_spacing:0.2e} [m]") - trans_prob = set_up_limited_matrix_diffusion(G, fracture_spacing, - matrix_porosity, - matrix_diffusivity) - # This doesn't change for the system. - # Transfer time diffusing between fracture blocks - transfer_time = fracture_spacing**2 / (2 * matrix_diffusivity) + # Check parameters for TDRW + if tdrw_flag: + _check_tdrw_params(matrix_porosity, matrix_diffusivity, + fracture_spacing, tdrw_model, tdrw_filename) + + if tdrw_model in ("roubinet", "dentz", "annulus", "from_file"): + self.print_log(f"--> Using limited matrix block size for TDRW") + self.print_log(f"--> Fracture spacing {fracture_spacing:0.2e} [m]") + transfer_time, trans_prob = _set_up_limited_matrix_diffusion(G = G, + tdrw_model = tdrw_model, + fracture_spacing = fracture_spacing, + matrix_porosity = matrix_porosity, + matrix_diffusivity = matrix_diffusivity, + tdrw_filename = tdrw_filename) + # np.savetxt("my_cdf.dat", np.c_[transfer_time,trans_prob]) + else: + trans_prob = None + transfer_time = None else: trans_prob = None transfer_time = None - ## main loop + if self.ncpu == 1: + self.print_log("--> Running in Serial") tic = timeit.default_timer() particles = [] for i in range(nparticles): if i % 1000 == 0: self.print_log(f"--> Starting particle {i} out of {nparticles}") - particle = Particle(i, ip[i], tdrw_flag, matrix_porosity, - matrix_diffusivity, fracture_spacing, - trans_prob, transfer_time, control_plane_flag, - control_planes, direction) + particle = Particle(particle_number = i, + ip = ip[i], + tdrw_flag = tdrw_flag, + tdrw_model = tdrw_model, + matrix_porosity = matrix_porosity, + matrix_diffusivity = matrix_diffusivity, + fracture_spacing = fracture_spacing, + transfer_time = transfer_time, + trans_prob = trans_prob, + cp_flag = control_plane_flag, + control_planes = control_planes, + direction = direction) + particle.track(G, nbrs_dict) particles.append(particle) @@ -388,23 +237,16 @@ def run_graph_transport(self, io.dump_trajectories(particles, 1) if self.ncpu > 1: - self.print_log(f"--> Using {self.ncpu} processors") - ## Prepare input data - inputs = [] - - tic = timeit.default_timer() - pool = mp.Pool(min(self.ncpu, nparticles)) - - particles = [] - - def gather_output(output): - particles.append(output) + self.print_log(f"--> Building initial data") + # build input data list for all particles + all_data = [] for i in range(nparticles): data = {} data["particle_number"] = i data["initial_position"] = ip[i] data["tdrw_flag"] = tdrw_flag + data["tdrw_model"] = tdrw_model data["matrix_porosity"] = matrix_porosity data["matrix_diffusivity"] = matrix_diffusivity data["fracture_spacing"] = fracture_spacing @@ -413,13 +255,27 @@ def gather_output(output): data["cp_flag"] = control_plane_flag data["control_planes"] = control_planes data["direction"] = direction - pool.apply_async(track_particle, - args=(data, ), - callback=gather_output) + all_data.append((data, verbose)) - pool.close() - pool.join() - pool.terminate() + tic = timeit.default_timer() + particles = [] + chunksize = max(1, nparticles // (4 * self.ncpu)) + + self.print_log("--> Starting main loop") + self.print_log(f"--> Using {self.ncpu} processors") + log_interval = max(1,int(nparticles * 0.025)) + + with mp.Pool(min(self.ncpu, nparticles)) as pool: + completed = 0 + next_log = log_interval + for i, particle in enumerate(pool.imap_unordered( + track_particle, all_data, chunksize=chunksize)): + particles.append(particle) + completed += 1 + if completed >= next_log: + percentage = 100 * completed/nparticles + self.print_log(f"--> Completed {completed} out of {nparticles} particles\t({percentage:0.2f}%)") + next_log += log_interval elapsed = timeit.default_timer() - tic self.print_log( diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport_setup_functions.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport_setup_functions.py new file mode 100644 index 00000000..d3158c72 --- /dev/null +++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport_setup_functions.py @@ -0,0 +1,183 @@ + +import networkx as nx +import numpy as np + +from pydfnworks.general.logging import local_print_log + +def _get_initial_posititions(G, initial_positions, nparticles): + """ Distributes initial particle positions + + Parameters + ---------- + + G : NetworkX graph + obtained from graph_flow + + initial_positions : str + distribution of initial conditions. options are uniform and flux (flux-weighted) + + nparticles : int + requested number of particles + + Returns + ------- + ip : numpy array + array nparticles long. Each element is the initial position for each particle + + """ + + inlet_nodes = [v for v in nx.nodes(G) if G.nodes[v]['inletflag']] + cnt = len(inlet_nodes) + local_print_log(f"--> There are {cnt} inlet nodes") + if cnt == 0: + error = "Error. There are no nodes in the inlet.\nExiting" + local_print_log(error, 'error') + + # Uniform Distribution for particles + if initial_positions == "uniform": + local_print_log("--> Using uniform initial positions.") + ip = np.zeros(nparticles).astype(int) + n = int(np.ceil(nparticles / cnt)) + local_print_log(f"--> {n} particles will be placed at every inflow node.\n") + ## this could be cleaned up using clever indexing + inflow_idx = 0 + inflow_cnt = 0 + for i in range(nparticles): + ip[i] = inlet_nodes[inflow_idx] + inflow_cnt += 1 + if inflow_cnt >= n: + inflow_idx += 1 + inflow_cnt = 0 + + ## flux weighted initial positions for particles + elif initial_positions == "flux": + local_print_log("--> Using flux-weighted initial positions.\n") + flux = np.zeros(cnt) + for i, u in enumerate(inlet_nodes): + for v in G.successors(u): + flux[i] += G.edges[u, v]['flux'] + flux /= flux.sum() + flux_cnts = [np.ceil(nparticles * i) for i in flux] + nparticles = int(sum(flux_cnts)) + ip = np.zeros(nparticles).astype(int) + ## Populate ip with Flux Cnts + ## this could be cleaned up using clever indexing + inflow_idx = 0 + inflow_cnt = 0 + for i in range(nparticles): + ip[i] = inlet_nodes[inflow_idx] + inflow_cnt += 1 + if inflow_cnt >= flux_cnts[inflow_idx]: + inflow_idx += 1 + inflow_cnt = 0 + + # Throw error if unknown initial position is provided + else: + error = f"Error. Unknown initial_positions input {initial_positions}. Options are uniform or flux \n" + local_print_log(error, 'error') + + return ip, nparticles + + +def _create_neighbor_list(G): + """ Create a list of downstream neighbor vertices for every vertex on NetworkX graph obtained after running graph_flow + + Parameters + ---------- + G: NetworkX graph + Directed Graph obtained from output of graph_flow + + Returns + ------- + dict : nested dictionary. + + Notes + ----- + dict[n]['child'] is a list of vertices downstream to vertex n + dict[n]['prob'] is a list of probabilities for choosing a downstream node for vertex n + """ + + nbrs_dict = {} + + for u in nx.nodes(G): + + if G.nodes[u]['outletflag']: + continue + + node_list = [] + prob_list = [] + nbrs_dict[u] = {} + + for v in G.successors(u): + node_list.append(v) + prob_list.append(G.edges[u, v]['vol_flow_rate']) + + if node_list: + nbrs_dict[u]['child'] = node_list + nbrs_dict[u]['prob'] = np.asarray(prob_list) / sum(prob_list) + else: + nbrs_dict[u]['child'] = None + nbrs_dict[u]['prob'] = None + + return nbrs_dict + +def _check_control_planes(control_planes, direction): + """ + Validate the control plane configuration and primary direction for a simulation. + + This function performs basic validation on the provided `control_planes` and `direction` + inputs used to define spatial control planes in a model or simulation. It logs errors + and diagnostic messages through `local_print_log` but does not raise exceptions directly. + The function returns a flag indicating whether the control planes list passed basic + validation checks. + + Parameters + ---------- + control_planes : list + A list defining the positions or indices of control planes along the specified + primary direction. Must be of type `list`. If validation succeeds, a flag is set + to indicate that control planes are active. + direction : str + The primary coordinate direction associated with the control planes. Must be one + of `'x'`, `'y'`, or `'z'`. Used to orient control plane processing. + + Returns + ------- + bool + `True` if the provided `control_planes` input is a list and passes validation, + otherwise `False`. + + Logging Behavior + ---------------- + Uses `local_print_log` to record: + * Errors if inputs are missing or invalid. + * Informational messages showing the provided control planes and direction. + + Notes + ----- + - The function does not modify the `control_planes` list. + - It is expected that higher-level routines handle program termination or recovery + in response to logged errors. + - The inclusion of `None` to mark the end of the control plane list should be + handled outside this function, if applicable. + """ + control_plane_flag = False + + # Validate control_planes + if not isinstance(control_planes, list): + local_print_log("Error: Provided control planes are not a list.", "error") + else: + control_plane_flag = True + + # Validate direction + allowed_dirs = {"x", "y", "z"} + if direction is None: + local_print_log("Error: Primary direction not provided. Required for control planes.", "error") + elif direction not in allowed_dirs: + local_print_log("Error: Primary direction is not known. Acceptable values are x, y, and z.", "error") + + # Diagnostics + local_print_log(f"Control planes: {control_planes}") + local_print_log(f"Direction: {direction}") + + return control_plane_flag \ No newline at end of file diff --git a/pydfnworks/pydfnworks/dfnGraph/particle_class.py b/pydfnworks/pydfnworks/dfnGraph/particle_class.py index 43bdbf98..4c10bbea 100644 --- a/pydfnworks/pydfnworks/dfnGraph/particle_class.py +++ b/pydfnworks/pydfnworks/dfnGraph/particle_class.py @@ -15,11 +15,16 @@ class Particle(): * frac_seq : Dictionary, contains information about fractures through which the particle went ''' - from pydfnworks.dfnGraph.graph_tdrw import unlimited_matrix_diffusion, limited_matrix_diffusion - - def __init__(self, particle_number, ip, tdrw_flag, matrix_porosity, - matrix_diffusivity, fracture_spacing, trans_prob, - transfer_time, cp_flag, control_planes, direction): + from pydfnworks.dfnGraph.tdrw_infinite import unlimited_matrix_diffusion + from pydfnworks.dfnGraph.tdrw_finite_dentz import limited_matrix_diffusion_dentz + from pydfnworks.dfnGraph.tdrw_finite_roubinet import limited_matrix_diffusion_roubinet + from pydfnworks.dfnGraph.tdrw_finite_annulus import limited_matrix_diffusion_annulus + + def __init__(self, particle_number, ip, tdrw_flag, tdrw_model, + matrix_porosity, + matrix_diffusivity, fracture_spacing, transfer_time, + trans_prob, cp_flag, control_planes, direction): + self.particle_number = particle_number self.ip = ip self.curr_node = ip @@ -38,7 +43,9 @@ def __init__(self, particle_number, ip, tdrw_flag, matrix_porosity, self.tdrw_flag = tdrw_flag self.matrix_porosity = matrix_porosity self.matrix_diffusivity = matrix_diffusivity + self.tau_D = None self.fracture_spacing = fracture_spacing + self.tdrw_model = tdrw_model self.trans_prob = trans_prob self.transfer_time = transfer_time self.cp_flag = cp_flag @@ -51,15 +58,20 @@ def __init__(self, particle_number, ip, tdrw_flag, matrix_porosity, self.cp_adv_time = [] self.cp_tdrw_time = [] self.cp_pathline_length = [] - # self.cp_x1 = [] - # self.cp_x2 = [] - self.velocity = [] self.lengths = [] self.times = [] self.coords = [] + if tdrw_model == "dentz": + self.tau_D = self.fracture_spacing**2 / (self.matrix_diffusivity) + + elif tdrw_model == "annulus": + # tau_D = tau1 = r1^2 / D where r1 = fracture_spacing is the outer + # (reflecting) radius of the cylindrical matrix block. Matches the + # tau1 = 1e0 reference timescale in the Laplace-domain solution. + self.tau_D = self.fracture_spacing**2 / (self.matrix_diffusivity) def initalize(self,G): self.frac = (G.nodes[self.curr_node]['frac'][1]) @@ -172,21 +184,12 @@ def cross_control_plane(self, G): l0 = self.interpolate_time(x0, l1, l2, x1, x2) self.cp_pathline_length.append(l0) - ##print(l1,l2,l0) - - # print(f"control plane: {x0:0.2f}, x1: {x1:0.2f}, x2:{x2:0.2f}, t1: {t1:0.2e}. t2: {t2:0.2e}, tau: {tau:0.2e}") if tau < 0: - self.print_log( f"control plane: {x0:0.2f}, x1: {x1:0.2f}, x2:{x2:0.2f}, t1: {t1:0.2e}. t2: {t2:0.2e}, tau: {tau:0.2e}" ) error = "Error. Interpolated negative travel time." self.print_log(error,'error') - # print(f"--> crossed control plane at {control_planes[cp_index]} {direction} at time {tau}") - # self.cp_adv_time.append(tau) - # self.cp_pathline_length.append(l0) - # self.cp_x1.append(x1) - # self.cp_x2.append(x2) if self.tdrw_flag: t1 = self.total_time @@ -197,7 +200,7 @@ def cross_control_plane(self, G): self.cp_tdrw_time.append(tau) self.cp_index += 1 - # if we're crossed all the control planes, turn off cp flag for this particle + # if we've crossed all the control planes, turn off cp flag for this particle if self.cp_index >= len(self.control_planes): self.cp_flag = False break @@ -228,7 +231,6 @@ def update(self): self.lengths.append(self.delta_l) self.times.append(self.delta_t) self.coords.append(self.curr_coords) - # self.fracs.append(self.frac) def cleanup_frac_seq(self): """ Cleanup @@ -269,15 +271,18 @@ def track(self, G, nbrs_dict): while not self.exit_flag: self.advect(G, nbrs_dict) if self.exit_flag: - # self.update() self.cleanup_frac_seq() break if self.tdrw_flag: - if self.fracture_spacing is None: + if self.tdrw_model == "roubinet": + self.limited_matrix_diffusion_roubinet(G) + elif self.tdrw_model == "dentz": + self.limited_matrix_diffusion_dentz(G) + elif self.tdrw_model == "annulus": + self.limited_matrix_diffusion_annulus(G) + elif self.tdrw_model == "infinite": self.unlimited_matrix_diffusion(G) - else: - self.limited_matrix_diffusion(G) if self.cp_flag: self.cross_control_plane(G) diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py new file mode 100644 index 00000000..abb14a93 --- /dev/null +++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py @@ -0,0 +1,176 @@ +import numpy as np +from scipy.special import ive, kve, iv, kv +from math import factorial + + +def _stehfest_coefficients(N=16): + # Compute Stehfest coefficients V_i for i = 1 ... N. N must be even. + if N % 2 != 0: + raise ValueError("Stehfest N must be even.") + M = N // 2 + V = np.zeros(N) + for i in range(1, N + 1): + total = 0.0 + for k in range(int((i + 1) // 2), min(i, M) + 1): + num = (k ** M) * factorial(2 * k) + den = (factorial(M - k) * factorial(k) * factorial(k - 1) + * factorial(i - k) * factorial(2 * k - i)) + total += num / den + V[i - 1] = ((-1) ** (i + M)) * total + return V + + +def Psi_star_annulus(s, eps, tau0, tau1): + # Laplace transform of the CDF for diffusion return time in a cylindrical annulus + # + # Geometry: cylindrical annulus with inner radius r0 (absorbing, fracture-matrix + # interface) and outer radius r1 (reflecting, block interior). + # Particle released at r' = r0*(1+eps), just inside the absorbing wall. + # + # Parameters: + # eps = (r' - r0) / r1 dimensionless release position + # tau0 = r0^2 / D inner-radius diffusion timescale + # tau1 = r1^2 / D outer-radius diffusion timescale + # + # Branch switch at |s*tau1| = 100, matching MATLAB cdfLaplace.m: + # small branch: unscaled iv/kv (safe for |s*tau1| <= 100) + # large branch: scaled ive/kve with exp prefactors (avoids overflow) + # + # ive(n,x) = exp(-x)*In(x) -> MATLAB besseli(n,x,1) + # kve(n,x) = exp(x) *Kn(x) -> MATLAB besselk(n,x,1) + # + # The branch switch is essential for correct small-t behavior. + # F_cdf(s) = F_psi(s) / s + + q0 = np.sqrt(s * tau0) + q0_eps = (1.0 + eps) * q0 + q1 = np.sqrt(s * tau1) + + large = s * tau1 > 100 + + # large-argument branch + I0e_eps = ive(0, q0_eps); K0e_eps = kve(0, q0_eps) + I0e_0 = ive(0, q0); K0e_0 = kve(0, q0) + I1e_1 = ive(1, q1); K1e_1 = kve(1, q1) + + exp1 = np.exp(np.clip((2 + eps) * q0 - 2 * q1, -500, 500)) + exp2 = np.exp(np.clip(-eps * q0, -500, 500)) + expd = np.exp(np.clip(2 * q0 - 2 * q1, -500, 500)) + + num_large = exp1 * I0e_eps * K1e_1 + exp2 * K0e_eps * I1e_1 + den_large = expd * I0e_0 * K1e_1 + K0e_0 * I1e_1 + + # small-argument branch + # errstate suppresses overflow warnings from elements that will be + # discarded by np.where (both branches are always evaluated) + with np.errstate(invalid='ignore', over='ignore'): + num_small = iv(0, q0_eps) * kv(1, q1) + kv(0, q0_eps) * iv(1, q1) + den_small = iv(0, q0) * kv(1, q1) + kv(0, q0) * iv(1, q1) + + num = np.where(large, num_large, num_small) + den = np.where(large, den_large, den_small) + + return num / (den * s) + + +def Psi_pdf_star_annulus(s, eps, tau0, tau1): + # Laplace transform of the first-passage time PDF + # psi*(s) = s * F_cdf(s) + return Psi_star_annulus(s, eps, tau0, tau1) * s + + +def _stehfest_invert(F, t_arr, V, **kwargs): + # Invert F at times t_arr using the Stehfest algorithm + # f(t) ≈ (ln2/t) * sum_i V_i * F(i*ln2/t) + ln2 = np.log(2.0) + f = np.zeros(len(t_arr)) + for idx, ti in enumerate(t_arr): + s_vals = np.array([(i + 1) * ln2 / ti for i in range(len(V))]) + f[idx] = (ln2 / ti) * np.dot(V, F(s_vals, **kwargs)) + return f + + +def _make_inverse_cdf_annulus(num_samples=100, eps=1e-2, tau0=1e-4, tau1=1e3, stehfest_n=16): + # Precompute the inverse CDF table for cylindrical annulus return-time sampling + # + # Parameters match Marco Dentz's InvLaplace.m defaults: + # eps = 1e-2 dimensionless release position (r'-r0)/r1 + # tau0 = 1e-4 inner-radius timescale r0^2/D + # tau1 = 1e3 outer-radius timescale r1^2/D + # + # Time range starts at 1e-10 to ensure Stehfest samples large enough s + # for correct early-time behavior (requires s >> (1/eps)^2 / tau0). + # + # Returns (times, cdf_vals) arrays for use with np.interp + + V = _stehfest_coefficients(stehfest_n) + times = np.logspace(-10, 2, num_samples) + + cdf_vals = _stehfest_invert(Psi_star_annulus, times, V, + eps=eps, tau0=tau0, tau1=tau1) + + # sort by cdf value to build the inverse CDF (cdf -> time) lookup + order = np.argsort(cdf_vals) + cdf_sorted = np.clip(cdf_vals[order], 0, 1) + t_sorted = np.maximum(times[order], 0) + + # remove duplicate cdf values to ensure monotone interpolation + cdf_unique, idx = np.unique(cdf_sorted, return_index=True) + t_unique = np.asarray(t_sorted[idx]) + + return t_unique, cdf_unique + + +def limited_matrix_diffusion_annulus(self, G): + """ Matrix diffusion with finite cylindrical annular block geometry + + Parameters + ---------- + G : NetworkX graph + graph obtained from graph_flow + + Returns + ------- + None + + Notes + ----- + Samples diffusion return times from a finite cylindrical annular + matrix block. The inner radius r0 is the absorbing fracture-matrix + interface; the outer radius r1 is the reflecting block interior. + A particle is released at r' = r0*(1 + eps), just inside the + absorbing wall. + + The Laplace-domain solution uses modified Bessel functions I0, I1, + K0, K1. Two numerical branches are used matching MATLAB cdfLaplace.m: + unscaled Bessels for |s*tau1| <= 100, and scaled Bessels (ive, kve) + with explicit exponential prefactors for |s*tau1| > 100. + + tau1 = r1^2 / D is stored on the particle as self.tau_D and used + to rescale dimensionless sampled return times to physical times [s]. + + The inverse CDF table (self.transfer_time, self.trans_prob) must + be precomputed via _make_inverse_cdf_annulus and attached to the + particle before transport begins. + + All other parameters are attached to the particle class. + """ + + eps = 1e-2 + b = G.edges[self.curr_node, self.next_node]['b'] + + # trapping rate into the cylindrical matrix block + gamma = (2 * self.matrix_porosity * self.matrix_diffusivity) / (b * eps * self.fracture_spacing) + + # average number of trapping events during this advective step + average_number_of_trapping_events = self.delta_t * gamma + + # sample number of trapping events from Poisson distribution + n = np.random.poisson(average_number_of_trapping_events) + + # sample uniform random variables and map to return times via inverse CDF + # self.tau_D = r1^2 / D rescales dimensionless times to physical times + xi = np.random.uniform(size=n) + tmp = self.tau_D * np.interp(xi, self.trans_prob, self.transfer_time) + + self.delta_t_md = tmp.sum() \ No newline at end of file diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py new file mode 100644 index 00000000..621dbf12 --- /dev/null +++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py @@ -0,0 +1,114 @@ +import numpy as np +from math import factorial + + +def _stehfest_coefficients(N=16): + # Compute Stehfest coefficients V_i for i = 1 ... N. N must be even. + if N % 2 != 0: + raise ValueError("Stehfest N must be even.") + M = N // 2 + V = np.zeros(N) + for i in range(1, N + 1): + total = 0.0 + for k in range(int((i + 1) // 2), min(i, M) + 1): + num = (k ** M) * factorial(2 * k) + den = (factorial(M - k) * factorial(k) * factorial(k - 1) + * factorial(i - k) * factorial(2 * k - i)) + total += num / den + V[i - 1] = ((-1) ** (i + M)) * total + return V + + +def Psi_star(s, eps): + # Laplace transform of the slab CDF Psi*_eps(s) + # cosh((1-eps)*sqrt(s)) / (s * cosh(sqrt(s))) + # + # Scaled form to avoid cosh overflow at large s: + # cosh((1-eps)*q) / cosh(q) = exp(-eps*q) * (1 + exp(-2*(1-eps)*q)) / (1 + exp(-2*q)) + # where q = sqrt(s). Both exp terms decay for large q so no overflow. + q = np.sqrt(s) + num = np.exp(-eps * q) * (1.0 + np.exp(-2.0 * (1.0 - eps) * q)) + den = 1.0 + np.exp(-2.0 * q) + return num / (den * s) + + +def Psi_pdf_star(s, eps): + # Laplace transform of the slab first-passage time PDF + # psi*(s) = s * Psi*(s) + return Psi_star(s, eps) * s + + +def _stehfest_invert(F, t_arr, V, **kwargs): + # Invert F at times t_arr using the Stehfest algorithm + # f(t) ≈ (ln2/t) * sum_i V_i * F(i*ln2/t) + ln2 = np.log(2.0) + f = np.zeros(len(t_arr)) + for idx, ti in enumerate(t_arr): + s_vals = np.array([(i + 1) * ln2 / ti for i in range(len(V))]) + f[idx] = (ln2 / ti) * np.dot(V, F(s_vals, **kwargs)) + return f + + +def _make_inverse_cdf_spline_for_times(num_samples=100, eps=1e-4, stehfest_n=16): + # Precompute the inverse CDF table for slab (Dentz) return-time sampling + # + # Replaces the original mpmath.invertlaplace implementation with + # Stehfest inversion in double precision -- orders of magnitude faster. + # + # Scaled cosh formulation avoids overflow for all s values. + # t_min=1e-10 ensures Stehfest samples s >> 1/eps^2 for correct + # early-time behavior (eps=1e-4 requires s ~ 1e8, t ~ ln2/1e8 ~ 1e-9). + # + # Returns (times, cdf_vals) arrays for use with np.interp + + V = _stehfest_coefficients(stehfest_n) + times = np.logspace(-10, 3, num_samples) + + cdf_vals = _stehfest_invert(Psi_star, times, V, eps=eps) + + # sort by cdf value to build the inverse CDF (cdf -> time) lookup + order = np.argsort(cdf_vals) + cdf_sorted = np.clip(cdf_vals[order], 0, 1) + t_sorted = np.maximum(times[order], 0) + + # remove duplicate cdf values to ensure monotone interpolation + cdf_unique, idx = np.unique(cdf_sorted, return_index=True) + t_unique = np.asarray(t_sorted[idx]) + + return t_unique, cdf_unique + + +def limited_matrix_diffusion_dentz(self, G): + """ Matrix diffusion with limited block size + + Parameters + ---------- + G : NetworkX graph + graph obtained from graph_flow + + Returns + ------- + None + + Notes + ----------- + All parameters are attached to the particle class + """ + + eps = 1e-4 + b = G.edges[self.curr_node, self.next_node]['b'] + + # trapping rate + gamma = (2 * self.matrix_porosity * self.matrix_diffusivity) / (b * eps * self.fracture_spacing) + + # average number of trapping events during this advective step + average_number_of_trapping_events = self.delta_t * gamma + + # sample number of trapping events from Poisson distribution + n = np.random.poisson(average_number_of_trapping_events) + + # sample uniform random variables and map to return times via inverse CDF + xi = np.random.uniform(size=n) + tmp = self.tau_D * np.interp(xi, self.trans_prob, self.transfer_time) + + self.delta_t_md = tmp.sum() \ No newline at end of file diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_from_file.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_from_file.py new file mode 100644 index 00000000..ecc7243a --- /dev/null +++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_from_file.py @@ -0,0 +1,73 @@ + + +from scipy.interpolate import PchipInterpolator +import numpy as np +import matplotlib.pyplot as plt +import mpmath as mp +import os +from pydfnworks.general.logging import local_print_log + +def _load_finite_time_cdf(tdrw_filename): + + data = np.genfromtxt(tdrw_filename) + + if data.ndim != 2 or data.shape[1] < 2: + local_print_log(f"Error. File {tdrw_filename} does not contain at least two columns.", "error") + + finite_md_times = data[:, 0] + finite_md_cdf = data[:, 1] + + return finite_md_times, finite_md_cdf + + +def limited_matrix_diffusion_from_file(self, G): + """ Matrix diffusion with limited block size + + Parameters + ---------- + G : NetworkX graph + graph obtained from graph_flow + + Returns + ------- + None + + Notes + ----------- + All parameters are attached to the particle class + """ + + # print(self.total_time,self.advect_time,self.matrix_diffusion_time) + # print("\nsegment limited sampling") + eps = 1e-4 + # print(self.advect_time, self.delta_t ) + # b = (2*self.delta_t) / self.beta + b = G.edges[self.curr_node, self.next_node]['b'] + # print(f"b: {b_eff}") + # # traping rate + gamma = (2*self.matrix_porosity*self.matrix_diffusivity)/(b * eps * self.fracture_spacing) + # gamma = (2*self.matrix_porosity*self.matrix_diffusivity)/(b * eps) + #print(f"gamma: {gamma}") + # average number of trapping events in gamma * advective time + average_number_of_trapping_events = self.delta_t * gamma + #print(f"average_number_of_trapping_events: {average_number_of_trapping_events}") + # number of trapping events in sampled from a poisson distribution + n = np.random.poisson(average_number_of_trapping_events) + # print(f"n: {n}") + xi = np.random.uniform(size = n) + # print(xi) + # + tmp = self.tau_D * np.interp(xi, self.trans_prob, self.transfer_time) + # print("*") + # print(self.transfer_time) + # print(self.trans_prob) + # print(tmp) + # exit(1) + # tmp = self.tau_D * self.inv_spline(xi) + self.delta_t_md = tmp.sum() + #print(self.delta_t_md) + #print("\n") + # self.total_time = self.advect_time + self.matrix_diffusion_time + #print(self.total_time,self.advect_time,self.matrix_diffusion_time) + + diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_roubinet.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_roubinet.py new file mode 100644 index 00000000..74703999 --- /dev/null +++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_roubinet.py @@ -0,0 +1,441 @@ + +import numpy as np +from scipy import special +import mpmath as mp + +from pydfnworks.general.logging import local_print_log +from pydfnworks.dfnGraph.tdrw_infinite import t_diff_unlimited + + +def get_fracture_segments(transfer_time, + fracture_length, + b, + velocity, + matrix_diffusivity, + matrix_porosity, + plim=0.01): + """ + This function calculates the maximum segment length that can be used in the TDRW model, following the protocol proposed by Roubinet et al., 2010. + + + Parameters + --------------------- + transfer_time : float + Time is takes for a particle to diffusive through the matrix across the fracture spacing. transfer_time = fracture_spacing**2 / (2 * matrix_diffusivity) + + fracture_length : float + Length of the current edge segment in the graph [m] + + b : float + Aperture of the current edge segment in the graph [m] + + velocity : float + particle velocity along the edge segment in the graph [m/s] + + matrix_diffusivity: float + Matrix Diffusivity [m^2/s] + + matrix_porosity: float + Matrix Porosity + + plim : float + Parameter used to break up the fracture + + Returns + --------------------- + segment_length : float + Length of the segment of the edge to compute the matrix diffusion time + + num_segments : int + Number of segments of length segment_length on the original edge/fracture. We ake the ceiling of the value. + + Notes + --------------------- + Roubinet use a plim value of 0.1. However, we are better able to match the results of Sudicky & Frind with a parameter values when plim = 0.01. + + """ + + # Compute the maximum segment length + segment_length = ((b * np.sqrt(transfer_time)) / + (matrix_porosity * np.sqrt(matrix_diffusivity)) + ) * special.erfcinv(1.0 - plim) * velocity + # if that's bigger than the edge size, just return the edge length as a single segment + if segment_length >= fracture_length: + return fracture_length, 1 + else: + # otherwise, get the number of segments of length segment_length on the edge + num_segments = int(np.ceil(fracture_length / segment_length)) + return segment_length, num_segments + +def transition_probability_cdf(t_min, t_max, frac_spacing, matrix_diffusivity, + num_pts): + """ + This function calculates the cummulative probability density that a particle changes fractures, given the time needed to reach a penetration depth and an associated time with each probabilty. + + Parameters + --------------- + t_min : float + Minumum value of diffusion time [s] + + t_max : float + Maximum value of diffusion time [s] + + frac_spacing : float + Spacing between fractures [m] + + matrix_diffusivity : float + Matrix Diffusivity value [m^2/s] + + num_pts : int + Number of points in the logspace array between t_min and t_max + + Returns + -------------- + times : np.array + Array of diffusion times + + prob_cdf : np.array + Array of cummulative probabilities. They only go to 0.5 + + + Notes + ------------ + Negative probabilities can sometimes appear due to numerical instabilities in our laplace inverse transform. These are removed in the section below marked CLEAN UP. We also ensure that our final distribution is monotonically increasing + + """ + local_print_log("--> Building transition probability cdf") + + times = np.logspace(np.log10(t_min), np.log10(t_max), num=num_pts) + prob_cdf = np.zeros(num_pts) + + ## Parameters for Roubinet function for the transfer probability function described + ## in Roubinet et al. WRR 2010. Equation number 5. + l1 = frac_spacing / 2 + l2 = -l1 + delta_l = l1 - l2 + + ## Roubinet et al. WRR 2010. Equation number 5. + ## s is the Laplace variable + laplace_trans_prob_function = lambda s: (mp.exp(l1*mp.sqrt(s/matrix_diffusivity))/s) \ + *((1.0 - mp.exp(-2*l2*mp.sqrt(s/matrix_diffusivity)))/(1.0 - mp.exp(2.0 *(delta_l)*mp.sqrt(s/matrix_diffusivity)))) + + for i, t in enumerate(times): + prob_cdf[i] = mp.invertlaplace(laplace_trans_prob_function, + t, + method='talbot', + dps=16, + degree=36) + + if prob_cdf[i] >= 0.5: + break + + # # CLEAN UP the solution. + # # Sometime small negative values will appear due to numerical instabilities. Remove these. + ind = np.asarray(np.where(prob_cdf >= 0)).flatten() + prob_cdf = prob_cdf[ind] + times = times[ind] + + # # Monotonically increasing + prob_cdf = np.maximum.accumulate(prob_cdf) + # unique values + prob_cdf, ind = np.unique(prob_cdf, return_index=True) + times = times[ind] + + return times, prob_cdf + + +def transfer_probabilities(b_min, + b_max, + tf_min, + tf_max, + matrix_porosity, + matrix_diffusivity, + frac_spacing, + eps=1e-16, + num_pts=100): + """ Returns the CDF of transfer probabilities and assocaited times. Lower and upper bounds are first estimated using the physical parameters of the system. The bounds are then tighted based on the returned probabilities wherein we find a range with probabilities greater than eps and 0.5 + + Parameters + ------------ + b_min : float + Minimum aperture in the network + + b_max : float + Maximum aperture in the network + + tf_min : float + Minimum advective travel time in the network + + tf_max : float + Maximum advective travel time in the network + + matrix_porosity: float + Matrix Porosity + + matrix_diffusivity : float + Matrix Diffusivity value [m^2/s] + + frac_spacing : float + Spacing between fractures [m] + + eps : float + Default - 1e-16 + + num_pts : int + Number of points in the logspace array between t_min and t_max + + Returns + ------------ + trans_prob : dictionary + Dictionary elements + times : np.array + Array of diffusion times + prob_cdf : np.array + Array of cummulative probabilities. They only go to 0.5 + + Notes + ------------ + The initial bounds are huge, don't worry. They get tighted up just fine. + + """ + # estimate lower bound + # get the minimum factor + a_min = (matrix_porosity * np.sqrt(matrix_diffusivity)) / b_max + # sample at the lowest value of tf, with sample ~ 0 + t_diff_lb = t_diff_unlimited(a_min, tf_min, eps) + # Below this value the inversse laplace transform does not convergence + if t_diff_lb < 1e-12: + local_print_log(f"lb too low {t_diff_lb:0.2e} changing to 1e-12", 'warning') + t_diff_lb = 1e-12 + + # estimate upper bound + # get the maximum factor + a_max = (matrix_porosity * np.sqrt(matrix_diffusivity)) / b_min + # sample at the largest value of tf, with sample ~ 1 + t_diff_ub = t_diff_unlimited(a_max, tf_max, 1 - eps) + # if t_diff_ub > 1e30: + # print(f"ub too high {t_diff_ub} changing to 1e30") + # t_diff_ub = 1e30 + + local_print_log( + f"--> Initial bounds for diffusion times. Min: {t_diff_lb:0.2e}, Max: {t_diff_ub:0.2e}" + ) + + # Compute the transition probabilities across the whole range of times. + times, trans_cdf = transition_probability_cdf(t_diff_lb, t_diff_ub, + frac_spacing, + matrix_diffusivity, num_pts) + + # Restricts the + # We end up with a lot of values we don't need, close to 0 probability. + # walk through the CDF and find the time corresponding to the + # first probability greater than epsilon + for i, val in enumerate(trans_cdf): + if val > eps: + t_diff_lb = times[i] + break + + # Likewise we end up with a lot of values we don't need, close to 0.5 probability. + # walk through the CDF and find the time corresponding to the + # first probability greater than 0.5 - epsilon + + for i, val in reversed(list(enumerate(trans_cdf))): + if val <= 0.5: + t_diff_ub = times[i] + break + + local_print_log( + f"--> Final bounds for diffusion times. Min: {t_diff_lb:0.2e}, Max: {t_diff_ub:0.2e}" + ) + + # Recompute the transition probabilities across the restricted range of times. + times, trans_cdf = transition_probability_cdf(t_diff_lb, t_diff_ub, + frac_spacing, + matrix_diffusivity, num_pts) + + #convert to as dictionary + trans_prob = {"times": times, "cdf": trans_cdf} + # print(times) + # print(trans_cdf) + return trans_prob + + +def segment_matrix_diffusion(trans_prob, matrix_porosity, matrix_diffusivity, + b, velocity, segment_length, num_segments): + """ Computes the time delay for a particle due to matrix diffusion on a given edge in the graph. The edge might already be broken into multiple segments depending on the parameters of the simulation. + + Parameters + -------------- + trans_prob : dictionary + Dictionary elements + times : np.array + Array of diffusion times + prob_cdf : np.array + Array of cummulative probabilities. They only go to 0.5 + + matrix_porosity: float + Matrix Porosity + + matrix_diffusivity : float + Matrix Diffusivity value [m^2/s] + + b : float + Aperture of the current edge segment in the graph [m] + + velocity : float + particle velocity along the edge segment in the graph [m/s] + + segment_length : float + Length of the segment of the edge to compute the matrix diffusion time + + num_segments : int + Number of segments of length segment_length on the original edge/fracture. We ake the ceiling of the value. + + + Returns + -------------- + t_diff : float + Total time delay due to retention via matrix diffusion + + Notes + ------------- + None + + """ + + a = (matrix_porosity * np.sqrt(matrix_diffusivity)) / b + + segment_time = segment_length / velocity + prob_min = min(trans_prob['cdf']) + prob_max = max(trans_prob['cdf']) + t_diff = 0 + for i in range(num_segments): + # sample a diffusion time from the unlimited scenario + xi = np.random.uniform(low=0, high=1) + t_diff_seg = t_diff_unlimited(a, segment_time, xi) + + # compare that time with limited times + # if it's less than the minimum time, there is 0 probability for + # transfer and we accept that time + if t_diff_seg < trans_prob['times'][0]: + limited_probability = 0 + + # if it's greater than the maximum time, there is probability 1 for + # transfer and we resample + elif t_diff_seg > trans_prob['times'][-1]: + limited_probability = 1 + + # otherwise, we get a new probability using the CDF based on the sampled time. + else: + limited_probability = 2 * np.interp( + t_diff_seg / 2, trans_prob['times'], trans_prob['cdf']) + + if limited_probability > 0: + # Draw a random number from U[0,1] for each particle. + # Particles transfer to a new fracture if this random number + # is less than the transfer probability. + # this is the inverse CDF method part of the algorithm. + xi = np.random.uniform(low=0, high=1) + if xi < limited_probability: + xi = np.random.uniform(low=prob_min, high=prob_max) + t_diff_seg = np.interp(xi, trans_prob['cdf'], + trans_prob['times']) + + + t_diff += t_diff_seg + return t_diff + + +def get_aperture_and_time_limits(G): + """ Walks through edges on the graph and returns min and max of aperture and advective travel times + + Parameters + --------------------- + G : networkX graph + Graph provided by graph_flow modules + + Returns + ------------- + b_min : float + Minimum value of aperture on the network + + b_max : float + Maximmum value of aperture on the network + + t_min : float + Minimum value of advective travel time on the network + + t_max : float + Maximmum value of advective travel time on the network + + Notes + -------------------- + Could be improved with nx.get_edge_attributes, maybe. JDH + + + """ + # get b_min, b_max, t_min, t_max + local_print_log("--> Getting b and t limits") + b_min = None + b_max = None + t_min = None + t_max = None + for u, v, d in G.edges(data=True): + if b_min is None: + b_min = d['b'] + elif d['b'] < b_min: + b_min = d['b'] + + if b_max is None: + b_max = d['b'] + elif d['b'] > b_max: + b_max = d['b'] + + if t_min is None: + t_min = d['time'] + elif d['time'] < t_min: + t_min = d['time'] + + if t_max is None: + t_max = d['time'] + elif d['time'] > t_max: + t_max = d['time'] + + local_print_log(f"--> b-min: {b_min:0.2e}, b-max: {b_max:0.2e}") + local_print_log(f"--> t-min: {t_min:0.2e}, t-max: {t_max:0.2e}") + return b_min, b_max, t_min, t_max + + +def limited_matrix_diffusion_roubinet(self, G): + """ Matrix diffusion with limited block size + + Parameters + ---------- + G : NetworkX graph + graph obtained from graph_flow + + Returns + ------- + None + + Notes + ----------- + All parameters are attached to the particle class + """ + + frac_length = G.edges[self.curr_node, self.next_node]['length'] + b = G.edges[self.curr_node, self.next_node]['b'] + velocity = G.edges[self.curr_node, self.next_node]['velocity'] + + segment_length, num_segments = get_fracture_segments( + self.transfer_time, frac_length, b, velocity, self.matrix_diffusivity, + self.matrix_porosity) + + self.delta_t_md = segment_matrix_diffusion(self.trans_prob, + self.matrix_porosity, + self.matrix_diffusivity, b, + velocity, segment_length, + num_segments) + + + diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_infinite.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_infinite.py new file mode 100644 index 00000000..6962e6fc --- /dev/null +++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_infinite.py @@ -0,0 +1,54 @@ + +from scipy import special +import numpy as np + +def t_diff_unlimited(a, tf, xi): + """ + This function returns one diffusion time sample for a and tf assuming an unlimited matrix block size. This is used throughout the limited block size sampling technique. + + Parameters + ------------------- + a : double + Constant parameter describing retention in the matrix. a = (matrix_porosity*matrix_diffusion)/aperture_length + + tf : double + Advective travel time + + xi : float + value between [0,1) + + Returns + -------------- + t_diff : float + Total time diffusing the matrix + + Notescx + ------------- + For a random sample, sample xi from U[0,1). + + """ + + return ((a * tf) / special.erfcinv(xi))**2 + +def unlimited_matrix_diffusion(self, G): + """ Matrix diffusion with unlimited block size + + Parameters + ---------- + G : NetworkX graph + graph obtained from graph_flow + + Returns + ---------- + None + + Notes + ----------- + All parameters are attached to the particle class + + """ + + b = G.edges[self.curr_node, self.next_node]['b'] + a_nondim = self.matrix_porosity * np.sqrt(self.matrix_diffusivity) / b + xi = np.random.uniform(size=1, low=0, high=1)[0] + self.delta_t_md = ((a_nondim * self.delta_t / special.erfcinv(xi))**2)