How to best solve partial differential equations in python accurately but quickly

Question

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 [s]
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[i] = D * (C[i + 1] - 2 * C[i] + 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 [s]

# -------------------------------
# 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 [s]

# 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 [s]
t_s_plot = np.linspace(0, total_duration_s, num=1000)      # Time vector for plotting [s]
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 [s]
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.

so 2 issues:

1. The results dont match nearly as much as I would like.
2. It's 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 painful 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.

So really I have just a few questions.

Anyone 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.

Answer 1

This may not be understood as a complete answer, as there is likely more potential, as e.g. to provide the Jacobian analytically (or e.g. by using casadi), but I doubled the speed roughly by three little things:

1. You do not need to evaluate 2000 time steps. The integrator takes care of the accuracy. The t_eval parameter only says how much you need to have reported:

   t_eval = np.linspace(0, total_emission_time_s, 50)  # 50 is enough
   

2. You can use np.empty_like to avoid initialisation, which is not necessary, as you overwrite all values subsequently:

   dCdt = np.empty_like(C)
   

3. This is the main improvement: Instead of the for-loop, do this:

   # Diffusion in the material
   dCdt[1:-1] = D * (C[2:] - 2 * C[1:-1] + C[:-2]) / dx ** 2
   

Now, I also get the same results as you (no surprise), but maybe the reference graphs are integrated less precise than yours - and therefore as well faster!?