### Install Dependencies and Run Tests

Source: https://github.com/tldr-group/taufactor/blob/main/CONTRIBUTING.md

Install the project in editable mode with development dependencies and run linters and tests. This verifies the setup and checks for code quality issues.

```bash
pip install -e .[dev]
ruff check .
pytest
```

--------------------------------

### Setup TauFactor Development Environment

Source: https://github.com/tldr-group/taufactor/blob/main/docs/installation.md

Create and activate a Conda environment using the provided environment.yml file, then install TauFactor in editable mode.

```default
conda env create -f environment.yml
conda activate taufactor
pip install -e .
```

--------------------------------

### Install PyTorch with CUDA

Source: https://github.com/tldr-group/taufactor/blob/main/README.md

Use this command to install PyTorch with CUDA support on a Linux machine. Ensure you have the correct PyTorch version installed.

```bash
conda install pytorch pytorch-cuda=11.7 -c pytorch -c nvidia
```

--------------------------------

### Install TauFactor via PyPI

Source: https://github.com/tldr-group/taufactor/blob/main/README.md

Install the TauFactor package using pip. This command should be run in your Python environment.

```bash
pip install taufactor
```

--------------------------------

### Run Unit Tests

Source: https://github.com/tldr-group/taufactor/blob/main/README.md

Navigate to the root directory of the project and run this command to execute the unit tests.

```bash
pytest
```

--------------------------------

### Initialize and Solve with Multi-Phase Solver

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/00-overview.ipynb

Initializes the `MultiPhaseSolver` with image data and conductivity values for multiple phases, setting an iteration limit before solving.

```python
# assign conductivity values, where key is segmented label in 'img'
# and value is conductivity
cond = {1:0.32, 2:0.44}

# create a multiphase solver object and set an iteration limit
s = tau.MultiPhaseSolver(img, cond=cond, iter_limit=1000)
# call solve function
s.solve()
```

--------------------------------

### Initialize and Solve with Periodic Solver

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/00-overview.ipynb

Initializes the `PeriodicSolver` with the image data and calls the solve method to compute tortuosity using periodic boundary conditions.

```python
s = tau.PeriodicSolver(img==0)
s.solve()
```

--------------------------------

### Visualize Flux Direction

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/00-overview.ipynb

Uses the utility function `flux_direction` to visualize the direction of transport in the image data, aiding in understanding the simulation setup.

```python
from taufactor.utils import flux_direction
figure = flux_direction(img)
```

--------------------------------

### Initialize and Solve with PeriodicElectrodeSolver (Default Convergence)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Initializes the PeriodicElectrodeSolver with the microstructure batch and solves for tortuosity using default convergence criteria. The spacing argument is crucial for correctly scaling the specific surface area.

```python
s = tau.PeriodicElectrodeSolver(batch, device='cuda', spacing=Lx/nx)
tau_pore = s.solve(iter_limit=50000, conv_crit=1e-3)
```

--------------------------------

### Calculate Raw Interfacial Areas for Multi-Label Structure

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Computes the total amount of counted faces for interfacial areas in a multi-label structure. Use `normalize=True` to get specific interfacial areas.

```python
tau.interfacial_areas(cubes, normalize = False)
```

--------------------------------

### Create and Activate Development Branch

Source: https://github.com/tldr-group/taufactor/blob/main/CONTRIBUTING.md

Navigate into the cloned directory and create a new branch for your development work. This isolates your changes.

```bash
cd taufactor
git checkout -b my_dev_branch
```

--------------------------------

### Solve for Tau using PeriodicSolver (Classic)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Initializes and solves for tau using the standard PeriodicSolver. This is a classic approach for comparison. Note the potential for non-convergence indicated by the output.

```python
t = tau.PeriodicSolver(batch, device='cuda')
tau_classic = t.solve(iter_limit=300000, conv_crit=1e-6, verbose='debug', plot_interval=100)
```

--------------------------------

### Load EIS Simulation Results

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Loads pre-computed Electrochemical Impedance Spectroscopy (EIS) simulation results from pickle files. Assumes files are named 'impedance_solve_[structure].pkl'.

```python
import pickle

Z_EIS = {}
for structure in structures:
    filename = "impedance_solve_" + structure + ".pkl"
    with open(filename, "rb") as f:
        data = pickle.load(f)
        Z_EIS[structure] = data.impedance
```

--------------------------------

### Stage, Commit, and Push Changes

Source: https://github.com/tldr-group/taufactor/blob/main/CONTRIBUTING.md

Add your modified files, commit them with a descriptive message, and push the branch to your remote fork. Include tests for new features.

```bash
git add my_changed_files
git commit -m "Detailed description of changes."
git push my_remote_fork my_dev_branch
```

--------------------------------

### Version Bump and Deployment

Source: https://github.com/tldr-group/taufactor/blob/main/CONTRIBUTING.md

Update the project version using bump2version and push changes and tags. Travis CI will handle deployment to PyPI if tests pass.

```bash
bump2version patch # possible: major / minor / patch
git push
git push --tags
```

--------------------------------

### Create and Activate Conda Environment

Source: https://github.com/tldr-group/taufactor/blob/main/CONTRIBUTING.md

Set up a dedicated Python environment for development using Conda. This ensures dependency isolation.

```bash
conda create --name myenv python=3.12
conda activate myenv
```

--------------------------------

### Simulate and Save Impedance (smooth-rough)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Simulates impedance for the 'smooth-rough' configuration in 'nyquist' mode. The solver is set with a standard iteration limit, and the resulting object is saved.

```python
e = tau.PeriodicImpedanceSolver(fields["smooth-rough"], mode='nyquist')
e.solve(iter_limit=500000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(e, "impedance_solve_smooth-rough")
```

--------------------------------

### Solve Tortuosity with Plotting and Custom Convergence

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/00-overview.ipynb

Initializes the solver and solves for tortuosity, enabling plotting of flux convergence and setting a custom convergence criterion. The output displays the convergence error and tau value at the final iteration.

```python
s = tau.Solver(img==0)
tau_x = s.solve(verbose='plot', plot_interval=1, conv_crit=1e-3)
```

--------------------------------

### Compare Performance of Specific Surface Area Calculation Methods

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Measures and prints the wall time for different specific surface area calculation methods ('face_counting', 'gradient', 'marching_cubes') with varying parameters ('smoothing', 'device'). Useful for understanding computational cost.

```python
import time
timer = []
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='face_counting', device='cuda'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='gradient', smoothing=False, device='cuda'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='gradient', smoothing=True, device='cuda'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='marching_cubes', smoothing=False, device='cpu'))
timer.append(time.time())
print(tau.specific_surface_area(structure, phases=labels, method='marching_cubes', smoothing=True, device='cpu'))
timer.append(time.time())

print(f"Face_counting: {timer[1]-timer[0]:.4f} s\n"
```

```python
      f"Gradient w/o smoothing: {timer[2]-timer[1]:.4f} s\n"
      f"Gradient w smoothing: {timer[3]-timer[2]:.4f} s\n"
      f"Marching w/o smoothing: {timer[4]-timer[3]:.4f} s\n"
      f"Marching w smoothing: {timer[5]-timer[4]:.4f} s"
```

--------------------------------

### Visualize Multi-Feature and Labeled Structures

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Creates test structures with multiple features and applies periodic labeling using `scipy.ndimage.generate_binary_structure` and `taufactor.metrics.label_periodic`. This is useful for benchmarking more complex scenarios and validating solver consistency.

```python
from taufactor.metrics import label_periodic
from scipy.ndimage import generate_binary_structure

fields = {}
features = 2
fields["blocks"] = create_stacked_blocks(Nx, features=features)
fields["2d_diagonals"] = create_2d_diagonals(Nx, features=features)
fields["2d_zigzag"] = create_2d_zigzag(Nx, features=features)
fields["3d_diagonals"] = create_3d_diagonals(Nx, features=features)

fields["blocks"][Nx//2:,:Nx//2,Nx//2:] *= 2
connectivity = generate_binary_structure(3, 1)
fields["2d_diagonals"] = label_periodic(fields["2d_diagonals"], 1, connectivity, periodic=(False, True, True))[0]
fields["2d_zigzag"] = label_periodic(fields["2d_zigzag"], 1, connectivity, periodic=(False, True, True))[0]
fields["3d_diagonals"] = label_periodic(fields["3d_diagonals"], 1, connectivity, periodic=(False, True, True))[0]

colors = ['white', 'gray', 'red', 'blue', 'green']
fig, axes = plt.subplots(nrows=1, ncols=4, figsize=(15, 4), dpi=100)
for i, structure in enumerate(fields.keys()):
    axes[i].remove()
    axes[i] = fig.add_subplot(1, 5, i + 1, projection='3d')
    for j in np.unique(fields[structure]):
        if j>0:
            axes[i].voxels(np.transpose(fields[structure], (2,1,0))==j, facecolors=colors[j], edgecolor='black', alpha=1.0)
    axes[i].set_zlim(0, Nx)
    axes[i].set_aspect('equal')
    axes[i].set_title(structure)
    axes[i].set_xlim(0, Nx)
    axes[i].set_ylim(0, Nx)
plt.tight_layout()
plt.show()
```

--------------------------------

### Running Benchmark with Custom Multi-Phase Structure

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Executes a benchmark study using the 'MultiPhaseSolver' with a custom 2D zigzag structure and varying diffusivities. The results are collected and stored in a DataFrame.

```python
from taufactor.utils import create_2d_zigzag

def multi_zigzag(Nx):
    geom = create_2d_zigzag(Nx, features=1)
    N2 = Nx // 2
    geom[:N2, :, :][geom[:N2, :, :] == 1] = 2
    return geom

all_rows = []
diffusivities = [0.1, 0.5, 1.0, 2.0, 10.0]

for c in diffusivities:
    rows = run_benchmark_study(
        Ns=[64],
        solver='MultiPhaseSolver',
        solver_kwargs={'diffusivities': {0:0, 1: 1.0, 2: c}},
        structure=multi_zigzag,
        write_file=False,
    )
    for r in rows:
        r['diffusivity'] = c
    all_rows.extend(rows)

df_ratio = pd.DataFrame(all_rows).sort_values('diffusivity').reset_index(drop=True)
df_ratio
```

--------------------------------

### Compare Discrete and Analytical Impedance Solutions

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/11-transmission-line-model.ipynb

Compares the impedance calculated using a discrete transmission line model with an analytical solution. It visualizes the results on a Nyquist plot and prints the percentage deviation at a reference frequency.

```python
L = 2
conductivity = 0.45 # [S/m] = [1/Ohm m]
capacitance = 0.19  # [F/m^2] = [C/V m^2]
area = 1 # cross-sectional area [m^2]
surf = 0.11 # specific surface area [1/m]

# Resistivity and capacitance per length 
r_x = 1/conductivity/area   # [Ohm / m]
c_x = capacitance*surf*area # [F / m]

freq_0 = 1 / r_x / c_x / L**2
freq = freq_0 * 2 ** np.arange(-3, 10, 0.5)
lambdas_0 = np.sqrt(1j * freq * r_x * c_x)
Z_pde_ref = r_x * np.cosh(lambdas_0*L) / np.sinh(lambdas_0*L) / lambdas_0

fig, ax = plt.subplots(figsize=(3, 4), dpi=100)
colors = ['b', 'g', 'r']
for i, N_el in enumerate([20, 50, 100]):
    # Discrete resistance and capacitance
    R_i = r_x * L / N_el
    C_i = c_x * L / N_el

    Z_mean = compute_impedance(R_i*np.ones(N_el), C_i*np.ones(N_el), freq)
    print(f"Deviation of {100*(Z_pde_ref[0].real-Z_mean[0].real)/Z_mean[0].real:.4f} %")
    plt.plot(Z_mean.real, -Z_mean.imag, '^-', color=colors[i], label=f'N={N_el}')
    ref_idx = np.argmin(np.abs(freq - freq_0))
    plt.scatter(Z_mean.real[ref_idx], -Z_mean.imag[ref_idx], c=colors[i], edgecolors='k', s=80)
ax.plot(Z_pde_ref.real, -Z_pde_ref.imag, label='Analytical', color='gray', linestyle='--')
plt.xlabel("Z' [Ohm]")
plt.ylabel("-Z'' [Ohm]")
# plt.title("Nyquist Plot")
# ax.set_aspect('equal')
plt.legend()
plt.grid()
plt.xlim([0, 1.2*Z_mean[0].real])
plt.ylim([0, 3*1.1*Z_mean[0].real])
plt.show()
```

--------------------------------

### Initialize and Solve with Anisotropic Solver

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/00-overview.ipynb

Initializes the `AnisotropicSolver` with image data and spacing information to account for non-cubic voxels, then solves for tortuosity.

```python
s = tau.AnisotropicSolver(img==0, spacing=(1,1,2))
s.solve()
```

--------------------------------

### Simulate and Save Impedance (wedge)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Executes an impedance simulation for the 'wedge' configuration with a high iteration limit. The simulation runs in 'nyquist' mode, and the solver object is persisted.

```python
f = tau.PeriodicImpedanceSolver(fields["wedge"], mode='nyquist')
f.solve(iter_limit=1000000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(f, "impedance_solve_wedge")
```

--------------------------------

### Running Benchmark with Periodic Multi-Phase Solver and Multi-Diagonal Structure

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Conducts a benchmark study using the 'PeriodicMultiPhaseSolver' with a custom multi-diagonal structure and varying diffusivities. Results are aggregated into a DataFrame.

```python
from taufactor.utils import create_2d_diagonals

def multi_diagonal(Nx):
    geom = create_2d_diagonals(Nx, features=2)
    geom = label_periodic(geom, 1, connectivity, periodic=(False, True, True))[0]
    return geom

all_rows = []
diffusivities = [0.1, 0.5, 1.0, 2.0, 10.0]

for c in diffusivities:
    rows = run_benchmark_study(
        Ns=[64],
        solver='PeriodicMultiPhaseSolver',
        solver_kwargs={'diffusivities': {0:0, 1:1, 2:c, 3:1, 4:c}},
        structure=multi_diagonal,
        write_file=False,
    )
    for r in rows:
        r['diffusivity'] = c
    all_rows.extend(rows)

df_ratio = pd.DataFrame(all_rows).sort_values('diffusivity').reset_index(drop=True)
df_ratio
```

--------------------------------

### Simulate and Save Impedance (Nyquist Mode)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Initializes and solves a `PeriodicImpedanceSolver` in 'nyquist' mode for 'diagonals' configuration. The simulation results are then saved using the `save_solver_object` function. Adjust `iter_limit` as needed for convergence.

```python
a = tau.PeriodicImpedanceSolver(fields['diagonals'], mode='nyquist')
a.solve(iter_limit=500000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(a, "impedance_solve_diagonals")
```

--------------------------------

### Simulate and Save Impedance (diag-pillars)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Solves the impedance for the 'diag-pillars' configuration using 'nyquist' mode. The solver is configured with a specific iteration limit and verbosity settings. Results are saved to a pickle file.

```python
b = tau.PeriodicImpedanceSolver(fields["diag-pillars"], mode='nyquist')
b.solve(iter_limit=500000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(b, "impedance_solve_diag-pillars")
```

--------------------------------

### Check Changes and Run Tests

Source: https://github.com/tldr-group/taufactor/blob/main/CONTRIBUTING.md

After making changes, use Ruff to check for code style issues and Pytest to ensure all tests pass. This should be done before committing.

```bash
ruff check .
pytest
```

--------------------------------

### Simulate and Save Impedance (small-large)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Runs an impedance simulation for the 'small-large' configuration, similar to the 'large-small' case but with potentially different field definitions. The simulation uses 'nyquist' mode and saves the results.

```python
d = tau.PeriodicImpedanceSolver(fields["small-large"], mode='nyquist')
d.solve(iter_limit=1000000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(d, "impedance_solve_small-large")
```

--------------------------------

### Solve Tortuosity and Effective Diffusivity

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/00-overview.ipynb

Loads a segmented image, initializes the solver, and calculates the tortuosity factor and effective diffusivity. The output shows the converged values and computation time.

```python
import taufactor as tau
import tifffile

# Load segmented image, in this case with
# labels {"pore":0, "NMC":85, "CBD":170}
img = tifffile.imread('electrode.tiff')
s = tau.Solver(img==0)
tau_x = s.solve()

print(f"tau = {s.tau[0]:.4f}, D_eff = {s.D_eff[0]:.4f}*D_0")
```

--------------------------------

### Visualize Simple Test Structures

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Generates and visualizes several predefined test structures (fcc, blocks, 2d_diagonals, 2d_zigzag, 3d_diagonals) using matplotlib for 3D voxel plotting. This is useful for initial inspection of the structures.

```python
from taufactor.utils import create_fcc_cube, \
    create_stacked_blocks, create_2d_diagonals, \
    create_3d_diagonals, create_2d_zigzag
import matplotlib.pyplot as plt
import numpy as np

Nx = 16
fields = {}

fields["fcc"] = create_fcc_cube(Nx)
fields["blocks"] = create_stacked_blocks(Nx, features=1)
fields["2d_diagonals"] = create_2d_diagonals(Nx, features=1)
fields["2d_zigzag"] = create_2d_zigzag(Nx, features=1)
fields["3d_diagonals"] = create_3d_diagonals(Nx, features=1)

fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(15, 4), dpi=100)
for i, structure in enumerate(fields.keys()):
    axes[i].remove()
    axes[i] = fig.add_subplot(1, 5, i + 1, projection='3d')
    axes[i].voxels(np.transpose(fields[structure], (2,1,0))==1, facecolors='gray', edgecolor='black', alpha=1.0)
    axes[i].set_zlim(0, Nx)
    axes[i].set_aspect('equal')
    axes[i].set_title(structure)
    axes[i].set_xlim(0, Nx)
    axes[i].set_ylim(0, Nx)
plt.tight_layout()
plt.show()
```

--------------------------------

### Load and Visualize Electrode Data

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Loads a 3D TIFF image of a battery electrode and displays the first slice. Ensure the 'electrode.tiff' file is in the correct directory. The labels dictionary maps phase names to their corresponding integer values.

```python
electrode = tifffile.imread('electrode.tiff')
print("Stack shape:", electrode.shape)

labels = {"pore":0, "NMC":85, "CBD":170}
boxsize = 100e-6 # physical domain in m
px = boxsize/electrode.shape[0] # pixel resolution in m

plt.imshow(electrode[:,:,0].T, cmap='gray_r', interpolation='none')
plt.show()
```

--------------------------------

### Clone TauFactor Repository

Source: https://github.com/tldr-group/taufactor/blob/main/CONTRIBUTING.md

Clone your forked repository locally to begin development. Replace 'your_name_here' with your GitHub username.

```bash
git clone git@github.com:your_name_here/taufactor.git
```

--------------------------------

### Run Benchmark Study for Solvers

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

This code iterates through different solver types ('PeriodicSolver', 'PeriodicMultiPhaseSolver'), runs a benchmark study with specified parameters (N=64, structure='diagonal2d', features=2), and collects the results. It's useful for comparing solver performance.

```python
all_rows = []
for solver in ['PeriodicSolver', 'PeriodicMultiPhaseSolver']:
    res = run_benchmark_study(
        Ns=[64],
        solver=solver,
        structure='diagonal2d',
        features=2,
        write_file=False,
    )
    all_rows.extend(res)
df = pd.DataFrame(all_rows).sort_values('N').reset_index(drop=True)
df
```

--------------------------------

### Solve with Debug Verbosity for Convergence Analysis

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Solves for tortuosity with 'debug' verbosity enabled, plotting the change in tortuosity for all batch elements. This helps visualize the convergence progress of each microstructure.

```python
s = tau.PeriodicElectrodeSolver(batch, device='cuda', spacing=Lx/nx)
tau_pore = s.solve(iter_limit=50000, conv_crit=1e-3, verbose='debug', plot_interval=1)
```

--------------------------------

### Simulate and Save Impedance (large-small)

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Performs an impedance simulation for the 'large-small' configuration with an increased iteration limit. The 'nyquist' mode is used, and the simulation output is saved.

```python
c = tau.PeriodicImpedanceSolver(fields["large-small"], mode='nyquist')
c.solve(iter_limit=1000000, conv_crit=-1, verbose='plot', plot_interval=100)
pkl_path = save_solver_object(c, "impedance_solve_large-small")
```

--------------------------------

### Plotting and Display Configuration

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Configures plot limits, legends, and grids for visualization. Ensures tight layout and displays the plot.

```python
ax[1,i].set_ylim(0.9, 2.2)
    ax[1,i].legend(loc='upper right')
    ax[2,i].legend(loc='upper right')
    ax[1,i].grid()
    ax[2,i].grid()
    ax[2,i].set_yticks([0,0.5, 1])
ax[2,0].set_ylim(0, 1)
ax[2,1].set_ylim(0.4, 1)


plt.tight_layout()
plt.show()
```

--------------------------------

### Import Libraries for Transmission Line Modelling

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/11-transmission-line-model.ipynb

Imports necessary libraries, NumPy for numerical operations and Matplotlib for plotting.

```python
import numpy as np
import matplotlib.pyplot as plt
```

--------------------------------

### Plotting Field Comparisons and Concentration Profiles

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Generates plots comparing different field structures and concentration profiles, including tau(x) and steady-state concentrations.

```python
structure = 'diag-pillars'
fig, ax = plt.subplots(nrows=3, ncols=2, figsize=(7, 6), dpi=300,
                       constrained_layout=True,
                       gridspec_kw={'height_ratios': [1.0, 0.5, 0.5]})

for i, suffix in enumerate(['_c_classic', '_c']):
    data1 = fields[structure+suffix][:, :150, 0].T
    data2 = fields[structure][:, :150, 0].T
    alpha = np.clip(data2==0, 0, 1).astype(float)
    ax[0, i].imshow(data1, cmap='turbo_r', interpolation='none', vmin=0, vmax=1)
    ax[0, i].imshow(data2, cmap='gray', interpolation='none', alpha=alpha)
    ax[0, i].set_aspect('equal')
    ax[0, i].set_xlim(0, nx)
    ax[0, i].set_ylim(0, 150)
    ax[0, i].set_xticks([])
    ax[0, i].set_yticks([])

dx = 1 / nx
x = np.arange(0, nx)+0.5
ax[1,0].plot(np.arange(1, nx), tau_x_classic[structure], label='$\tau(x)$', color='green', linestyle='-')
ax[1,0].plot(np.arange(0, nx), 1.5*np.ones(nx), label='$\tau_c$', color='limegreen', linestyle='--')
ax[1,1].plot(np.arange(0, nx), tau_x[structure], label='$\tau(x)$', color='green', linestyle='-')
ax[1,1].plot(np.arange(0, nx), 1.5*np.ones(nx), label='$\tau_\text{EIS}$', color='limegreen', linestyle='--')

for i, suffix in enumerate(['_c_classic', '_c']):
    ax[2,i].plot(np.arange(0, 250), c_x[structure+suffix], label='$\bar{c}(x)$', color='blue', linestyle='-')

c_tau = solve_steady_state_c(0.5*(vol_x[structure][:-1]+vol_x[structure][1:])/tau_x_classic[structure])
ax[2,0].plot(x, c_tau, label='$c^*(\tau(x))$', color='cornflowerblue', linestyle='--')
c_tau = solve_steady_state_c(0.5/1.5*np.ones(nx-1))
ax[2,0].plot(x, c_tau, label='$c^*(\tau_c)$', color='black', linestyle='--')
k0 = np.mean(vol_x[structure]) / np.mean(a_x[structure]) #/ nx**2
c_tau = solve_steady_state_reactive_zero_flux(vol_x[structure]/tau_x[structure], k0*a_x[structure])
ax[2,1].plot(np.arange(0, 251), c_tau, label='$c^*(\tau(x))$', color='cornflowerblue', linestyle='--')
c_tau = solve_steady_state_reactive_zero_flux(0.5/1.5*np.ones(nx), k0*a_x[structure])
ax[2,1].plot(np.arange(0, 251), c_tau, label='$c^*(\tau_\text{EIS})$', color='black', linestyle='--')

for i, suffix in enumerate(['_c_classic', '_c']):
    ax[1,i].set_xlim(0, 250)
    ax[2,i].set_xlim(0, 250)
    ax[1, i].set_xticks([])
    ax[2, i].set_xticks([])
```

--------------------------------

### Compare Performance of Specific Surface Area on Segmented Particles

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Measures and prints the wall time for different specific surface area calculation methods when applied to individual segmented particles. This highlights performance differences on a more complex, segmented dataset.

```python
timer = []
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='face_counting', device='cuda')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='gradient', smoothing=False, device='cuda')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='gradient', smoothing=True, device='cuda')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='marching_cubes', smoothing=False, device='cpu')
timer.append(time.time())
surf = tau.specific_surface_area(snow_labels, method='marching_cubes', smoothing=True, device='cpu')
timer.append(time.time())

print(f"Face_counting: {timer[1]-timer[0]:.4f} s\n"
      f"Gradient w/o smoothing: {timer[2]-timer[1]:.4f} s\n"
      f"Gradient w smoothing: {timer[3]-timer[2]:.4f} s\n"
      f"Marching w/o smoothing: {timer[4]-timer[3]:.4f} s\n"
      f"Marching w smoothing: {timer[5]-timer[4]:.4f} s")
```

--------------------------------

### Calculate and Print Volume Fractions

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Computes the volume fraction for the 'gray' phase of each structure using the `tau.volume_fraction` function and prints the results.

```python
for key, struc in structures.items():
    print("Volume fraction of "+key+f": {tau.volume_fraction(struc, labels)['gray']:.4f}")
```

--------------------------------

### Visualize Generated Microstructures

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Visualizes the six generated microstructure patterns using matplotlib. Each pattern is displayed as a 2D grayscale image.

```python
fig, axes = plt.subplots(nrows=1, ncols=6, figsize=(16, 8), dpi=300)
for i, structure in enumerate(structures):
    data2d = fields[structure][:, :, 0].T
    axes[i].imshow(data2d, cmap='gray', interpolation='none')
    axes[i].set_aspect('equal')
    axes[i].set_title(structure)
    axes[i].set_xlim(0, nx)
    axes[i].set_ylim(0, nx)

plt.tight_layout()
plt.show()
```

--------------------------------

### Prepare Microstructure Batch for Solver

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Stacks the generated microstructure fields into a single numpy array for batch processing by the taufactor solver. This prepares the data for efficient computation on the GPU.

```python
batch = np.stack([fields[structure] for structure in structures], axis=0)
batch.shape
```

--------------------------------

### Load Structure and Define Labels

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Loads a 3D structure from a TIFF file and defines labels for different phases. This is a prerequisite for surface area calculations on segmented images.

```python
structure = tifffile.imread('electrode.tiff')
print("Stack shape:", structure.shape)
labels = {"pore":0, "NMC":85, "CBD":170}
```

--------------------------------

### Calculate Volume Fractions and Surface Areas

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Computes the volume fractions and specific surface areas for each phase in the electrode. The `spacing` parameter should be provided as a tuple of physical dimensions per pixel for each axis.

```python
vol_frac = tau.metrics.volume_fraction(electrode , labels)
print("volume fractions: " + str(vol_frac))

surface_areas = tau.metrics.specific_surface_area(electrode, phases=labels, spacing=(px,px,px), method='face_counting')
print("specific surface areas: " + str(surface_areas))
```

--------------------------------

### Plotting Hardware Reference Lines

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Adds reference lines and text labels to a plot to indicate hardware specifications for different GPU models.

```python
for x, y, txt in [(3e4, 4, 'RTX 500 Ada (4 GB)'), (3e5, 48, 'RTX A6000 (48 GB)'), (3e6, 80, 'A100 (80 GB)')]:
    ax2.axhline(y, color='black', linestyle='--', linewidth=1)
    ax2.text(x, y, txt, fontsize=8, color='black')
```

--------------------------------

### Plot Solver Comparison Metrics

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Generates three subplots comparing different solvers: solve time vs. voxel count, RAM usage (torch_max) vs. voxel count, and Tau Factor vs. voxel count. Uses predefined color and marker styles.

```python
colors = {
    'Solver': '#1f77b4',
    'PeriodicSolver': '#2ca02c',
    'AnisotropicSolver': '#ff7f0e',
    'MultiPhaseSolver': '#d62728',
}
markers = {
    'Solver': 'o',
    'PeriodicSolver': 's',
    'AnisotropicSolver': '^',
    'MultiPhaseSolver': 'D',
}

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4), dpi=200)
for label, g in df_compare.groupby('solver_label'):
    g = g.sort_values('N')
    x = g['N'] ** 3
    ax1.loglog(
        x, g['solve_time'], marker=markers.get(label, 'o'), linewidth=2,
        color=colors.get(label), label=label
    )
    ax2.loglog(
        x, g['torch_max']/1024, marker=markers.get(label, 'o'), linewidth=2,
        color=colors.get(label), label=label
    )
    ax3.semilogx(
        x, g['taufactor'], marker=markers.get(label, 'o'), linewidth=2,
        color=colors.get(label), label=label, linestyle='--'
    )
```

--------------------------------

### Run Benchmark Study for Solver

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Runs a benchmark study for a specific solver ('Solver') across various voxel counts (Ns) and a fixed convergence criterion. Results are stored in a pandas DataFrame.

```python
from taufactor.benchmark import run_benchmark_study
import pandas as pd

Ns = [32, 64, 100, 128, 200, 256, 300]
res = run_benchmark_study(
    Ns=Ns,
    devices=['cuda'],
    conv_crit_values=[1e-3],
    structure='fcc',
    solver='Solver',
    write_file=False,
)

df_solver = pd.DataFrame(res).sort_values('N').reset_index(drop=True)
df_solver
```

--------------------------------

### Solve Tortuosity for Battery Electrode

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Initializes and solves for tortuosity of a battery electrode using CUDA acceleration. Suitable for detailed electrode simulations.

```python
s = tau.ElectrodeSolver(electrode, device='cuda', spacing=px,
                        conductive_label=labels["connected_pore"], reactive_label=labels["NMC"])
tau_electrode = s.solve(iter_limit=50000, conv_crit=1e-4, plot_interval=10, verbose='plot')
```

--------------------------------

### Calculate Specific Surface Area and Interfacial Areas for Binary Structure

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Compares the specific surface area of a binary structure with its normalized interfacial area. Useful for validating calculations on simple structures.

```python
print(tau.specific_surface_area(sphere, method='face_counting'))
print(tau.interfacial_areas(sphere, normalize=True))
```

--------------------------------

### Calculate Effective Diffusivity and Tortuosity

Source: https://github.com/tldr-group/taufactor/blob/main/README.md

Load image data, create a Solver object, and call the solve function to compute effective diffusivity and tortuosity. Ensure image data consists of 1s for the conductive phase and 0s otherwise.

```python
import taufactor as tau
import tifffile

# load image
img = tifffile.imread('path/filename')
# ensure 1s for conductive phase and 0s otherwise.

# create a solver object with loaded image
s = tau.Solver(img)

# call solve function
s.solve()

# view effective diffusivity and tau
print(s.D_eff, s.tau)

```

--------------------------------

### Import Libraries for Surface Area Computation

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Imports necessary libraries including matplotlib, numpy, taufactor.metrics, tifffile, and scipy.ndimage for data visualization, numerical operations, surface area calculations, image handling, and image rotation.

```python
import matplotlib.pyplot as plt
import numpy as np
import taufactor.metrics as tau
import tifffile
from scipy.ndimage import rotate
```

--------------------------------

### Import Libraries for Tortuosity Analysis

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Imports necessary libraries including numpy for numerical operations, tifffile for image handling, matplotlib for plotting, and taufactor for the core calculations.

```python
import numpy as np
import tifffile
import matplotlib.pyplot as plt
import taufactor as tau
```

--------------------------------

### Compute Impedance using TLM Formula

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/10-tau-x-validation.ipynb

Computes impedance for both heterogeneous and homogenized cases using the recursive TLM formula. Requires pre-computed R, C, and frequency data.

```python
from taufactor.utils import compute_impedance

frequency_tau_x = {}
Z_TLM_tau_x = {}
Z_TLM_tau_mean = {}
for i, structure in enumerate(structures):
    # Compute impedance for heterogeneous distribution
    R = tau_x[structure] / vol_x[structure]
    R[vol_x[structure] == 0] = 1e30
    R[np.isnan(tau_x[structure])] = 1e30
    C = a_x[structure]

    freq_0 = np.mean(vol_x[structure]) / np.mean(C) / nx**2
    freq_range = freq_0 * 2 ** np.arange(-3, 10, 0.5)
    freq = freq_range[::-1].copy()
    frequency_tau_x[structure] = freq
    Z_TLM_tau_x[structure] = compute_impedance(R, C, freq)

    # Compute impedance for homogenized case
    R_mean = np.full(nx, 1/np.mean(vol_x[structure]))
    C_mean = np.full(nx, np.mean(C))
    Z_TLM_tau_mean[structure] = compute_impedance(R_mean, C_mean, freq)
```

--------------------------------

### Benchmark Multiple Solvers

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/12-solver-benchmark.ipynb

Iterates through a list of solver specifications, running a benchmark study for each and collecting results into a single DataFrame. Results are written to 'fcc_benchmark_results.txt'.

```python
solver_specs = [
    {'solver': 'Solver', 'label': 'Solver', 'solver_kwargs': {}},
    {'solver': 'PeriodicSolver', 'label': 'PeriodicSolver', 'solver_kwargs': {}},
    {'solver': 'AnisotropicSolver', 'label': 'AnisotropicSolver', 'solver_kwargs': {'spacing': (1.0, 1.0, 1.0)}},
    {'solver': 'MultiPhaseSolver', 'label': 'MultiPhaseSolver', 'solver_kwargs': {}},
]

all_rows = []
for spec in solver_specs:
    print(f"Running {spec['label']}...")
    rows = run_benchmark_study(
        Ns=Ns,
        devices=['cuda'],
        conv_crit_values=[1e-3],
        structure='fcc',
        solver=spec['solver'],
        solver_kwargs=spec['solver_kwargs'],
        outfile='fcc_benchmark_results.txt',
        write_file=True,
    )
    for r in rows:
        r['solver_label'] = spec['label']
    all_rows.extend(rows)

df_compare = pd.DataFrame(all_rows)
# df_compare
```

--------------------------------

### Visualize Downsampled Tortuosity Data

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/01-tortuosity.ipynb

Generates a plot comparing the original tortuosity, volume fraction, and interfacial area profiles with downsampled versions. This visualization helps in understanding the effect of downsampling on the data and assessing the smoothness of the profiles.

```python
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(8, 3), dpi=100)
x_orig = np.linspace(0, 1, s.vol_x[0].shape[0])
axes[0].plot(x_orig, s.vol_x[0], 'r--', label=f'Original ({x_orig.size} points)')
axes[0].set_ylim(0.42, 0.48)
axes[1].plot(x_orig, s.a_x[0], 'r--')
axes[1].set_ylim(1e5, 1.7e5)
axes[2].plot(x_orig, s.tau_x[0], 'r--')
axes[2].set_ylim(0, 4)

downsample_and_plot(s, factor=2, color='blue')
downsample_and_plot(s, factor=4, color='green')

titles = ['Volume fraction $\varepsilon(x)$', 'Interfacial area $a(x)$', 'Tortuosity $\tau(x)$']
for ax, title in zip(axes, titles):
    ax.set_title(title)
    ax.set_xlim(0, 1)
    ax.grid(True, alpha=0.25)
fig.legend(loc='upper center', ncol=3, frameon=False, bbox_to_anchor=(0.5, 0))
plt.tight_layout()
plt.show()
```

--------------------------------

### Visualize Geometric Structures

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Renders the created 3D structures using matplotlib's 3D voxel plotting. This helps visualize the shapes before surface area calculations.

```python
fig, axes = plt.subplots(nrows=1, ncols=5, figsize=(16, 4), dpi=100)

i = 0
for key, value in structures.items():
    axes[i].remove()
    axes[i] = fig.add_subplot(1, 5, i + 1, projection='3d')

    axes[i].voxels(value==1, facecolors='gray', edgecolor='black', alpha=1.0)
    axes[i].voxels(value==2, facecolors='blue', edgecolor='black', alpha=1.0)
    axes[i].voxels(value==3, facecolors='green', edgecolor='black', alpha=1.0)
    axes[i].set_aspect('equal')
    axes[i].set_title(key)
    axes[i].set_xlim(0, nx)
    axes[i].set_ylim(0, ny)
    axes[i].set_zlim(0, nz)
    i = i + 1
plt.tight_layout()
plt.show()
```

--------------------------------

### Calculate Specific Surface Areas

Source: https://github.com/tldr-group/taufactor/blob/main/docs/notebooks/02-surface-areas.ipynb

Calculates the specific surface area for the 'gray' phase of each structure using 'face_counting', 'gradient', and 'marching_cubes' methods. The 'gradient' and 'marching_cubes' methods utilize smoothing and specify the device.

```python
a_faces = {}
a_marching = {}
a_gradient = {}

for key, struc in structures.items():
    a_faces[key] = tau.specific_surface_area(struc,
                                             method='face_counting',
                                             phases=labels)['gray']
    a_gradient[key] = tau.specific_surface_area(struc,
                                                method='gradient',
                                                phases=labels,
                                                smoothing=False)['gray']
    a_marching[key] = tau.specific_surface_area(struc,
                                                method='marching_cubes',
                                                phases=labels,
                                                smoothing=False,
                                                device='cpu')['gray']

print("Spec. surf. area:     faces         \t|      gradient      \t|   marching")
for key, struc in structures.items():
    print(key+f":   \t{a_faces[key]:.5f} ({((a_faces[key]-a_theo[key])/a_theo[key]*100):.2f} %)\t|"
        f"  {a_gradient[key]:.5f} ({(a_gradient[key]-a_theo[key])/a_theo[key]*100:.2f} %)\t|"
        f"  {a_marching[key]:.5f} ({(a_marching[key]-a_theo[key])/a_theo[key]*100:.2f} %)")
```