OxyGenie.diffusion

Module contents

OxyGenie.diffusion.FromCustom	A class to initialize concentration, diffusion coefficients, and reaction rates from an custom distributions.
OxyGenie.diffusion.FromIMG	A class to initialize concentration, diffusion coefficients, and reaction rates from an image.
OxyGenie.diffusion.FromPGVNet	A class to initialize concentration, diffusion coefficients, and reaction rates from a vascular network generated by the pgvnet sub-module.
OxyGenie.diffusion.Simu_plot
OxyGenie.diffusion.SimulationParams	A class to define and manage the parameters for a 2D simulation.
OxyGenie.diffusion.run_simulation	Runs the simulation of the 2D diffusion equation:
OxyGenie.diffusion.compute_step_rk4	Performs one step of the 4th-order Runge-Kutta (RK4) method to update the concentration field.
OxyGenie.diffusion.compute_step_euler	Performs one step of the Euler method to update the concentration field.
OxyGenie.diffusion.crit_stab	Evaluates the Courant–Friedrichs–Lewy (CFL) condition for numerical stability with the maximum values of alphax and alphay.
OxyGenie.diffusion.F	Computes the rate of change of the concentration field C based on the diffusion and reaction terms.

class OxyGenie.diffusion.FromCustom(C_init, D_init, k_init)

A class to initialize concentration, diffusion coefficients, and reaction rates from an custom distributions.

Warning

• C_init and D_init values must be normalized to [0, 1].

• k_init values must be normalized to [-1, 1].

Parameters:

• C_init (2D np.ndarray) – Initial concentration distribution scaled in the range [0, 1].

• D_init (2D np.ndarray) – Initial diffusion coefficient distribution scaled in the range [0, 1].

• k_init (2D np.ndarray) – Initial reaction rate distribution scaled in the range [-1, 1].

Returns:

• C2D np.ndarray

  Initialized concentration field based on grid values.

• D2D np.ndarray

  Diffusion coefficient field initialized based on grid values.

• k2D np.ndarray

  Absorbtion rate field initialized based on grid values.

• alphax2D np.ndarray

  Prefactor for diffusion in the x-direction.

• alphay2D np.ndarray

  Prefactor for diffusion in the y-direction.

• betax2D np.ndarray

  Prefactor for the diffusion flux term in the x-direction.

• betay2D np.ndarray

  Prefactor for the diffusion flux term in the y-direction

Return type:

tuple

from OxyGenie.diffusion import *

params = {
    "D": 1e-5, "k": 8, "Lx": 0.05, "Ly": 0.05, "T": 0.5,
    "nt": 2500, "nx": 256 * 2, "ny": 256 * 2,
    "initial_concentration": 100.0, "step": 10,
}
simparams = SimulationParams(**params)
C = np.zeros((simparams.nx,simparams.ny))
C[120:130, 120:130] = 1

X, Y = np.meshgrid(np.linspace(-1, 1, simparams.ny), np.linspace(-1, 1, simparams.nx))
D = np.exp(-((X)**2+(Y-1)**2)/1.5)

k = np.ones((simparams.nx, simparams.ny))

from_Custom = FromCustom(C, D, k)
C, D, k, alphax, alphay, betax, betay = from_Custom(params)

class OxyGenie.diffusion.FromIMG(vasc_path)

A class to initialize concentration, diffusion coefficients, and reaction rates from an image.

Note

• The provided image is resized to match the grid size defined by params.nx and params.ny.

• The concentration (C), diffusion coefficient (D), and reaction rate (k) fields are initialized based on chanel values in the grid.

• Return via __call__

Parameters:

vasc_path (str) – Path to the image file containing spatial information. The image is used to define regions with different diffusion and reaction properties based on the channel pixel values.

Returns:

• C2D np.ndarray

  Initialized concentration field based on grid values.

• D2D np.ndarray

  Diffusion coefficient field initialized based on grid values.

• k2D np.ndarray

  Absorbtion rate field initialized based on grid values.

• alphax2D np.ndarray

  Prefactor for diffusion in the x-direction.

• alphay2D np.ndarray

  Prefactor for diffusion in the y-direction.

• betax2D np.ndarray

  Prefactor for the diffusion flux term in the x-direction.

• betay2D np.ndarray

  Prefactor for the diffusion flux term in the y-direction

Return type:

tuple

from OxyGenie.diffusion import *

params = {
    "D": 1e-5, "k": 8, "Lx": 0.05, "Ly": 0.05, "T": 0.5,
    "nt": 2500, "nx": 256 * 2, "ny": 256 * 2,
    "initial_concentration": 100.0, "step": 10,
}
simparams = SimulationParams(**params)
from_IMG = FromIMG("/path.jpg")
C, D, k, alphax, alphay, betax, betay = from_IMG(params)

class OxyGenie.diffusion.FromPGVNet(grid)

A class to initialize concentration, diffusion coefficients, and reaction rates from a vascular network generated by the pgvnet sub-module.

Note

• The provided grid is resized to match the grid size defined by params.nx and params.ny.

• The concentration (C), diffusion coefficient (D), and reaction rate (k) fields are initialized based on values in the grid.

• Diffusion coefficient and absorbtion rate are adjusted based on pixel values, where:

  • Areas with higher grid values represent regions with higher diffusion coefficients.

  • Regions with pixel values below a threshold represent areas with production or absorption rates.

• Return via __call__

Parameters:

grid (2D np.ndarray) – Vascular network generated by the pgvnet sub-module. The grid values are used to define regions with different diffusion and absorbtion properties based on thresholds.

Returns:

• C2D np.ndarray

  Initialized concentration field based on grid values.

• D2D np.ndarray

  Diffusion coefficient field initialized based on grid values.

• k2D np.ndarray

  Absorbtion rate field initialized based on grid values.

• alphax2D np.ndarray

  Prefactor for diffusion in the x-direction.

• alphay2D np.ndarray

  Prefactor for diffusion in the y-direction.

• betax2D np.ndarray

  Prefactor for the diffusion flux term in the x-direction.

• betay2D np.ndarray

  Prefactor for the diffusion flux term in the y-direction

Return type:

tuple

from OxyGenie.diffusion import *
import OxyGenie.pgvnet as pgv

grid = pgv.simple_generation(grid_size=1280)[0]
params = {
    "D": 1e-5, "k": 8, "Lx": 0.05, "Ly": 0.05, "T": 0.5,
    "nt": 2500, "nx": 256 * 2, "ny": 256 * 2,
    "initial_concentration": 100.0, "step": 10,
}
simparams = SimulationParams(**params)
from_pgvnet = FromPGVNet(grid)
C, D, k, alphax, alphay, betax, betay = from_pgvnet(params)

class OxyGenie.diffusion.Simu_plot

classmethod anim(params, L_result, anim=True)

Creates an animated visualization showing the evolution of the concentration field, along with graphs of related metrics such as total concentration over time.

classmethod anim_vect(params, L_result, s_ech=5)

Produces an animated vector field representation of concentration gradients over time, combined with concentration maps.

classmethod champ_vect(params, L_result, i, s_ech=5)

Visualizes the concentration gradient as a vector field overlay on a concentration map for a specific time step.

classmethod simple(params, C_result, D=None)

Generates a static plot showing the final concentration field. Optionally includes a plot of the diffusion coefficient if D is provided.

class OxyGenie.diffusion.SimulationParams(**kwargs)

A class to define and manage the parameters for a 2D simulation.

This class handles both user-defined and derived parameters needed for numerical simulations, including spatial and temporal discretizations. It ensures consistency among related parameters (e.g., grid size and spacing) and recalculates derived parameters when key attributes are updated.

Parameters:

**kwargs (dict) –

Arbitrary keyword arguments to initialize the simulation parameters:

• Dfloat, optional

  Diffusion coefficient (default: 1e-5).

• kfloat, optional

  Absorbtion rate constant (default: 0).

• Lxfloat, optional

  Length of the simulation domain in the x-direction (default: 0.01).

• Lyfloat, optional

  Length of the simulation domain in the y-direction (default: 0.01).

• Tfloat, optional

  Total simulation time (default: 0.01).

• initial_concentrationfloat, optional

  Initial concentration (default: 100.0).

• speedfloat, optional

  Speed of the simulation (default: 1).

• stepint, optional

  Simulation step parameter (default: 1).

• ntint, optional

  Number of time steps (default: 10000).

• nxint, optional

  Number of spatial steps in the x-direction (default: 100).

• nyint, optional

  Number of spatial steps in the y-direction (default: 100).

• dtfloat, optional

  Time step size (computed if not provided).

• dxfloat, optional

  Grid spacing in the x-direction (computed if not provided).

• dyfloat, optional

  Grid spacing in the y-direction (computed if not provided).

params = {
    "D": 1e-5, "k": 8, "Lx": 0.05, "Ly": 0.05, "T": 0.5, "nt": 2500, "nx": 256 * 2, "ny": 256 * 2,
    "initial_concentration": 100.0, "speed": 1, "step": 10,
}
simparams = SimulationParams(**params)
print(simparams)

property Lx

property Ly

property T

property dt

property dx

property dy

property nt

property nx

property ny

OxyGenie.diffusion.F(C, k, cax, cay, cbx, cby)

Computes the rate of change of the concentration field C based on the diffusion and reaction terms.

This function calculates the right-hand side of the diffusion-reaction equation using finite differences for the diffusion terms and a simple linear reaction term.

Warning

This function is only working for RK4 compute step scheme

OxyGenie.diffusion.N_tot(C, dx, dy)

Calculates the total number of particles in a 2D concentration field.

Warning

-Avogadro’s number is used to convert molar concentration to the number of particles.

-The boundary of width l = 1 is excluded from the sum to avoid edge effects from the boundary condintions.

Parameters:

• C (2D np.ndarray) – The 2D array representing the concentration field (in moles per square meter). The boundary region is excluded from the calculation.

• dx (float) – The grid spacing in the x-direction (in meters).

• dy (float) – The grid spacing in the y-direction (in meters).

Returns:

The total number of particles within the central region of the grid. The result is dimensionless.

Return type:

float

OxyGenie.diffusion.compute_step_euler(C, D, k, cax, cay, cbx, cby, dt)

Performs one step of the Euler method to update the concentration field.

The update is based on the following form of the diffusion equation:

\(\frac{\partial C(\textbf{r}, t)}{\partial t} = \text{div}(D(\textbf{r}) \cdot \overrightarrow{\text{grad}}(C(\textbf{r}, t)) - k(\textbf{r}) \cdot C(\textbf{r}, t)\)

where the diffusion terms are discretized using finite differences (forward for flux terms and central for diffusion terms).

Note

The boundary conditions are set to be “no flux” (Neumann).

Warning

The stability criterion needs to be respected.

Parameters:

• C (2D np.ndarray) – The current concentration field.

• D (2D np.ndarray) – The diffusion coefficient field, with the same shape as C.

• k (2D np.ndarray) – The absorbtion rate field, with the same shape as C.

• cax (2D np.ndarray) – The precomputed prefactor for the diffusion term in the x-direction.

• cay (2D np.ndarray) – The precomputed prefactor for the diffusion term in the y-direction.

• cbx (2D np.ndarray) – The precomputed prefactor for the diffusion flux term in the x-direction.

• cby (2D np.ndarray) – The precomputed prefactor for the diffusion flux term in the y-direction.

• dt (float) – The time step (dt) used for the RK4 method to update the concentration.

Returns:

The updated concentration field C_new after one RK4 step, with the same shape as C.

Return type:

2D np.ndarray

OxyGenie.diffusion.compute_step_rk4(C, D, k, cax, cay, cbx, cby, dt)

Performs one step of the 4th-order Runge-Kutta (RK4) method to update the concentration field.

The RK4 method is a numerical technique for solving ordinary differential equations (ODEs). This function applies it to the concentration field C in a 2D grid, accounting for diffusion and absorbtion.

Note

The boundary conditions are set to be “no flux” (Neumann).

Parameters:

• C (2D np.ndarray) – The current concentration field.

• D (2D np.ndarray) – The diffusion coefficient field, with the same shape as C.

• k (2D np.ndarray) – The absorbtion rate field, with the same shape as C.

• cax (2D np.ndarray) – The precomputed prefactor for the diffusion term in the x-direction.

• cay (2D np.ndarray) – The precomputed prefactor for the diffusion term in the y-direction.

• cbx (2D np.ndarray) – The precomputed prefactor for the diffusion flux term in the x-direction.

• cby (2D np.ndarray) – The precomputed prefactor for the diffusion flux term in the y-direction.

• dt (float) – The time step (dt) used for the RK4 method to update the concentration.

Returns:

The updated concentration field C_new after one RK4 step, with the same shape as C.

Return type:

2D np.ndarray

OxyGenie.diffusion.crit_stab(alphax, alphay)

Evaluates the Courant–Friedrichs–Lewy (CFL) condition for numerical stability with the maximum values of alphax and alphay.

The function checks if the sum of the maximum values of alphax and alphay exceeds the stability threshold (0.5). If the criterion is not met, a ValueError is raised.

Parameters:

• alphax (2D np.ndarray) – Here: \(\alpha_x = \frac{D \cdot \Delta t}{(\Delta x)^2}\)

• alphay (2D np.ndarray) – Here: \(\alpha_y = \frac{D \cdot \Delta t}{(\Delta y)^2}\)

Raises:

ValueError – If the stability criterion is not respected (stab_coef >= 0.5). The error message includes the computed stability coefficient.

Returns:

Prints a message indicating whether the stability criterion is respected.

Return type:

None.

OxyGenie.diffusion.run_simulation(params, init_method, C_0_cst=True, save_last_only=False, C0=None)

Runs the simulation of the 2D diffusion equation:

\[\frac{\partial C(\mathbf{r}, t)}{\partial t} = \nabla \cdot \big(D(\mathbf{r}) \cdot \nabla C(\mathbf{r}, t)\big) - k(\mathbf{r}) \cdot C(\mathbf{r}, t)\]

Parameters:

• params (object) – Contains the simulation parameters, including time step (dt), grid resolution (dx, dy), and number of time steps (nt).

• init_method (callable) – A function that initializes the state variables, returning C, D, k, alphax, alphay, betax, and betay.

• C_0_cst (bool, optional) – If True, maintains the initial constant value of C throughout the simulation (default: True).

• save_last_only (bool, optional) – If True, only saves the final state of the simulation (default: False).

• C0 (2D np.ndarray, optional) – An optional initial concentration array to override the default initialization.

Returns:

A list containing the concentration fields at each saved time step. If save_last_only is True, only the final state is saved.

Return type:

list of 2D np.ndarray

Raises:

ValueError – If the stability criterion is not respected.

Notes

The function uses an explicit Euler method to solve the diffusion equation. The stability criterion is checked before the simulation begins using the crit_stab function. The total amount of substance in the system is printed before and after the simulation.

Examples

To run a simulation with default parameters:

from OxyGenie.diffusion import *

params = {
    "D": 1e-5, "k": 8, "Lx": 0.05, "Ly": 0.05, "T": 0.5,
    "nt": 2500, "nx": 256 * 2, "ny": 256 * 2,
    "initial_concentration": 100.0, "step": 10,
}
simparams = SimulationParams(**params)

C = np.zeros((simparams.nx,simparams.ny))
C[120:130, 120:130] = 1
X, Y = np.meshgrid(np.linspace(-1, 1, simparams.ny), np.linspace(-1, 1, simparams.nx))
D = np.exp(-((X)**2+(Y-1)**2)/1.5)
k = np.ones((simparams.nx, simparams.ny))

results = run_simulation(simparams, FromCustom(C, D, k), C_0_cst=True)

Simu_plot.simple(simparams, results[-1])