From a59c34b1cdc697b57debacc13511497ab69ca165 Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Mon, 24 Nov 2025 13:07:51 -0500
Subject: [PATCH 01/19] refactor of current methods into multiple files

---
 pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py  | 548 +++---------------
 .../pydfnworks/dfnGraph/graph_transport.py    | 254 ++------
 .../pydfnworks/dfnGraph/particle_class.py     |  28 +-
 3 files changed, 133 insertions(+), 697 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
index 31ccba13..3ea61449 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
@@ -3,438 +3,75 @@
 import mpmath as mp
 
 from pydfnworks.general.logging import local_print_log
+import pydfnworks.dfnGraph.tdrw_finite_roubinet as roubinet  
 
-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]
+def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model):
+    """ Check that the provided tdrw values are physiscal
 
-        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
-
-
-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
+    ----------
+        G: NetworkX graph 
+            Directed Graph obtained from output of graph_flow
 
     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
+    -------
+        dict : nested dictionary.
 
     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
+    -----
+        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
     
-    Notes
-    -------------
-        None
-
     """
+    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}")
+
+    # --- 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 ("roubinet", "dentz"):
+        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'])
-
-
-        t_diff += t_diff_seg
-    return t_diff
-
+            local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e}")
+    else:
+        local_print_log(
+            f"Error: Unknown TDRW model '{tdrw_model}'. Valid options are 'infinite', 'roubinet', or 'dentz' (case-sensitive).",
+            "error"
+        )
 
-def get_aperture_and_time_limits(G):
-    """ Walks through edges on the graph and returns min and max of aperture and advective travel times 
+    local_print_log("")  # clean trailing newline   
     
-    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 set_up_limited_matrix_diffusion(G,
-                                    frac_spacing,
+def _set_up_limited_matrix_diffusion(G,tdrw_model,
+                                    fracture_spacing,
                                     matrix_porosity,
                                     matrix_diffusivity,
                                     eps=1e-16,
@@ -476,65 +113,20 @@ def set_up_limited_matrix_diffusion(G,
         None
 
     """
+    if tdrw_model == 'roubinet':
+        # local_print_log("Using Roubinet Model for finite Matrix Diffusion")
+        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)
+    elif tdrw_model == "dentz":
+        local_print_log("Using Dentz Model for finite Matrix Diffusion") 
+        transfer_time = None
+        trans_prob = 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
-
-    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)
-
-
-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 
-
-    """
+    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..a623a27b 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -13,10 +13,12 @@
 
 # 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):
     """ Tracks a single particle through the graph
@@ -47,12 +49,18 @@ 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"])
+    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"])
 
     # # get current process information
     global nbrs_dict
@@ -65,181 +73,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,6 +82,7 @@ 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,
@@ -320,57 +154,63 @@ 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)
+        
+        if tdrw_model in ("roubinet", "dentz"):
+            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,
+                                                        tdrw_model,
+                                                        fracture_spacing,
+                                                        matrix_porosity,
+                                                        matrix_diffusivity) 
     else:
         trans_prob = None
         transfer_time = None
+
     ## main loop
+    print("Here")
     if self.ncpu == 1:
         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)
 
diff --git a/pydfnworks/pydfnworks/dfnGraph/particle_class.py b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
index 43bdbf98..af3d5bab 100644
--- a/pydfnworks/pydfnworks/dfnGraph/particle_class.py
+++ b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
@@ -15,11 +15,15 @@ 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
+    from pydfnworks.dfnGraph.tdrw_infinite import unlimited_matrix_diffusion
+    from pydfnworks.dfnGraph.tdrw_finite_dentz import unlimited_matrix_diffusion_dentz
+    from pydfnworks.dfnGraph.tdrw_finite_roubinet import limited_matrix_diffusion_roubinet 
 
-    def __init__(self, particle_number, ip, tdrw_flag, matrix_porosity,
-                 matrix_diffusivity, fracture_spacing, trans_prob,
-                 transfer_time, cp_flag, control_planes, direction):
+
+    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
@@ -39,6 +43,7 @@ def __init__(self, particle_number, ip, tdrw_flag, matrix_porosity,
         self.matrix_porosity = matrix_porosity
         self.matrix_diffusivity = matrix_diffusivity
         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,9 +56,6 @@ 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 = []
@@ -273,11 +275,13 @@ def track(self, G, nbrs_dict):
                 self.cleanup_frac_seq()
                 break
 
-            if self.tdrw_flag:
-                if self.fracture_spacing is None:
-                    self.unlimited_matrix_diffusion(G)
-                else:
-                    self.limited_matrix_diffusion(G)
+            if self.tdrw_flag: 
+                    if self.tdrw_model == "roubinet":
+                        self.limited_matrix_diffusion_roubinet(G)
+                    elif self.tdrw_model == "dentz":
+                        self.limited_matrix_diffusion_roubinet(G)
+                    elif self.tdrw_model == "infinite":
+                        self.unlimited_matrix_diffusion(G)
 
             if self.cp_flag:
                 self.cross_control_plane(G)

From ddab7cffaf62fd4c53a2c9c444a95307b6c8f7ff Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Mon, 24 Nov 2025 13:08:13 -0500
Subject: [PATCH 02/19] refactor fioes

---
 .../graph_transport_setup_functions.py        | 183 ++++++++
 .../pydfnworks/dfnGraph/tdrw_finite_dentz.py  |   0
 .../dfnGraph/tdrw_finite_roubinet.py          | 441 ++++++++++++++++++
 .../pydfnworks/dfnGraph/tdrw_infinite.py      |  54 +++
 4 files changed, 678 insertions(+)
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/graph_transport_setup_functions.py
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_roubinet.py
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_infinite.py

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/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
new file mode 100644
index 00000000..e69de29b
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)

From 257a3847111d9516f505fa7f59a3d8dc9a087a3d Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Mon, 24 Nov 2025 14:29:12 -0500
Subject: [PATCH 03/19] fixed parallel issues. Seems to be working for
 roubinet. Dentz is off

---
 pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py  |  13 +--
 .../pydfnworks/dfnGraph/graph_transport.py    |  41 ++++---
 .../pydfnworks/dfnGraph/particle_class.py     |  24 ++--
 .../pydfnworks/dfnGraph/tdrw_finite_dentz.py  | 109 ++++++++++++++++++
 4 files changed, 153 insertions(+), 34 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
index 3ea61449..bcdd0010 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
@@ -4,6 +4,7 @@
 
 from pydfnworks.general.logging import local_print_log
 import pydfnworks.dfnGraph.tdrw_finite_roubinet as roubinet  
+import pydfnworks.dfnGraph.tdrw_finite_dentz as dentz
 
 def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model):
     """ Check that the provided tdrw values are physiscal
@@ -35,7 +36,7 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
             "error"
         )
     else:
-        local_print_log(f"--> Matrix porosity value: {matrix_porosity:.2f}")
+        local_print_log(f"--> Matrix porosity value: {matrix_porosity:.2f} [-]")
 
     # --- Matrix Diffusivity ---
     if matrix_diffusivity is None:
@@ -43,7 +44,7 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
     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}")
+        local_print_log(f"--> Matrix diffusivity value: {matrix_diffusivity:.2e} [m^2/s]")
 
     # --- TDRW Model Selection ---
     if tdrw_model is None:
@@ -61,7 +62,7 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
         elif fracture_spacing <= 0:
             local_print_log(f"Error: Non-positive fracture_spacing={fracture_spacing}. Exiting program.", "error")
         else:
-            local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e}")
+            local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e} [m]")
     else:
         local_print_log(
             f"Error: Unknown TDRW model '{tdrw_model}'. Valid options are 'infinite', 'roubinet', or 'dentz' (case-sensitive).",
@@ -121,10 +122,8 @@ def _set_up_limited_matrix_diffusion(G,tdrw_model,
                                             fracture_spacing, eps, num_pts)
         transfer_time = fracture_spacing**2 / (2 * matrix_diffusivity)
     elif tdrw_model == "dentz":
-        local_print_log("Using Dentz Model for finite Matrix Diffusion") 
-        transfer_time = None
-        trans_prob = None 
-
+        # local_print_log("Using Dentz Model for finite Matrix Diffusion") 
+        transfer_time, trans_prob = dentz._make_inverse_cdf_spline_for_times()
     else:
         trans_prob = None
         transfer_time = None
diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
index a623a27b..06470ba9 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -49,18 +49,21 @@ def track_particle(data, verbose=False):
             f"--> Particle {data['particle_number']} is starting on worker {cpu_id}"
         )
 
-    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"])
+    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
@@ -88,7 +91,8 @@ def run_graph_transport(self,
                         fracture_spacing=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
@@ -185,14 +189,16 @@ def run_graph_transport(self,
                                                         tdrw_model,
                                                         fracture_spacing,
                                                         matrix_porosity,
-                                                        matrix_diffusivity) 
+                                                        matrix_diffusivity)
+        else:
+            trans_prob = None
+            transfer_time = None           
     else:
         trans_prob = None
         transfer_time = None
 
-    ## main loop
-    print("Here")
     if self.ncpu == 1:
+        self.print_log("--> Running in Serial")
         tic = timeit.default_timer()
         particles = []
         for i in range(nparticles):
@@ -245,6 +251,7 @@ def gather_output(output):
             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
@@ -254,7 +261,7 @@ def gather_output(output):
             data["control_planes"] = control_planes
             data["direction"] = direction
             pool.apply_async(track_particle,
-                             args=(data, ),
+                             args=(data, verbose),
                              callback=gather_output)
 
         pool.close()
diff --git a/pydfnworks/pydfnworks/dfnGraph/particle_class.py b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
index af3d5bab..7ac646c8 100644
--- a/pydfnworks/pydfnworks/dfnGraph/particle_class.py
+++ b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
@@ -16,11 +16,11 @@ class Particle():
     '''
 
     from pydfnworks.dfnGraph.tdrw_infinite import unlimited_matrix_diffusion
-    from pydfnworks.dfnGraph.tdrw_finite_dentz import unlimited_matrix_diffusion_dentz
+    from pydfnworks.dfnGraph.tdrw_finite_dentz import limited_matrix_diffusion_dentz
     from pydfnworks.dfnGraph.tdrw_finite_roubinet import limited_matrix_diffusion_roubinet 
 
-
-    def __init__(self, particle_number, ip, tdrw_flag, tdrw_model, matrix_porosity,
+    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):
         
@@ -42,6 +42,7 @@ def __init__(self, particle_number, ip, tdrw_flag, tdrw_model, 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
@@ -62,6 +63,9 @@ def __init__(self, particle_number, ip, tdrw_flag, tdrw_model, matrix_porosity,
         self.times = []
         self.coords = []
 
+        if tdrw_model == "dentz":
+            # self.tau_D = self.fracture_spacing**2/(self.matrix_diffusivity*self.matrix_porosity) 
+            self.tau_D = self.fracture_spacing**2/(self.matrix_diffusivity) 
 
     def initalize(self,G):
         self.frac = (G.nodes[self.curr_node]['frac'][1])
@@ -276,13 +280,13 @@ def track(self, G, nbrs_dict):
                 break
 
             if self.tdrw_flag: 
-                    if self.tdrw_model == "roubinet":
-                        self.limited_matrix_diffusion_roubinet(G)
-                    elif self.tdrw_model == "dentz":
-                        self.limited_matrix_diffusion_roubinet(G)
-                    elif self.tdrw_model == "infinite":
-                        self.unlimited_matrix_diffusion(G)
-
+                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 == "infinite":
+                    self.unlimited_matrix_diffusion(G)
+                
             if self.cp_flag:
                 self.cross_control_plane(G)
             self.update()
diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
index e69de29b..cd6c1787 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
@@ -0,0 +1,109 @@
+
+
+from scipy.interpolate import PchipInterpolator
+import numpy as np
+import matplotlib.pyplot as plt
+import mpmath as mp
+
+
+def Psi_star(lmbda, eps):
+    # Laplace transform of the CDF Ψ*_ε(λ)
+    # eps must satisfy 0 < eps < 1 for a slab of unit width in the note
+    sqrt_l = mp.sqrt(lmbda)
+    return mp.cosh((1 - eps) * sqrt_l) / (lmbda * mp.cosh(sqrt_l))
+
+def Psi_cdf(t, eps, method="dehoog"):
+    # Inverse Laplace transform of Ψ*_ε at time t
+    # method can be "dehoog" or "talbot"
+    if t <= 0:
+        return mp.mpf("0.0")
+    F = lambda s: Psi_star(s, eps)
+    return mp.invertlaplace(F, t, method=method)
+
+def Psi_pdf(t, eps, method="dehoog", h=None):
+    # Density ψ_ε(t) = d/dt Ψ_ε(t)
+    # You can either invert ψ*_ε directly or differentiate the CDF numerically
+    # This version inverts ψ*_ε for accuracy
+    if t <= 0:
+        return mp.mpf("0.0")
+    psi_star = lambda s: mp.cosh((1 - eps) * mp.sqrt(s)) / mp.cosh(mp.sqrt(s))
+    return mp.invertlaplace(psi_star, t, method=method)
+
+
+def _make_inverse_cdf_spline_for_times(num_samples = 50, eps = 1e-3, precision = 40):
+    # --- compute CDF and PDF values ---
+    mp.mp.dps = precision
+    num_samples = num_samples
+    eps = mp.mpf(eps)
+    times = [10**x for x in np.linspace(-8, 1, num_samples)]
+
+    cdf_vals = np.zeros(num_samples, dtype=float)
+
+    for i, t in enumerate(times):
+        cdf_vals[i] = float(Psi_cdf(t, eps, method="dehoog"))
+
+    # --- build inverse CDF spline (t vs cdf) ---
+    order = np.argsort(cdf_vals)
+    cdf_sorted = np.clip(cdf_vals[order], 0, 1)
+    t_sorted = np.maximum(np.array(times)[order], 0)
+
+    # remove duplicates
+    cdf_unique, idx = np.unique(cdf_sorted, return_index=True)
+    t_unique = np.asarray(t_sorted[idx])
+
+    finite_md_times = t_unique
+    finite_md_cdf = cdf_unique
+
+    return finite_md_times, finite_md_cdf 
+
+
+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 
+    """
+
+    # print(self.total_time,self.advect_time,self.matrix_diffusion_time)
+    # print("\nsegment limited sampling")
+    eps = 1e-3
+    # 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)
+
+

From d0d15db7d39b4ac6969102324cf455edf4280494 Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Wed, 26 Nov 2025 09:24:34 -0500
Subject: [PATCH 04/19] fixed epsilon and print statements

---
 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
index cd6c1787..489d318b 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
@@ -30,12 +30,12 @@ def Psi_pdf(t, eps, method="dehoog", h=None):
     return mp.invertlaplace(psi_star, t, method=method)
 
 
-def _make_inverse_cdf_spline_for_times(num_samples = 50, eps = 1e-3, precision = 40):
+def _make_inverse_cdf_spline_for_times(num_samples = 100, eps = 1e-4, precision = 40):
     # --- compute CDF and PDF values ---
     mp.mp.dps = precision
     num_samples = num_samples
     eps = mp.mpf(eps)
-    times = [10**x for x in np.linspace(-8, 1, num_samples)]
+    times = [10**x for x in np.linspace(-8, 3, num_samples)]
 
     cdf_vals = np.zeros(num_samples, dtype=float)
 
@@ -76,7 +76,7 @@ def limited_matrix_diffusion_dentz(self,G):
 
     # print(self.total_time,self.advect_time,self.matrix_diffusion_time)
     # print("\nsegment limited sampling")
-    eps = 1e-3
+    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']
@@ -93,11 +93,12 @@ def limited_matrix_diffusion_dentz(self,G):
     # 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)
+    # print(tmp)
     # exit(1) 
     # tmp = self.tau_D * self.inv_spline(xi)
     self.delta_t_md = tmp.sum() 

From 3d7566f75e2b444add3e1233be2119a5b561d889 Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Mon, 24 Nov 2025 13:07:51 -0500
Subject: [PATCH 05/19] refactor of current methods into multiple files

---
 pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py  | 548 +++---------------
 .../pydfnworks/dfnGraph/graph_transport.py    | 254 ++------
 .../pydfnworks/dfnGraph/particle_class.py     |  28 +-
 3 files changed, 133 insertions(+), 697 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
index 31ccba13..3ea61449 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
@@ -3,438 +3,75 @@
 import mpmath as mp
 
 from pydfnworks.general.logging import local_print_log
+import pydfnworks.dfnGraph.tdrw_finite_roubinet as roubinet  
 
-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]
+def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model):
+    """ Check that the provided tdrw values are physiscal
 
-        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
-
-
-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
+    ----------
+        G: NetworkX graph 
+            Directed Graph obtained from output of graph_flow
 
     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
+    -------
+        dict : nested dictionary.
 
     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
+    -----
+        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
     
-    Notes
-    -------------
-        None
-
     """
+    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}")
+
+    # --- 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 ("roubinet", "dentz"):
+        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'])
-
-
-        t_diff += t_diff_seg
-    return t_diff
-
+            local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e}")
+    else:
+        local_print_log(
+            f"Error: Unknown TDRW model '{tdrw_model}'. Valid options are 'infinite', 'roubinet', or 'dentz' (case-sensitive).",
+            "error"
+        )
 
-def get_aperture_and_time_limits(G):
-    """ Walks through edges on the graph and returns min and max of aperture and advective travel times 
+    local_print_log("")  # clean trailing newline   
     
-    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 set_up_limited_matrix_diffusion(G,
-                                    frac_spacing,
+def _set_up_limited_matrix_diffusion(G,tdrw_model,
+                                    fracture_spacing,
                                     matrix_porosity,
                                     matrix_diffusivity,
                                     eps=1e-16,
@@ -476,65 +113,20 @@ def set_up_limited_matrix_diffusion(G,
         None
 
     """
+    if tdrw_model == 'roubinet':
+        # local_print_log("Using Roubinet Model for finite Matrix Diffusion")
+        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)
+    elif tdrw_model == "dentz":
+        local_print_log("Using Dentz Model for finite Matrix Diffusion") 
+        transfer_time = None
+        trans_prob = 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
-
-    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)
-
-
-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 
-
-    """
+    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..a623a27b 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -13,10 +13,12 @@
 
 # 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):
     """ Tracks a single particle through the graph
@@ -47,12 +49,18 @@ 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"])
+    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"])
 
     # # get current process information
     global nbrs_dict
@@ -65,181 +73,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,6 +82,7 @@ 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,
@@ -320,57 +154,63 @@ 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)
+        
+        if tdrw_model in ("roubinet", "dentz"):
+            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,
+                                                        tdrw_model,
+                                                        fracture_spacing,
+                                                        matrix_porosity,
+                                                        matrix_diffusivity) 
     else:
         trans_prob = None
         transfer_time = None
+
     ## main loop
+    print("Here")
     if self.ncpu == 1:
         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)
 
diff --git a/pydfnworks/pydfnworks/dfnGraph/particle_class.py b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
index 43bdbf98..af3d5bab 100644
--- a/pydfnworks/pydfnworks/dfnGraph/particle_class.py
+++ b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
@@ -15,11 +15,15 @@ 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
+    from pydfnworks.dfnGraph.tdrw_infinite import unlimited_matrix_diffusion
+    from pydfnworks.dfnGraph.tdrw_finite_dentz import unlimited_matrix_diffusion_dentz
+    from pydfnworks.dfnGraph.tdrw_finite_roubinet import limited_matrix_diffusion_roubinet 
 
-    def __init__(self, particle_number, ip, tdrw_flag, matrix_porosity,
-                 matrix_diffusivity, fracture_spacing, trans_prob,
-                 transfer_time, cp_flag, control_planes, direction):
+
+    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
@@ -39,6 +43,7 @@ def __init__(self, particle_number, ip, tdrw_flag, matrix_porosity,
         self.matrix_porosity = matrix_porosity
         self.matrix_diffusivity = matrix_diffusivity
         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,9 +56,6 @@ 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 = []
@@ -273,11 +275,13 @@ def track(self, G, nbrs_dict):
                 self.cleanup_frac_seq()
                 break
 
-            if self.tdrw_flag:
-                if self.fracture_spacing is None:
-                    self.unlimited_matrix_diffusion(G)
-                else:
-                    self.limited_matrix_diffusion(G)
+            if self.tdrw_flag: 
+                    if self.tdrw_model == "roubinet":
+                        self.limited_matrix_diffusion_roubinet(G)
+                    elif self.tdrw_model == "dentz":
+                        self.limited_matrix_diffusion_roubinet(G)
+                    elif self.tdrw_model == "infinite":
+                        self.unlimited_matrix_diffusion(G)
 
             if self.cp_flag:
                 self.cross_control_plane(G)

From 485ac596b06447db8dd785bfc29e60c27404ee45 Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Mon, 24 Nov 2025 13:08:13 -0500
Subject: [PATCH 06/19] refactor fioes

---
 .../graph_transport_setup_functions.py        | 183 ++++++++
 .../pydfnworks/dfnGraph/tdrw_finite_dentz.py  |   0
 .../dfnGraph/tdrw_finite_roubinet.py          | 441 ++++++++++++++++++
 .../pydfnworks/dfnGraph/tdrw_infinite.py      |  54 +++
 4 files changed, 678 insertions(+)
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/graph_transport_setup_functions.py
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_roubinet.py
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_infinite.py

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/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
new file mode 100644
index 00000000..e69de29b
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)

From 70a2fb7ae79b42d391831ff0fcbdfafefe8edff7 Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Mon, 24 Nov 2025 14:29:12 -0500
Subject: [PATCH 07/19] fixed parallel issues. Seems to be working for
 roubinet. Dentz is off

---
 pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py  |  13 +--
 .../pydfnworks/dfnGraph/graph_transport.py    |  41 ++++---
 .../pydfnworks/dfnGraph/particle_class.py     |  24 ++--
 .../pydfnworks/dfnGraph/tdrw_finite_dentz.py  | 109 ++++++++++++++++++
 4 files changed, 153 insertions(+), 34 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
index 3ea61449..bcdd0010 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
@@ -4,6 +4,7 @@
 
 from pydfnworks.general.logging import local_print_log
 import pydfnworks.dfnGraph.tdrw_finite_roubinet as roubinet  
+import pydfnworks.dfnGraph.tdrw_finite_dentz as dentz
 
 def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model):
     """ Check that the provided tdrw values are physiscal
@@ -35,7 +36,7 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
             "error"
         )
     else:
-        local_print_log(f"--> Matrix porosity value: {matrix_porosity:.2f}")
+        local_print_log(f"--> Matrix porosity value: {matrix_porosity:.2f} [-]")
 
     # --- Matrix Diffusivity ---
     if matrix_diffusivity is None:
@@ -43,7 +44,7 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
     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}")
+        local_print_log(f"--> Matrix diffusivity value: {matrix_diffusivity:.2e} [m^2/s]")
 
     # --- TDRW Model Selection ---
     if tdrw_model is None:
@@ -61,7 +62,7 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
         elif fracture_spacing <= 0:
             local_print_log(f"Error: Non-positive fracture_spacing={fracture_spacing}. Exiting program.", "error")
         else:
-            local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e}")
+            local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e} [m]")
     else:
         local_print_log(
             f"Error: Unknown TDRW model '{tdrw_model}'. Valid options are 'infinite', 'roubinet', or 'dentz' (case-sensitive).",
@@ -121,10 +122,8 @@ def _set_up_limited_matrix_diffusion(G,tdrw_model,
                                             fracture_spacing, eps, num_pts)
         transfer_time = fracture_spacing**2 / (2 * matrix_diffusivity)
     elif tdrw_model == "dentz":
-        local_print_log("Using Dentz Model for finite Matrix Diffusion") 
-        transfer_time = None
-        trans_prob = None 
-
+        # local_print_log("Using Dentz Model for finite Matrix Diffusion") 
+        transfer_time, trans_prob = dentz._make_inverse_cdf_spline_for_times()
     else:
         trans_prob = None
         transfer_time = None
diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
index a623a27b..06470ba9 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -49,18 +49,21 @@ def track_particle(data, verbose=False):
             f"--> Particle {data['particle_number']} is starting on worker {cpu_id}"
         )
 
-    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"])
+    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
@@ -88,7 +91,8 @@ def run_graph_transport(self,
                         fracture_spacing=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
@@ -185,14 +189,16 @@ def run_graph_transport(self,
                                                         tdrw_model,
                                                         fracture_spacing,
                                                         matrix_porosity,
-                                                        matrix_diffusivity) 
+                                                        matrix_diffusivity)
+        else:
+            trans_prob = None
+            transfer_time = None           
     else:
         trans_prob = None
         transfer_time = None
 
-    ## main loop
-    print("Here")
     if self.ncpu == 1:
+        self.print_log("--> Running in Serial")
         tic = timeit.default_timer()
         particles = []
         for i in range(nparticles):
@@ -245,6 +251,7 @@ def gather_output(output):
             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
@@ -254,7 +261,7 @@ def gather_output(output):
             data["control_planes"] = control_planes
             data["direction"] = direction
             pool.apply_async(track_particle,
-                             args=(data, ),
+                             args=(data, verbose),
                              callback=gather_output)
 
         pool.close()
diff --git a/pydfnworks/pydfnworks/dfnGraph/particle_class.py b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
index af3d5bab..7ac646c8 100644
--- a/pydfnworks/pydfnworks/dfnGraph/particle_class.py
+++ b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
@@ -16,11 +16,11 @@ class Particle():
     '''
 
     from pydfnworks.dfnGraph.tdrw_infinite import unlimited_matrix_diffusion
-    from pydfnworks.dfnGraph.tdrw_finite_dentz import unlimited_matrix_diffusion_dentz
+    from pydfnworks.dfnGraph.tdrw_finite_dentz import limited_matrix_diffusion_dentz
     from pydfnworks.dfnGraph.tdrw_finite_roubinet import limited_matrix_diffusion_roubinet 
 
-
-    def __init__(self, particle_number, ip, tdrw_flag, tdrw_model, matrix_porosity,
+    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):
         
@@ -42,6 +42,7 @@ def __init__(self, particle_number, ip, tdrw_flag, tdrw_model, 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
@@ -62,6 +63,9 @@ def __init__(self, particle_number, ip, tdrw_flag, tdrw_model, matrix_porosity,
         self.times = []
         self.coords = []
 
+        if tdrw_model == "dentz":
+            # self.tau_D = self.fracture_spacing**2/(self.matrix_diffusivity*self.matrix_porosity) 
+            self.tau_D = self.fracture_spacing**2/(self.matrix_diffusivity) 
 
     def initalize(self,G):
         self.frac = (G.nodes[self.curr_node]['frac'][1])
@@ -276,13 +280,13 @@ def track(self, G, nbrs_dict):
                 break
 
             if self.tdrw_flag: 
-                    if self.tdrw_model == "roubinet":
-                        self.limited_matrix_diffusion_roubinet(G)
-                    elif self.tdrw_model == "dentz":
-                        self.limited_matrix_diffusion_roubinet(G)
-                    elif self.tdrw_model == "infinite":
-                        self.unlimited_matrix_diffusion(G)
-
+                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 == "infinite":
+                    self.unlimited_matrix_diffusion(G)
+                
             if self.cp_flag:
                 self.cross_control_plane(G)
             self.update()
diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
index e69de29b..cd6c1787 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
@@ -0,0 +1,109 @@
+
+
+from scipy.interpolate import PchipInterpolator
+import numpy as np
+import matplotlib.pyplot as plt
+import mpmath as mp
+
+
+def Psi_star(lmbda, eps):
+    # Laplace transform of the CDF Ψ*_ε(λ)
+    # eps must satisfy 0 < eps < 1 for a slab of unit width in the note
+    sqrt_l = mp.sqrt(lmbda)
+    return mp.cosh((1 - eps) * sqrt_l) / (lmbda * mp.cosh(sqrt_l))
+
+def Psi_cdf(t, eps, method="dehoog"):
+    # Inverse Laplace transform of Ψ*_ε at time t
+    # method can be "dehoog" or "talbot"
+    if t <= 0:
+        return mp.mpf("0.0")
+    F = lambda s: Psi_star(s, eps)
+    return mp.invertlaplace(F, t, method=method)
+
+def Psi_pdf(t, eps, method="dehoog", h=None):
+    # Density ψ_ε(t) = d/dt Ψ_ε(t)
+    # You can either invert ψ*_ε directly or differentiate the CDF numerically
+    # This version inverts ψ*_ε for accuracy
+    if t <= 0:
+        return mp.mpf("0.0")
+    psi_star = lambda s: mp.cosh((1 - eps) * mp.sqrt(s)) / mp.cosh(mp.sqrt(s))
+    return mp.invertlaplace(psi_star, t, method=method)
+
+
+def _make_inverse_cdf_spline_for_times(num_samples = 50, eps = 1e-3, precision = 40):
+    # --- compute CDF and PDF values ---
+    mp.mp.dps = precision
+    num_samples = num_samples
+    eps = mp.mpf(eps)
+    times = [10**x for x in np.linspace(-8, 1, num_samples)]
+
+    cdf_vals = np.zeros(num_samples, dtype=float)
+
+    for i, t in enumerate(times):
+        cdf_vals[i] = float(Psi_cdf(t, eps, method="dehoog"))
+
+    # --- build inverse CDF spline (t vs cdf) ---
+    order = np.argsort(cdf_vals)
+    cdf_sorted = np.clip(cdf_vals[order], 0, 1)
+    t_sorted = np.maximum(np.array(times)[order], 0)
+
+    # remove duplicates
+    cdf_unique, idx = np.unique(cdf_sorted, return_index=True)
+    t_unique = np.asarray(t_sorted[idx])
+
+    finite_md_times = t_unique
+    finite_md_cdf = cdf_unique
+
+    return finite_md_times, finite_md_cdf 
+
+
+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 
+    """
+
+    # print(self.total_time,self.advect_time,self.matrix_diffusion_time)
+    # print("\nsegment limited sampling")
+    eps = 1e-3
+    # 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)
+
+

From ad3b2dfa8c08962cacfa64012bf102d5b7ef123d Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Wed, 26 Nov 2025 09:24:34 -0500
Subject: [PATCH 08/19] fixed epsilon and print statements

---
 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
index cd6c1787..489d318b 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
@@ -30,12 +30,12 @@ def Psi_pdf(t, eps, method="dehoog", h=None):
     return mp.invertlaplace(psi_star, t, method=method)
 
 
-def _make_inverse_cdf_spline_for_times(num_samples = 50, eps = 1e-3, precision = 40):
+def _make_inverse_cdf_spline_for_times(num_samples = 100, eps = 1e-4, precision = 40):
     # --- compute CDF and PDF values ---
     mp.mp.dps = precision
     num_samples = num_samples
     eps = mp.mpf(eps)
-    times = [10**x for x in np.linspace(-8, 1, num_samples)]
+    times = [10**x for x in np.linspace(-8, 3, num_samples)]
 
     cdf_vals = np.zeros(num_samples, dtype=float)
 
@@ -76,7 +76,7 @@ def limited_matrix_diffusion_dentz(self,G):
 
     # print(self.total_time,self.advect_time,self.matrix_diffusion_time)
     # print("\nsegment limited sampling")
-    eps = 1e-3
+    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']
@@ -93,11 +93,12 @@ def limited_matrix_diffusion_dentz(self,G):
     # 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)
+    # print(tmp)
     # exit(1) 
     # tmp = self.tau_D * self.inv_spline(xi)
     self.delta_t_md = tmp.sum() 

From a6e86df56a35afa8a48b61601db04f63763a0ec5 Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Tue, 10 Mar 2026 15:10:02 -0400
Subject: [PATCH 09/19] matrix diffusion from file option

---
 pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py  | 21 ++++++++++++++++---
 .../pydfnworks/dfnGraph/graph_transport.py    | 17 ++++++++-------
 2 files changed, 28 insertions(+), 10 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
index bcdd0010..d5dbcf28 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
@@ -1,12 +1,15 @@
 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_from_file as from_file
 
-def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model):
+
+def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model, tdrw_filename):
     """ Check that the provided tdrw values are physiscal
 
 
@@ -63,6 +66,14 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
             local_print_log(f"Error: Non-positive fracture_spacing={fracture_spacing}. Exiting program.", "error")
         else:
             local_print_log(f"--> Fracture spacing value: {fracture_spacing:.2e} [m]")
+    elif tdrw_model == "from_file":
+        if tdrw_filename is None:
+            local_print_log(f"Error. TDRW Filename not provides", "error")
+ 
+        if not os.path.isfile(tdrw_filename):
+            local_print_log(f"Error. File not found: {tdrw_filename}", "error")
+
+        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', or 'dentz' (case-sensitive).",
@@ -71,12 +82,14 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
 
     local_print_log("")  # clean trailing newline   
     
-def _set_up_limited_matrix_diffusion(G,tdrw_model,
+def _set_up_limited_matrix_diffusion(G,
+                                    tdrw_model,
                                     fracture_spacing,
                                     matrix_porosity,
                                     matrix_diffusivity,
                                     eps=1e-16,
-                                    num_pts=100):
+                                    num_pts=100,
+                                    tdrw_filename = None):
     """ Sets up transition probabilities for limited block size matrix diffusion
     
     Parameters
@@ -124,6 +137,8 @@ def _set_up_limited_matrix_diffusion(G,tdrw_model,
     elif tdrw_model == "dentz":
         # local_print_log("Using Dentz Model for finite Matrix Diffusion") 
         transfer_time, trans_prob = dentz._make_inverse_cdf_spline_for_times()
+    elif tdrw_model == "from_file":
+        transfer_time, trans_prob = from_file._load_finite_time_cdf(tdrw_filename)
     else:
         trans_prob = None
         transfer_time = None
diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
index 06470ba9..b12e71a4 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -89,6 +89,7 @@ def run_graph_transport(self,
                         matrix_porosity=None,
                         matrix_diffusivity=None,
                         fracture_spacing=None,
+                        tdrw_filename = None, 
                         control_planes=None,
                         direction=None,
                         cp_filename='control_planes',
@@ -180,16 +181,18 @@ def run_graph_transport(self,
     # Check parameters for TDRW
     if tdrw_flag:
         _check_tdrw_params(matrix_porosity, matrix_diffusivity,
-                          fracture_spacing, tdrw_model)
+                          fracture_spacing, tdrw_model, tdrw_filename)
         
-        if tdrw_model in ("roubinet", "dentz"):
+        if tdrw_model in ("roubinet", "dentz", "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,
-                                                        tdrw_model,
-                                                        fracture_spacing,
-                                                        matrix_porosity,
-                                                        matrix_diffusivity)
+            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           

From 9e6e91231b6524eebb2e48de11866fc2e1330f44 Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Tue, 10 Mar 2026 15:10:15 -0400
Subject: [PATCH 10/19] file for MD from file

---
 .../dfnGraph/tdrw_finite_from_file.py         | 73 +++++++++++++++++++
 1 file changed, 73 insertions(+)
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_from_file.py

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

From 1cc951c7a989bf3340a5f3159cf32edcf7cc00ac Mon Sep 17 00:00:00 2001
From: Jeffrey Hyman <jhyman@lanl.gov>
Date: Tue, 10 Mar 2026 15:33:15 -0400
Subject: [PATCH 11/19] annulus runs

---
 pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py  | 135 +++++++++++-------
 .../pydfnworks/dfnGraph/graph_transport.py    |   2 +-
 .../pydfnworks/dfnGraph/particle_class.py     |  29 ++--
 3 files changed, 99 insertions(+), 67 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
index d5dbcf28..db8f786b 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_tdrw.py
@@ -6,27 +6,44 @@
 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 _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model, tdrw_filename):
-    """ Check that the provided tdrw values are physiscal
+# Valid finite matrix diffusion models that require fracture_spacing
+_FINITE_MODELS_REQUIRING_SPACING = ("roubinet", "dentz", "annulus")
+
 
+def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, tdrw_model, tdrw_filename):
+    """ Check that the provided tdrw values are physical
 
     Parameters
     ----------
-        G: NetworkX graph 
-            Directed Graph obtained from output of graph_flow
+        matrix_porosity : float
+            Porosity of the rock matrix [-]
+
+        matrix_diffusivity : float
+            Effective diffusivity of the rock matrix [m^2/s]
+
+        fracture_spacing : float
+            Characteristic matrix block half-width / fracture spacing [m]
+            Required for finite matrix models.
+
+        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
     -------
-        dict : nested dictionary.
+        None
 
     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
-    
+        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.")
 
@@ -52,9 +69,11 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
     # --- 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 ("roubinet", "dentz"):
+
+    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:
@@ -66,79 +85,97 @@ def _check_tdrw_params(matrix_porosity, matrix_diffusivity, fracture_spacing, td
             local_print_log(f"Error: Non-positive fracture_spacing={fracture_spacing}. Exiting program.", "error")
         else:
             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.")
+
     elif tdrw_model == "from_file":
         if tdrw_filename is None:
-            local_print_log(f"Error. TDRW Filename not provides", "error")
- 
-        if not os.path.isfile(tdrw_filename):
+            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}.")
 
-        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', or 'dentz' (case-sensitive).",
+            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):
+                                     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
-    
+
     Parameters
-    ---------------------
-        G : networkX graph  
+    ----------
+        G : networkX graph
             Graph provided by graph_flow modules
 
-        fracture_length : float
-            Length of the current edge segment in the graph [m]
+        tdrw_model : str
+            Name of the TDRW matrix diffusion model. Valid options:
+            'roubinet', 'dentz', 'annulus', 'from_file'.
+
+        fracture_spacing : float
+            Characteristic matrix block half-width [m].
+            Required for 'roubinet', 'dentz', and 'annulus' models.
 
-        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.
+
+        tdrw_filename : str or None
+            Path to CDF file, required when tdrw_model == 'from_file'.
 
     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
+    -------
+        transfer_time : np.ndarray or None
+            Array of dimensionless diffusion return times for inverse CDF lookup.
 
-    Notes
-    --------------------
-        None
+        trans_prob : np.ndarray or None
+            Array of cumulative probabilities corresponding to transfer_time.
 
+    Notes
+    -----
+        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':
-        # local_print_log("Using Roubinet Model for finite Matrix Diffusion")
         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)
+                                                     matrix_porosity, matrix_diffusivity,
+                                                     fracture_spacing, eps, num_pts)
         transfer_time = fracture_spacing**2 / (2 * matrix_diffusivity)
+
     elif tdrw_model == "dentz":
-        # local_print_log("Using Dentz Model for finite Matrix Diffusion") 
         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)
+
     elif tdrw_model == "from_file":
         transfer_time, trans_prob = from_file._load_finite_time_cdf(tdrw_filename)
+
     else:
         trans_prob = None
         transfer_time = None
diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
index b12e71a4..ab9e3187 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -183,7 +183,7 @@ def run_graph_transport(self,
         _check_tdrw_params(matrix_porosity, matrix_diffusivity,
                           fracture_spacing, tdrw_model, tdrw_filename)
         
-        if tdrw_model in ("roubinet", "dentz", "from_file"):
+        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,
diff --git a/pydfnworks/pydfnworks/dfnGraph/particle_class.py b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
index 7ac646c8..4aec9080 100644
--- a/pydfnworks/pydfnworks/dfnGraph/particle_class.py
+++ b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
@@ -17,7 +17,8 @@ class Particle():
 
     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_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,
@@ -64,8 +65,11 @@ def __init__(self, particle_number, ip, tdrw_flag, tdrw_model,
         self.coords = []
 
         if tdrw_model == "dentz":
-            # self.tau_D = self.fracture_spacing**2/(self.matrix_diffusivity*self.matrix_porosity) 
-            self.tau_D = self.fracture_spacing**2/(self.matrix_diffusivity) 
+            self.tau_D = self.fracture_spacing**2 / (self.matrix_diffusivity)
+
+        elif tdrw_model == "annulus":
+            # tauD = w^2 / D where w is the annular gap width (fracture_spacing)
+            self.tau_D = self.fracture_spacing**2 / (self.matrix_diffusivity)
 
     def initalize(self,G):
         self.frac = (G.nodes[self.curr_node]['frac'][1])
@@ -178,21 +182,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
@@ -203,7 +198,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
@@ -234,7 +229,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
@@ -275,18 +269,19 @@ 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.tdrw_flag:
                 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)
-                
+
             if self.cp_flag:
                 self.cross_control_plane(G)
             self.update()

From d2e8bd07a2822ed4bf431782871eb5ad00ca7475 Mon Sep 17 00:00:00 2001
From: Jeffrey Hyman <jhyman@lanl.gov>
Date: Tue, 10 Mar 2026 15:34:08 -0400
Subject: [PATCH 12/19] annulus file

---
 .../dfnGraph/tdrw_finite_annulus.py           | 108 ++++++++++++++++++
 1 file changed, 108 insertions(+)
 create mode 100644 pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py

diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
new file mode 100644
index 00000000..5ba20dbc
--- /dev/null
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
@@ -0,0 +1,108 @@
+
+import numpy as np
+import mpmath as mp
+
+
+def Psi_star_annulus(lmbda, eps):
+    # Laplace transform of the CDF for diffusion return time in an annulus
+    # Geometry: absorbing wall at one end, reflecting wall at the other
+    # eps = r'/w, dimensionless release position from the absorbing wall (0 < eps <= 1)
+    # F(s) = [ exp(-eps*q) + exp(-(2-eps)*q) ] / [ (1 + exp(-2q)) * s ]
+    # where q = sqrt(s * tauD);  tauD absorbed into s at call time (s -> s*tauD)
+    sqrt_l = mp.sqrt(lmbda)
+    numerator   = mp.exp(-eps * sqrt_l) + mp.exp(-(2 - eps) * sqrt_l)
+    denominator = (1 + mp.exp(-2 * sqrt_l)) * lmbda
+    return numerator / denominator
+
+
+def Psi_cdf_annulus(t, eps, method="dehoog"):
+    # Inverse Laplace transform of Psi_star_annulus at time t
+    # Returns cumulative return-time distribution P(T_return <= t)
+    if t <= 0:
+        return mp.mpf("0.0")
+    F = lambda s: Psi_star_annulus(s, eps)
+    return mp.invertlaplace(F, t, method=method)
+
+
+def Psi_pdf_annulus(t, eps, method="dehoog"):
+    # Density psi(t) = d/dt Psi(t), the first-passage time PDF
+    # Inverts the PDF Laplace transform directly for accuracy
+    # psi*(s) = [ exp(-eps*q) + exp(-(2-eps)*q) ] / (1 + exp(-2q))
+    if t <= 0:
+        return mp.mpf("0.0")
+    psi_star = lambda s: (mp.exp(-eps * mp.sqrt(s)) + mp.exp(-(2 - eps) * mp.sqrt(s))) / (1 + mp.exp(-2 * mp.sqrt(s)))
+    return mp.invertlaplace(psi_star, t, method=method)
+
+
+def _make_inverse_cdf_annulus(num_samples=100, eps=1e-4, precision=40):
+    # Precompute the inverse CDF table for annulus return-time sampling
+    # Returns (times, cdf_vals) arrays for use with np.interp
+    mp.mp.dps = precision
+    eps = mp.mpf(eps)
+    times = [10**x for x in np.linspace(-8, 3, num_samples)]
+
+    cdf_vals = np.zeros(num_samples, dtype=float)
+
+    for i, t in enumerate(times):
+        cdf_vals[i] = float(Psi_cdf_annulus(t, eps, method="dehoog"))
+
+    # 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(np.array(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 annular block geometry
+
+    Parameters
+    ----------
+        G : NetworkX graph
+            graph obtained from graph_flow
+
+    Returns
+    -------
+        None
+
+    Notes
+    -----
+        Samples diffusion return times from a finite annular matrix block
+        with one absorbing boundary (fracture-matrix interface) and one
+        reflecting boundary (block interior).
+
+        The dimensionless release position eps = r'/w is fixed at 1e-4,
+        placing the particle just inside the absorbing wall.
+
+        tauD = w^2 / D_matrix sets the diffusion timescale; it is stored
+        on the particle as self.tau_D and used to rescale sampled times.
+
+        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-4
+    b = G.edges[self.curr_node, self.next_node]['b']
+
+    # trapping rate into the 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
+    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()

From 86c049c486ec6f2eea3a8011e56c98d69551737f Mon Sep 17 00:00:00 2001
From: Jeffrey Hyman <jhyman@lanl.gov>
Date: Tue, 10 Mar 2026 17:02:13 -0400
Subject: [PATCH 13/19] update

---
 .gitignore                                    |  2 +-
 .../pydfnworks/dfnGraph/graph_transport.py    | 49 +++++++++----------
 2 files changed, 24 insertions(+), 27 deletions(-)

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_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
index ab9e3187..52b31922 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -20,27 +20,27 @@
 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("-")
@@ -238,17 +238,9 @@ def run_graph_transport(self,
 
     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)
 
+        # build input data list for all particles
+        all_data = []
         for i in range(nparticles):
             data = {}
             data["particle_number"] = i
@@ -263,13 +255,18 @@ 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, verbose),
-                             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))
+
+        with mp.Pool(min(self.ncpu, nparticles)) as pool:
+            for i, particle in enumerate(pool.imap_unordered(
+                    track_particle, all_data, chunksize=chunksize)):
+                particles.append(particle)
+                if i % 1000 == 0:
+                    self.print_log(f"--> Completed {i} out of {nparticles} particles")
 
         elapsed = timeit.default_timer() - tic
         self.print_log(

From 8896662b91ef06a6246b48df2dab88dcbf4833de Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Tue, 10 Mar 2026 17:07:26 -0400
Subject: [PATCH 14/19] some print statements

---
 pydfnworks/pydfnworks/dfnGraph/graph_transport.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
index 52b31922..498ea338 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -237,7 +237,7 @@ def run_graph_transport(self,
             io.dump_trajectories(particles, 1)
 
     if self.ncpu > 1:
-        self.print_log(f"--> Using {self.ncpu} processors")
+        self.print_log(f"--> Building initial data")
 
         # build input data list for all particles
         all_data = []
@@ -261,6 +261,8 @@ def run_graph_transport(self,
         particles = []
         chunksize = max(1, nparticles // (4 * self.ncpu))
 
+        self.print_log("--> Starting main loop")
+        self.print_log(f"--> Using {self.ncpu} processors")
         with mp.Pool(min(self.ncpu, nparticles)) as pool:
             for i, particle in enumerate(pool.imap_unordered(
                     track_particle, all_data, chunksize=chunksize)):

From 1fa29ca02d8a1759fcfdea8587142cdb9dd70d0a Mon Sep 17 00:00:00 2001
From: Jeffrey Hyman <jhyman@lanl.gov>
Date: Wed, 11 Mar 2026 14:41:51 -0400
Subject: [PATCH 15/19] the correct annulus functions now

---
 .../pydfnworks/dfnGraph/particle_class.py     |   4 +-
 .../dfnGraph/tdrw_finite_annulus.py           | 148 +++++++++++++-----
 2 files changed, 114 insertions(+), 38 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/particle_class.py b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
index 4aec9080..4c10bbea 100644
--- a/pydfnworks/pydfnworks/dfnGraph/particle_class.py
+++ b/pydfnworks/pydfnworks/dfnGraph/particle_class.py
@@ -68,7 +68,9 @@ def __init__(self, particle_number, ip, tdrw_flag, tdrw_model,
             self.tau_D = self.fracture_spacing**2 / (self.matrix_diffusivity)
 
         elif tdrw_model == "annulus":
-            # tauD = w^2 / D where w is the annular gap width (fracture_spacing)
+            # 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):
diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
index 5ba20dbc..80986bd2 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
@@ -3,51 +3,117 @@
 import mpmath as mp
 
 
-def Psi_star_annulus(lmbda, eps):
-    # Laplace transform of the CDF for diffusion return time in an annulus
-    # Geometry: absorbing wall at one end, reflecting wall at the other
-    # eps = r'/w, dimensionless release position from the absorbing wall (0 < eps <= 1)
-    # F(s) = [ exp(-eps*q) + exp(-(2-eps)*q) ] / [ (1 + exp(-2q)) * s ]
-    # where q = sqrt(s * tauD);  tauD absorbed into s at call time (s -> s*tauD)
-    sqrt_l = mp.sqrt(lmbda)
-    numerator   = mp.exp(-eps * sqrt_l) + mp.exp(-(2 - eps) * sqrt_l)
-    denominator = (1 + mp.exp(-2 * sqrt_l)) * lmbda
-    return numerator / denominator
-
-
-def Psi_cdf_annulus(t, eps, method="dehoog"):
+def Psi_star_annulus(lmbda, 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).
+    # A particle is released at r' = r0*(1+eps), just inside the absorbing wall.
+    #
+    # Parameters:
+    #   eps  = (r' - r0) / r1   dimensionless release position  (0 < eps << 1)
+    #   tau0 = r0^2 / D         inner-radius diffusion timescale
+    #   tau1 = r1^2 / D         outer-radius diffusion timescale
+    #
+    # Bessel functions: mpmath besseli / besselk
+    #   mpmath.besseli(n, x) = In(x)    (modified Bessel, first kind)
+    #   mpmath.besselk(n, x) = Kn(x)    (modified Bessel, second kind)
+    #
+    # Stability switch: when |s*tau1| > 100 the unscaled Bessels overflow.
+    # Use explicitly scaled forms:
+    #   Ie(n,x) = exp(-x) * In(x)   ->   mp.besseli(n,x) * mp.exp(-x)
+    #   Ke(n,x) = exp(x)  * Kn(x)   ->   mp.besselk(n,x) * mp.exp(x)
+    # with compensating exponential prefactors so the ratio stays finite.
+    #
+    # F_cdf(s) = F_psi(s) / s
+
+    eps  = mp.mpf(eps)
+    tau0 = mp.mpf(tau0)
+    tau1 = mp.mpf(tau1)
+
+    arg0     = mp.sqrt(lmbda * tau0)       # sqrt(s * tau0)
+    arg0_eps = (1 + eps) * arg0            # (1+eps) * sqrt(s * tau0)
+    arg1     = mp.sqrt(lmbda * tau1)       # sqrt(s * tau1)
+
+    if abs(complex(lmbda * tau1)) > 100:
+        # scaled Bessels with explicit exponential prefactors
+        # matches MATLAB besseli(n,x,1) and besselk(n,x,1)
+        I0e_eps = mp.besseli(0, arg0_eps) * mp.exp(-arg0_eps)
+        K0e_eps = mp.besselk(0, arg0_eps) * mp.exp( arg0_eps)
+        I0e_0   = mp.besseli(0, arg0)     * mp.exp(-arg0)
+        K0e_0   = mp.besselk(0, arg0)     * mp.exp( arg0)
+        I1e_1   = mp.besseli(1, arg1)     * mp.exp(-arg1)
+        K1e_1   = mp.besselk(1, arg1)     * mp.exp( arg1)
+
+        # numerator: exp((2+eps)*sqrt(s*tau0) - 2*sqrt(s*tau1)) * I0e_eps * K1e_1
+        #          + exp(-eps*sqrt(s*tau0)) * K0e_eps * I1e_1
+        exp_num1 = mp.exp((2 + eps) * arg0 - 2 * arg1)
+        exp_num2 = mp.exp(-eps * arg0)
+        numerator = exp_num1 * I0e_eps * K1e_1 + exp_num2 * K0e_eps * I1e_1
+
+        # denominator: exp(2*sqrt(s*tau0) - 2*sqrt(s*tau1)) * I0e_0 * K1e_1
+        #            + K0e_0 * I1e_1
+        exp_den     = mp.exp(2 * arg0 - 2 * arg1)
+        denominator = exp_den * I0e_0 * K1e_1 + K0e_0 * I1e_1
+
+    else:
+        # standard unscaled Bessel functions, safe for |s*tau1| <= 100
+        I0_eps = mp.besseli(0, arg0_eps)
+        K0_eps = mp.besselk(0, arg0_eps)
+        I0_0   = mp.besseli(0, arg0)
+        K0_0   = mp.besselk(0, arg0)
+        I1_1   = mp.besseli(1, arg1)
+        K1_1   = mp.besselk(1, arg1)
+
+        numerator   = I0_eps * K1_1 + K0_eps * I1_1
+        denominator = I0_0   * K1_1 + K0_0   * I1_1
+
+    return numerator / (denominator * lmbda)
+
+
+def Psi_pdf_star_annulus(lmbda, eps, tau0, tau1):
+    # Laplace transform of the first-passage time PDF
+    # psi*(s) = s * F_cdf(s)  (removes the 1/s factor)
+    return Psi_star_annulus(lmbda, eps, tau0, tau1) * lmbda
+
+
+def Psi_cdf_annulus(t, eps, tau0, tau1, method="dehoog"):
     # Inverse Laplace transform of Psi_star_annulus at time t
     # Returns cumulative return-time distribution P(T_return <= t)
     if t <= 0:
         return mp.mpf("0.0")
-    F = lambda s: Psi_star_annulus(s, eps)
+    F = lambda s: Psi_star_annulus(s, eps, tau0, tau1)
     return mp.invertlaplace(F, t, method=method)
 
 
-def Psi_pdf_annulus(t, eps, method="dehoog"):
-    # Density psi(t) = d/dt Psi(t), the first-passage time PDF
-    # Inverts the PDF Laplace transform directly for accuracy
-    # psi*(s) = [ exp(-eps*q) + exp(-(2-eps)*q) ] / (1 + exp(-2q))
+def Psi_pdf_annulus(t, eps, tau0, tau1, method="dehoog"):
+    # Inverse Laplace transform of the first-passage time density psi(t)
+    # Inverts psi*(s) directly for accuracy
     if t <= 0:
         return mp.mpf("0.0")
-    psi_star = lambda s: (mp.exp(-eps * mp.sqrt(s)) + mp.exp(-(2 - eps) * mp.sqrt(s))) / (1 + mp.exp(-2 * mp.sqrt(s)))
-    return mp.invertlaplace(psi_star, t, method=method)
+    F = lambda s: Psi_pdf_star_annulus(s, eps, tau0, tau1)
+    return mp.invertlaplace(F, t, method=method)
 
 
-def _make_inverse_cdf_annulus(num_samples=100, eps=1e-4, precision=40):
-    # Precompute the inverse CDF table for annulus return-time sampling
+def _make_inverse_cdf_annulus(num_samples=100, eps=1e-2, tau0=1e-4, tau1=1e0, precision=40):
+    # Precompute the inverse CDF table for cylindrical annulus return-time sampling
+    #
+    # Default parameters match InvLaplace.m:
+    #   eps  = 1e-2   dimensionless release position (r'-r0)/r1
+    #   tau0 = 1e-4   inner-radius timescale r0^2/D  (dimensionless, tau0/tau1)
+    #   tau1 = 1e0    outer-radius timescale r1^2/D  (reference, = 1)
+    #
     # Returns (times, cdf_vals) arrays for use with np.interp
     mp.mp.dps = precision
-    eps = mp.mpf(eps)
-    times = [10**x for x in np.linspace(-8, 3, num_samples)]
 
+    times    = [10**x for x in np.linspace(-8, 3, num_samples)]
     cdf_vals = np.zeros(num_samples, dtype=float)
 
     for i, t in enumerate(times):
-        cdf_vals[i] = float(Psi_cdf_annulus(t, eps, method="dehoog"))
+        cdf_vals[i] = float(Psi_cdf_annulus(t, eps, tau0, tau1, method="dehoog"))
 
     # sort by cdf value to build the inverse CDF (cdf -> time) lookup
-    order = np.argsort(cdf_vals)
+    order      = np.argsort(cdf_vals)
     cdf_sorted = np.clip(cdf_vals[order], 0, 1)
     t_sorted   = np.maximum(np.array(times)[order], 0)
 
@@ -59,7 +125,7 @@ def _make_inverse_cdf_annulus(num_samples=100, eps=1e-4, precision=40):
 
 
 def limited_matrix_diffusion_annulus(self, G):
-    """ Matrix diffusion with finite annular block geometry
+    """ Matrix diffusion with finite cylindrical annular block geometry
 
     Parameters
     ----------
@@ -72,15 +138,22 @@ def limited_matrix_diffusion_annulus(self, G):
 
     Notes
     -----
-        Samples diffusion return times from a finite annular matrix block
-        with one absorbing boundary (fracture-matrix interface) and one
-        reflecting boundary (block interior).
+        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 from the cylindrical geometry, unlike the slab (Dentz) model
+        which uses hyperbolic functions.
 
-        The dimensionless release position eps = r'/w is fixed at 1e-4,
-        placing the particle just inside the absorbing wall.
+        Two diffusion timescales parameterize the geometry:
+            tau0 = r0^2 / D   (inner radius timescale)
+            tau1 = r1^2 / D   (outer radius timescale)
 
-        tauD = w^2 / D_matrix sets the diffusion timescale; it is stored
-        on the particle as self.tau_D and used to rescale sampled times.
+        tau1 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
@@ -89,10 +162,10 @@ def limited_matrix_diffusion_annulus(self, G):
         All other parameters are attached to the particle class.
     """
 
-    eps = 1e-4
+    eps = 1e-2
     b = G.edges[self.curr_node, self.next_node]['b']
 
-    # trapping rate into the matrix block
+    # 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
@@ -102,7 +175,8 @@ def limited_matrix_diffusion_annulus(self, G):
     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)
+    # 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()

From e0547cbf8e861c7eff17f4c3817dd5b2aefa35ec Mon Sep 17 00:00:00 2001
From: Jeffrey Hyman <jhyman@lanl.gov>
Date: Mon, 16 Mar 2026 15:46:03 -0400
Subject: [PATCH 16/19] new form

---
 .../dfnGraph/tdrw_finite_annulus.py           | 181 ++++++++----------
 1 file changed, 85 insertions(+), 96 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
index 80986bd2..4fe7ae6c 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
@@ -1,9 +1,28 @@
 
 import numpy as np
-import mpmath as mp
-
-
-def Psi_star_annulus(lmbda, eps, tau0, tau1):
+from scipy.special import ive, kve
+from math import factorial
+
+
+def _stehfest_coefficients(N=16):
+    # Compute Stehfest coefficients V_i for i = 1 ... N
+    # N must be even. N=16 recommended for double precision.
+    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
@@ -11,111 +30,83 @@ def Psi_star_annulus(lmbda, eps, tau0, tau1):
     # A particle is released at r' = r0*(1+eps), just inside the absorbing wall.
     #
     # Parameters:
-    #   eps  = (r' - r0) / r1   dimensionless release position  (0 < eps << 1)
-    #   tau0 = r0^2 / D         inner-radius diffusion timescale
-    #   tau1 = r1^2 / D         outer-radius diffusion timescale
+    #   eps  = (r' - r0) / r1   dimensionless release position  (0 < eps < 1)
+    #   tau0 = r0^2 / D         inner-radius timescale ratio r0^2/r1^2
+    #   tau1 = r1^2 / D         outer-radius timescale (reference = 1)
     #
-    # Bessel functions: mpmath besseli / besselk
-    #   mpmath.besseli(n, x) = In(x)    (modified Bessel, first kind)
-    #   mpmath.besselk(n, x) = Kn(x)    (modified Bessel, second kind)
+    # Uses scipy scaled Bessel functions throughout:
+    #   ive(n, x) = exp(-x) * In(x)   matches MATLAB besseli(n, x, 1)
+    #   kve(n, x) = exp(x)  * Kn(x)   matches MATLAB besselk(n, x, 1)
     #
-    # Stability switch: when |s*tau1| > 100 the unscaled Bessels overflow.
-    # Use explicitly scaled forms:
-    #   Ie(n,x) = exp(-x) * In(x)   ->   mp.besseli(n,x) * mp.exp(-x)
-    #   Ke(n,x) = exp(x)  * Kn(x)   ->   mp.besselk(n,x) * mp.exp(x)
-    # with compensating exponential prefactors so the ratio stays finite.
+    # Exponential compensation factors normalize by exp(2*q0 - 2*q1) so
+    # the ratio goes correctly to 0 as s -> inf:
+    #   F_psi(s->0)   = 1
+    #   F_psi(s->inf) = 0
     #
     # F_cdf(s) = F_psi(s) / s
 
-    eps  = mp.mpf(eps)
-    tau0 = mp.mpf(tau0)
-    tau1 = mp.mpf(tau1)
-
-    arg0     = mp.sqrt(lmbda * tau0)       # sqrt(s * tau0)
-    arg0_eps = (1 + eps) * arg0            # (1+eps) * sqrt(s * tau0)
-    arg1     = mp.sqrt(lmbda * tau1)       # sqrt(s * tau1)
-
-    if abs(complex(lmbda * tau1)) > 100:
-        # scaled Bessels with explicit exponential prefactors
-        # matches MATLAB besseli(n,x,1) and besselk(n,x,1)
-        I0e_eps = mp.besseli(0, arg0_eps) * mp.exp(-arg0_eps)
-        K0e_eps = mp.besselk(0, arg0_eps) * mp.exp( arg0_eps)
-        I0e_0   = mp.besseli(0, arg0)     * mp.exp(-arg0)
-        K0e_0   = mp.besselk(0, arg0)     * mp.exp( arg0)
-        I1e_1   = mp.besseli(1, arg1)     * mp.exp(-arg1)
-        K1e_1   = mp.besselk(1, arg1)     * mp.exp( arg1)
-
-        # numerator: exp((2+eps)*sqrt(s*tau0) - 2*sqrt(s*tau1)) * I0e_eps * K1e_1
-        #          + exp(-eps*sqrt(s*tau0)) * K0e_eps * I1e_1
-        exp_num1 = mp.exp((2 + eps) * arg0 - 2 * arg1)
-        exp_num2 = mp.exp(-eps * arg0)
-        numerator = exp_num1 * I0e_eps * K1e_1 + exp_num2 * K0e_eps * I1e_1
-
-        # denominator: exp(2*sqrt(s*tau0) - 2*sqrt(s*tau1)) * I0e_0 * K1e_1
-        #            + K0e_0 * I1e_1
-        exp_den     = mp.exp(2 * arg0 - 2 * arg1)
-        denominator = exp_den * I0e_0 * K1e_1 + K0e_0 * I1e_1
-
-    else:
-        # standard unscaled Bessel functions, safe for |s*tau1| <= 100
-        I0_eps = mp.besseli(0, arg0_eps)
-        K0_eps = mp.besselk(0, arg0_eps)
-        I0_0   = mp.besseli(0, arg0)
-        K0_0   = mp.besselk(0, arg0)
-        I1_1   = mp.besseli(1, arg1)
-        K1_1   = mp.besselk(1, arg1)
-
-        numerator   = I0_eps * K1_1 + K0_eps * I1_1
-        denominator = I0_0   * K1_1 + K0_0   * I1_1
-
-    return numerator / (denominator * lmbda)
-
-
-def Psi_pdf_star_annulus(lmbda, eps, tau0, tau1):
-    # Laplace transform of the first-passage time PDF
-    # psi*(s) = s * F_cdf(s)  (removes the 1/s factor)
-    return Psi_star_annulus(lmbda, eps, tau0, tau1) * lmbda
+    q0     = np.sqrt(s * tau0)
+    q0_eps = (1.0 + eps) * q0
+    q1     = np.sqrt(s * tau1)
 
+    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)
 
-def Psi_cdf_annulus(t, eps, tau0, tau1, method="dehoog"):
-    # Inverse Laplace transform of Psi_star_annulus at time t
-    # Returns cumulative return-time distribution P(T_return <= t)
-    if t <= 0:
-        return mp.mpf("0.0")
-    F = lambda s: Psi_star_annulus(s, eps, tau0, tau1)
-    return mp.invertlaplace(F, t, method=method)
+    # exponential compensation factors, clipped to avoid overflow
+    e1 = np.exp(np.clip( eps * q0,               -500, 500))
+    e2 = np.exp(np.clip(-eps * q0 - 2*q0 + 2*q1, -500, 500))
+    e3 = np.exp(np.clip(-(2*q0 - 2*q1),           -500, 500))
 
+    numerator   = e1 * I0e_eps * K1e_1 + e2 * K0e_eps * I1e_1
+    denominator = I0e_0 * K1e_1        + e3 * K0e_0   * I1e_1
 
-def Psi_pdf_annulus(t, eps, tau0, tau1, method="dehoog"):
-    # Inverse Laplace transform of the first-passage time density psi(t)
-    # Inverts psi*(s) directly for accuracy
-    if t <= 0:
-        return mp.mpf("0.0")
-    F = lambda s: Psi_pdf_star_annulus(s, eps, tau0, tau1)
-    return mp.invertlaplace(F, t, method=method)
+    return numerator / (denominator * 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 _make_inverse_cdf_annulus(num_samples=100, eps=1e-2, tau0=1e-4, tau1=1e0, precision=40):
+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=0.5, tau0=0.5, tau1=1.0, stehfest_n=16):
     # Precompute the inverse CDF table for cylindrical annulus return-time sampling
     #
-    # Default parameters match InvLaplace.m:
-    #   eps  = 1e-2   dimensionless release position (r'-r0)/r1
-    #   tau0 = 1e-4   inner-radius timescale r0^2/D  (dimensionless, tau0/tau1)
-    #   tau1 = 1e0    outer-radius timescale r1^2/D  (reference, = 1)
+    # Parameters:
+    #   eps  = (r'-r0)/r1   dimensionless release position  (0 < eps < 1)
+    #   tau0 = r0^2/r1^2    inner-radius timescale ratio
+    #   tau1 = 1.0          outer-radius reference timescale
+    #
+    # Defaults are physically meaningful values that produce a well-behaved CDF.
+    # Update once Marco Dentz confirms the physical parameter mapping.
     #
     # Returns (times, cdf_vals) arrays for use with np.interp
-    mp.mp.dps = precision
 
-    times    = [10**x for x in np.linspace(-8, 3, num_samples)]
-    cdf_vals = np.zeros(num_samples, dtype=float)
+    V     = _stehfest_coefficients(stehfest_n)
+    times = np.logspace(-3, 3, num_samples)
 
-    for i, t in enumerate(times):
-        cdf_vals[i] = float(Psi_cdf_annulus(t, eps, tau0, tau1, method="dehoog"))
+    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(np.array(times)[order], 0)
+    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)
@@ -146,14 +137,12 @@ def limited_matrix_diffusion_annulus(self, G):
 
         The Laplace-domain solution uses modified Bessel functions I0, I1,
         K0, K1 from the cylindrical geometry, unlike the slab (Dentz) model
-        which uses hyperbolic functions.
-
-        Two diffusion timescales parameterize the geometry:
-            tau0 = r0^2 / D   (inner radius timescale)
-            tau1 = r1^2 / D   (outer radius timescale)
+        which uses hyperbolic functions. Scaled scipy Bessel functions
+        (ive, kve) are used with explicit exponential compensation factors
+        to ensure numerical stability for all parameter values.
 
-        tau1 is stored on the particle as self.tau_D and used to rescale
-        dimensionless sampled return times to physical times [s].
+        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
@@ -162,7 +151,7 @@ def limited_matrix_diffusion_annulus(self, G):
         All other parameters are attached to the particle class.
     """
 
-    eps = 1e-2
+    eps = 0.5
     b = G.edges[self.curr_node, self.next_node]['b']
 
     # trapping rate into the cylindrical matrix block

From fdfc8811d805051f2f155cf95e97787e9355d064 Mon Sep 17 00:00:00 2001
From: Jeffrey Hyman <jhyman@lanl.gov>
Date: Tue, 17 Mar 2026 16:15:00 -0400
Subject: [PATCH 17/19] new annulus code

---
 .../dfnGraph/tdrw_finite_annulus.py           | 89 ++++++++++---------
 1 file changed, 47 insertions(+), 42 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
index 4fe7ae6c..abb14a93 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_annulus.py
@@ -1,12 +1,10 @@
-
 import numpy as np
-from scipy.special import ive, kve
+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. N=16 recommended for double precision.
+    # 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
@@ -27,44 +25,52 @@ def Psi_star_annulus(s, eps, tau0, tau1):
     #
     # Geometry: cylindrical annulus with inner radius r0 (absorbing, fracture-matrix
     # interface) and outer radius r1 (reflecting, block interior).
-    # A particle is released at r' = r0*(1+eps), just inside the absorbing wall.
+    # Particle released at r' = r0*(1+eps), just inside the absorbing wall.
     #
     # Parameters:
-    #   eps  = (r' - r0) / r1   dimensionless release position  (0 < eps < 1)
-    #   tau0 = r0^2 / D         inner-radius timescale ratio r0^2/r1^2
-    #   tau1 = r1^2 / D         outer-radius timescale (reference = 1)
+    #   eps  = (r' - r0) / r1   dimensionless release position
+    #   tau0 = r0^2 / D         inner-radius diffusion timescale
+    #   tau1 = r1^2 / D         outer-radius diffusion timescale
     #
-    # Uses scipy scaled Bessel functions throughout:
-    #   ive(n, x) = exp(-x) * In(x)   matches MATLAB besseli(n, x, 1)
-    #   kve(n, x) = exp(x)  * Kn(x)   matches MATLAB besselk(n, x, 1)
+    # 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)
     #
-    # Exponential compensation factors normalize by exp(2*q0 - 2*q1) so
-    # the ratio goes correctly to 0 as s -> inf:
-    #   F_psi(s->0)   = 1
-    #   F_psi(s->inf) = 0
+    #   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)
 
-    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)
+    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))
 
-    # exponential compensation factors, clipped to avoid overflow
-    e1 = np.exp(np.clip( eps * q0,               -500, 500))
-    e2 = np.exp(np.clip(-eps * q0 - 2*q0 + 2*q1, -500, 500))
-    e3 = 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
 
-    numerator   = e1 * I0e_eps * K1e_1 + e2 * K0e_eps * I1e_1
-    denominator = I0e_0 * K1e_1        + e3 * 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)
 
-    return numerator / (denominator * s)
+    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):
@@ -84,21 +90,21 @@ def _stehfest_invert(F, t_arr, V, **kwargs):
     return f
 
 
-def _make_inverse_cdf_annulus(num_samples=100, eps=0.5, tau0=0.5, tau1=1.0, stehfest_n=16):
+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:
-    #   eps  = (r'-r0)/r1   dimensionless release position  (0 < eps < 1)
-    #   tau0 = r0^2/r1^2    inner-radius timescale ratio
-    #   tau1 = 1.0          outer-radius reference timescale
+    # 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
     #
-    # Defaults are physically meaningful values that produce a well-behaved CDF.
-    # Update once Marco Dentz confirms the physical parameter mapping.
+    # 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(-3, 3, num_samples)
+    times = np.logspace(-10, 2, num_samples)
 
     cdf_vals = _stehfest_invert(Psi_star_annulus, times, V,
                                 eps=eps, tau0=tau0, tau1=tau1)
@@ -136,10 +142,9 @@ def limited_matrix_diffusion_annulus(self, G):
         absorbing wall.
 
         The Laplace-domain solution uses modified Bessel functions I0, I1,
-        K0, K1 from the cylindrical geometry, unlike the slab (Dentz) model
-        which uses hyperbolic functions. Scaled scipy Bessel functions
-        (ive, kve) are used with explicit exponential compensation factors
-        to ensure numerical stability for all parameter values.
+        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].
@@ -151,7 +156,7 @@ def limited_matrix_diffusion_annulus(self, G):
         All other parameters are attached to the particle class.
     """
 
-    eps = 0.5
+    eps = 1e-2
     b = G.edges[self.curr_node, self.next_node]['b']
 
     # trapping rate into the cylindrical matrix block
@@ -168,4 +173,4 @@ def limited_matrix_diffusion_annulus(self, G):
     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()
+    self.delta_t_md = tmp.sum()
\ No newline at end of file

From a2b1c813f2e145648ab2914efb2ee9bc0705dd23 Mon Sep 17 00:00:00 2001
From: Jeffrey Hyman <jhyman@lanl.gov>
Date: Tue, 17 Mar 2026 16:37:25 -0400
Subject: [PATCH 18/19] fixed inv laplace issue

---
 .../pydfnworks/dfnGraph/tdrw_finite_dentz.py  | 164 +++++++++---------
 1 file changed, 84 insertions(+), 80 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
index 489d318b..621dbf12 100644
--- a/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
+++ b/pydfnworks/pydfnworks/dfnGraph/tdrw_finite_dentz.py
@@ -1,63 +1,84 @@
-
-
-from scipy.interpolate import PchipInterpolator
 import numpy as np
-import matplotlib.pyplot as plt
-import mpmath as mp
-
-
-def Psi_star(lmbda, eps):
-    # Laplace transform of the CDF Ψ*_ε(λ)
-    # eps must satisfy 0 < eps < 1 for a slab of unit width in the note
-    sqrt_l = mp.sqrt(lmbda)
-    return mp.cosh((1 - eps) * sqrt_l) / (lmbda * mp.cosh(sqrt_l))
-
-def Psi_cdf(t, eps, method="dehoog"):
-    # Inverse Laplace transform of Ψ*_ε at time t
-    # method can be "dehoog" or "talbot"
-    if t <= 0:
-        return mp.mpf("0.0")
-    F = lambda s: Psi_star(s, eps)
-    return mp.invertlaplace(F, t, method=method)
-
-def Psi_pdf(t, eps, method="dehoog", h=None):
-    # Density ψ_ε(t) = d/dt Ψ_ε(t)
-    # You can either invert ψ*_ε directly or differentiate the CDF numerically
-    # This version inverts ψ*_ε for accuracy
-    if t <= 0:
-        return mp.mpf("0.0")
-    psi_star = lambda s: mp.cosh((1 - eps) * mp.sqrt(s)) / mp.cosh(mp.sqrt(s))
-    return mp.invertlaplace(psi_star, t, method=method)
-
-
-def _make_inverse_cdf_spline_for_times(num_samples = 100, eps = 1e-4, precision = 40):
-    # --- compute CDF and PDF values ---
-    mp.mp.dps = precision
-    num_samples = num_samples
-    eps = mp.mpf(eps)
-    times = [10**x for x in np.linspace(-8, 3, num_samples)]
-
-    cdf_vals = np.zeros(num_samples, dtype=float)
-
-    for i, t in enumerate(times):
-        cdf_vals[i] = float(Psi_cdf(t, eps, method="dehoog"))
-
-    # --- build inverse CDF spline (t vs cdf) ---
-    order = np.argsort(cdf_vals)
+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(np.array(times)[order], 0)
+    t_sorted   = np.maximum(times[order], 0)
 
-    # remove duplicates
+    # 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])
 
-    finite_md_times = t_unique
-    finite_md_cdf = cdf_unique
+    return t_unique, cdf_unique
 
-    return finite_md_times, finite_md_cdf 
 
-
-def limited_matrix_diffusion_dentz(self,G):
+def limited_matrix_diffusion_dentz(self, G):
     """ Matrix diffusion with limited block size
 
     Parameters
@@ -74,37 +95,20 @@ def limited_matrix_diffusion_dentz(self,G):
         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
+
+    # 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
-    #print(f"average_number_of_trapping_events: {average_number_of_trapping_events}")
-    # number of trapping events in sampled from a poisson distribution
+
+    # sample number of trapping events from 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)
 
+    # 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

From 32730087c08e0a52a83b0c653f600db2d7dc0e7f Mon Sep 17 00:00:00 2001
From: Jeffrey De'Haven Hyman <jhyman@lanl.gov>
Date: Tue, 17 Mar 2026 16:49:44 -0400
Subject: [PATCH 19/19] new counter for graph transport

---
 pydfnworks/pydfnworks/dfnGraph/graph_transport.py | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
index 498ea338..3cd1a3f6 100644
--- a/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
+++ b/pydfnworks/pydfnworks/dfnGraph/graph_transport.py
@@ -263,12 +263,19 @@ def run_graph_transport(self,
 
         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)
-                if i % 1000 == 0:
-                    self.print_log(f"--> Completed {i} out of {nparticles} particles")
+                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(