[Python] How to best solve partial differential equations in Python accurately but quickly

I am trying to reproduce 4.1.3 Emission from Solid Materials as per this PDF.

https://www.rivm.nl/bibliotheek/rapporten/2017-0197.pdf

Produced this version of the model and this is the best I have gotten.

# Import necessary libraries
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# -------------------------------
# Constants
# -------------------------------
R_J_per_mol_K = 8.314 # Universal gas constant [J/(mol·K)]

# -------------------------------
# Scenario Selection
# -------------------------------
selected_scenario = 1 # Change this to select scenarios (1, 2, or 3)

# -------------------------------
# Input Parameters
# -------------------------------
# General Parameters (Common across all scenarios)
body_weight_kg = 60 # Body weight [kg]
frequency_per_year = 1 # Frequency of exposure [per year]

# Scenario Specific Parameters
if selected_scenario == 1:
# Scenario 1 (Base Scenario)
product_surface_area_cm2 = 10 # Product surface area [cm²]
product_thickness_mm = 5 # Product thickness [mm]
product_density_g_per_cm3 = 0.1 # Product density [g/cm³]
diffusion_coefficient_m2_per_s = 1E-08 # Diffusion coefficient [m²/s]
weight_fraction_substance = 0.6 # Weight fraction of the substance [fraction]
product_air_partition_coefficient = 0.5 # Product/air partition coefficient [-]

room_volume_m3 = 0.5 # Room volume [m³]
ventilation_rate_per_hour = 58 # Ventilation rate [per hour]
inhalation_rate_m3_per_hour = 10 # Inhalation rate [m³/hr]

start_exposure_h = 1 # Start exposure time [hours]
exposure_duration_h = 2 # Exposure duration [hours]
mass_transfer_coefficient_m_per_h = 10 # Mass transfer coefficient [m/h]

#Absorption
absorption_fraction = 1 # Absorption fraction [fraction]

# -------------------------------
# Unit Conversions
# -------------------------------
product_surface_area_m2 = product_surface_area_cm2 / 10000 # Convert cm² to m²
product_thickness_m = product_thickness_mm / 1000 # Convert mm to m
product_density_kg_per_m3 = product_density_g_per_cm3 * 1000 # Convert g/cm³ to kg/m³
mass_transfer_coefficient_m_per_s = mass_transfer_coefficient_m_per_h / 3600 # Convert m/h to m/s
ventilation_rate_per_s = ventilation_rate_per_hour / 3600 # Convert per hour to per second
start_exposure_s = start_exposure_h * 3600 # Convert hours to seconds
exposure_duration_s = exposure_duration_h * 3600 # Convert hours to seconds
total_emission_time_s = start_exposure_s + exposure_duration_s # Total time to simulate
frequency_per_day = frequency_per_year / 365 # Frequency per day
inhalation_rate_m3_per_s = inhalation_rate_m3_per_hour / 3600 # Convert inhalation rate from m³/h to m³/s

# -------------------------------
# Initial Concentrations
# -------------------------------
volume_material_m3 = product_surface_area_m2 * product_thickness_m # Volume of material [m³]
mass_material_kg = volume_material_m3 * product_density_kg_per_m3 # Mass of material [kg]
mass_substance_initial_kg = mass_material_kg * weight_fraction_substance # Initial mass of substance [kg]
concentration_initial_kg_per_m3 = mass_substance_initial_kg / volume_material_m3 # Initial concentration [kg/m³]
N = 1000 # Reduced number of layers for faster computation
dx = product_thickness_m / (N - 1) # Spatial step size [m]
t_eval = np.linspace(0, total_emission_time_s, 2000) # Time points for evaluation

# -------------------------------
# Emission Model Definition
# -------------------------------
# Precompute constants
D = diffusion_coefficient_m2_per_s
hm = mass_transfer_coefficient_m_per_s
K = product_air_partition_coefficient
S = product_surface_area_m2
V = room_volume_m3
q = ventilation_rate_per_s

def emission_model(t, y):
C = y[:-1] # Concentration in the material layers [kg/m³]
y_air = y[-1] # Air concentration [kg/m³]

dCdt = np.zeros_like(C)

# No-flux boundary condition at x=0 (back surface)
dCdt[0] = D * (C[1] - C[0]) / dx**2 # Diffusion at back surface [kg/(m²·s²)]

# Diffusion in the material
for i in range(1, N - 1):
dCdt = D * (C[i + 1] - 2 * C + C[i - 1]) / dx**2 # Diffusion in layers [kg/(m²·s²)]

# Flux boundary condition at x=L (front surface)
J_diff = -D * (C[N - 1] - C[N - 2]) / dx # Diffusion flux at surface [kg/(m²·s)]
J_air = hm * (C[N - 1] / K - y_air) # Mass transfer to air [kg/(m²·s)]

dCdt[N - 1] = D * (C[N - 2] - C[N - 1]) / dx**2 - J_air / dx # Change in surface layer concentration [kg/m³/s]

dy_air_dt = (S / V) * hm * (C[N - 1] / K - y_air) - q * y_air # Change in air concentration [kg/m³/s]

return np.concatenate((dCdt, [dy_air_dt]))

# -------------------------------
# Initial Conditions
# -------------------------------
C0 = np.full(N, concentration_initial_kg_per_m3) # Initial concentration in material layers [kg/m³]
y_air_0 = 0.0 # Initial air concentration [kg/m³]
y0 = np.concatenate((C0, [y_air_0])) # Initial condition vector

# -------------------------------
# Solve the ODE system with BDF solver for speed and accuracy
# -------------------------------
solution = solve_ivp(emission_model, [0, total_emission_time_s], y0, t_eval=t_eval, method='BDF', vectorized=False)

# Extract results
C_material = solution.y[:-1, :] # Material concentrations over time [kg/m³]
y_air = solution.y[-1, :] # Air concentrations over time [kg/m³]
time = solution.t # Time vector

# -------------------------------
# Extract Air Concentration during Exposure
# -------------------------------
exposure_indices = np.where((time >= start_exposure_s) & (time <= start_exposure_s + exposure_duration_s))[0]
y_air_exposure = y_air[exposure_indices] # Air concentrations during exposure period [kg/m³]
time_exposure = time[exposure_indices] # Time during exposure period

# Check for valid exposure period
if len(exposure_indices) == 0:
print("No exposure period found. Check start and duration times.")
exit()

# -------------------------------
# Results Formatting
# -------------------------------
def format_value(value):
"""Format a number in scientific notation with one decimal place."""
return f"{value:.1e}"

# -------------------------------
# Display Results
# -------------------------------
C_air_mg_per_m3 = y_air_exposure * 1e6 # Convert air concentration to mg/m³
C_mean_event_mg_per_m3 = np.trapz(C_air_mg_per_m3, time_exposure) / (time_exposure[-1] - time_exposure[0]) # Mean event concentration [mg/m³]
C_mean_day_mg_per_m3 = (C_mean_event_mg_per_m3 * exposure_duration_s) / (24 * 3600) # Mean concentration on day of exposure [mg/m³]
C_year_avg_mg_per_m3 = C_mean_day_mg_per_m3 * frequency_per_day # Year average concentration [mg/m³]
Inhalation_volume_event_m3 = inhalation_rate_m3_per_s * exposure_duration_s # Inhalation volume during exposure [m³]
Dose_external_event_mg_per_kg_bw = (C_mean_event_mg_per_m3 * Inhalation_volume_event_m3) / body_weight_kg # External event dose [mg/kg bw]
Dose_external_day_mg_per_kg_bw = Dose_external_event_mg_per_kg_bw # External dose on day of exposure [mg/kg bw]
Dose_internal_event_mg_per_kg_bw = Dose_external_event_mg_per_kg_bw * absorption_fraction # Internal event dose [mg/kg bw]
Dose_internal_day_mg_per_kg_bw = Dose_internal_event_mg_per_kg_bw # Internal dose on day of exposure [mg/kg bw]
Dose_internal_year_avg_mg_per_kg_bw_per_day = Dose_internal_day_mg_per_kg_bw * frequency_per_day # Internal year average dose [mg/kg bw/day]

# Display Results
print(f"Mean Event Concentration: {format_value(C_mean_event_mg_per_m3)} mg/m³")
print(f"Mean Concentration on Day of Exposure: {format_value(C_mean_day_mg_per_m3)} mg/m³")
print(f"Year Average Concentration: {format_value(C_year_avg_mg_per_m3)} mg/m³")
print(f"External Event Dose: {format_value(Dose_external_event_mg_per_kg_bw)} mg/kg bw")
print(f"External Dose on Day of Exposure: {format_value(Dose_external_day_mg_per_kg_bw)} mg/kg bw")
print(f"Internal Event Dose: {format_value(Dose_internal_event_mg_per_kg_bw)} mg/kg bw")
print(f"Internal Dose on Day of Exposure: {format_value(Dose_internal_day_mg_per_kg_bw)} mg/kg bw/day")
print(f"Internal Year Average Dose: {format_value(Dose_internal_year_avg_mg_per_kg_bw_per_day)} mg/kg bw/day")

# -------------------------------
# Additional Calculations for Plots
# -------------------------------
total_duration_s = start_exposure_s + exposure_duration_s # Total duration
t_s_plot = np.linspace(0, total_duration_s, num=1000) # Time vector for plotting
C_t_kg_per_m3 = np.interp(t_s_plot, time, y_air) # Interpolated air concentration over time [kg/m³]
C_t_mg_per_m3 = C_t_kg_per_m3 * 1e6 # Convert air concentration to mg/m³
delta_t_s = t_s_plot[1] - t_s_plot[0] # Time step
delta_inhaled_volume_m3 = inhalation_rate_m3_per_s * delta_t_s # Inhaled volume per time step [m³]
exposure_mask = t_s_plot >= start_exposure_s # Mask for exposure period
incremental_external_dose_mg_per_kg_bw = np.zeros_like(t_s_plot) # Initialize dose array
incremental_external_dose_mg_per_kg_bw[exposure_mask] = (C_t_mg_per_m3[exposure_mask] * delta_inhaled_volume_m3) / body_weight_kg # Incremental dose [mg/kg bw]
cumulative_external_dose_mg_per_kg_bw = np.cumsum(incremental_external_dose_mg_per_kg_bw) # Cumulative external dose [mg/kg bw]
cumulative_internal_dose_mg_per_kg_bw = cumulative_external_dose_mg_per_kg_bw * absorption_fraction # Cumulative internal dose [mg/kg bw]

# -------------------------------
# Visualization
# -------------------------------
fig, axs = plt.subplots(1, 3, figsize=(24, 6)) # Create a figure with three subplots

# 1. Air Concentration Over Time
axs[0].plot(t_s_plot / 3600, C_t_mg_per_m3, label='Air Concentration (mg/m³)', color='blue')
axs[0].set_title('Air Concentration Over Time')
axs[0].set_xlabel('Time (hours)')
axs[0].set_ylabel('Concentration (mg/m³)')
axs[0].grid(True)
axs[0].legend()
axs[0].axvline(x=start_exposure_s / 3600, color='red', linestyle='--', label='Exposure Start')
axs[0].axvline(x=(start_exposure_s + 15 * 60) / 3600, color='green', linestyle='--', label='15 min after Exposure Start')
axs[0].legend()

# 2. Cumulative External Dose Over Time
axs[1].plot(t_s_plot / 3600, cumulative_external_dose_mg_per_kg_bw, label='Cumulative External Dose (mg/kg bw)', color='orange')
axs[1].set_title('External Dose Over Time')
axs[1].set_xlabel('Time (hours)')
axs[1].set_ylabel('Dose (mg/kg bw)')
axs[1].grid(True)
axs[1].legend()
axs[1].axvline(x=start_exposure_s / 3600, color='red', linestyle='--', label='Exposure Start')

# 3. Cumulative Internal Dose Over Time (Adjusted for absorption)
axs[2].plot(t_s_plot / 3600, cumulative_internal_dose_mg_per_kg_bw, label='Cumulative Internal Dose (mg/kg bw)', color='green')
axs[2].set_title('Internal Dose Over Time')
axs[2].set_xlabel('Time (hours)')
axs[2].set_ylabel('Dose (mg/kg bw)')
axs[2].grid(True)
axs[2].legend()
axs[2].axvline(x=start_exposure_s / 3600, color='red', linestyle='--', label='Exposure Start')

plt.tight_layout() # Adjust layout for clarity
plt.show() # Show the plots

# -------------------------------
# End of Script
# -------------------------------


RESULTS:

Mean Event Concentration: 1.3e-01 mg/m³
Mean Concentration on Day of Exposure: 1.1e-02 mg/m³
Year Average Concentration: 2.9e-05 mg/m³
External Event Dose: 4.3e-02 mg/kg bw
External Dose on Day of Exposure: 4.3e-02 mg/kg bw
Internal Event Dose: 4.3e-02 mg/kg bw
Internal Dose on Day of Exposure: 4.3e-02 mg/kg bw/day
Internal Year Average Dose: 1.2e-04 mg/kg bw/day




EXPECTED:







Mean event concentration
1.2 × 10⁻¹ mg/m³
average air concentration on exposure event. Note: depends strongly on chosen exposure duration
Mean concentration on day of exposure
9.6 × 10⁻³ mg/m³
average air concentration over the day (accounts for the number of events on one day)
Year average concentration
2.6 × 10⁻⁵ mg/m³
mean daily air concentration averaged over a year
External event dose
3.8 × 10⁻² mg/kg bw
the amount that can potentially be absorbed per kg body weight during one event
External dose on day of exposure
3.8 × 10⁻² mg/kg bw
the amount that can potentially be absorbed per kg body weight during one day
Internal event dose
3.8 × 10⁻² mg/kg bw
absorbed dose per kg body weight during one exposure event
Internal dose on day of exposure
3.8 × 10⁻² mg/kg bw/day
absorbed dose per kg body weight during one day. Note: these can be higher than the ‘event dose’ for exposure frequencies larger than 1 per day.
Internal year average dose
1.1 × 10⁻⁴ mg/kg bw/day
daily absorbed dose per kg body weight averaged over a year.


I have two issues:

1. The results don't match nearly as much as I would like.
2. It can be very slow to load. I have adapted it to be quicker but still if I increase the layers to get a more accurate results it's painfully slow and not all that much more accurate.

When I compare this to the web based application I am using to get the expected results these are almost instant.

For readers familiar with PDE/ODE and Mass transfer coefficients, what is the best way to solve this?

I feel this script is perhaps a very inefficient way to reproduce the model.

Continue reading...