Worked Example in Python

10th Aug 2021, Ian Fu

An worked example of sfeprapy.mcs0 module running a built-in exmaple input file.

# Dependencies, helper functions
from IPython.display import display, HTML  # Used to display DataFrame in notebook

Inputs

from sfeprapy.mcs0 import EXAMPLE_INPUT_DF  # Example input
inputs = EXAMPLE_INPUT_DF.copy()  # Make a copy

# Make some changes to the default inputs
inputs.loc['n_simulations'] = 5_000  # limit to 1,000 iterations for each case
inputs.loc['solver_thickness_ubound'] = 0.045  # Change protection max thickness to 0.045 m

display(HTML(inputs.to_html()))  # Print in HTML

	Standard Case 1	Standard Case 2 (with teq_phi)	Standard Case 3 (with timber)
case_name
n_simulations	5000	5000	5000
fire_time_step	10	10	10
fire_time_duration	18000	18000	18000
fire_hrr_density:dist	uniform_	uniform_	uniform_
fire_hrr_density:lbound	0.249	0.249	0.249
fire_hrr_density:ubound	0.251	0.251	0.251
fire_load_density:dist	gumbel_r_	gumbel_r_	gumbel_r_
fire_load_density:lbound	10	10	10
fire_load_density:ubound	1500	1500	1500
fire_load_density:mean	420	420	420
fire_load_density:sd	126	126	126
fire_spread_speed:dist	uniform_	uniform_	uniform_
fire_spread_speed:lbound	0.0035	0.0035	0.0035
fire_spread_speed:ubound	0.019	0.019	0.019
fire_nft_limit:dist	norm_	norm_	norm_
fire_nft_limit:lbound	623.15	623.15	623.15
fire_nft_limit:ubound	1473.15	1473.15	1473.15
fire_nft_limit:mean	1323.15	1323.15	1323.15
fire_nft_limit:sd	93	93	93
fire_combustion_efficiency:dist	uniform_	uniform_	uniform_
fire_combustion_efficiency:lbound	0.8	0.8	0.8
fire_combustion_efficiency:ubound	1.0	1.0	1.0
window_open_fraction:dist	lognorm_mod_	lognorm_mod_	lognorm_mod_
window_open_fraction:ubound	0.9999	0.9999	0.9999
window_open_fraction:lbound	0.0001	0.0001	0.0001
window_open_fraction:mean	0.2	0.2	0.2
window_open_fraction:sd	0.2	0.2	0.2
phi_teq:dist	constant_	lognorm_	constant_
phi_teq:ubound	1	3	1
phi_teq:lbound	1	0.00001	1
phi_teq:mean	0	1	0
phi_teq:sd	0	0.25	0
beam_cross_section_area	0.017	0.017	0.017
beam_position_horizontal:dist	uniform_	uniform_	uniform_
beam_position_horizontal:lbound	18.75	18.75	18.75
beam_position_horizontal:ubound	28.125	28.125	28.125
beam_position_vertical	3.2	3.2	3.2
beam_rho	7850	7850	7850
fire_mode	3	3	3
fire_gamma_fi_q	1	1	1
fire_t_alpha	300	300	300
fire_tlim	0.333	0.333	0.333
protection_c	1700	1700	1700
protection_k	0.2	0.2	0.2
protection_protected_perimeter	2.14	2.14	2.14
protection_rho	800	800	800
room_breadth	16	16	16
room_depth	31.25	31.25	31.25
room_height	3.3	3.3	3.3
room_wall_thermal_inertia	720	720	720
solver_temperature_goal	823.15	823.15	823.15
solver_max_iter	20	20	20
solver_thickness_lbound	0.0001	0.0001	0.0001
solver_thickness_ubound	0.045	0.045	0.045
solver_tol	1.0	1.0	1.0
window_height	2.8	2.8	2.8
window_width	72	72	72
window_open_fraction_permanent	0	0	0
timber_exposed_area	0	0	500.0
timber_charring_rate	0.7	0.7	0.7
timber_hc	13.2	13.2	13.2
timber_density	400	400	400
timber_solver_ilim	20	20	20
timber_solver_tol	1	1	1
p1	0.0	0.0	0.0
p2	0.1	0.1	0.1
p3	0.25	0.25	0.25
p4	0.09	0.09	0.09
general_room_floor_area	500	500	500

Run MCS

from sfeprapy.mcs0 import MCS0  # MCS0 class, used to take inputs, run simulation etc.

my_mcs = MCS0()
my_mcs.inputs = inputs
my_mcs.n_threads = 4
my_mcs.run_mcs()

CASE                    : Standard Case 1
NO. OF THREADS          : 4
NO. OF SIMULATIONS      : 5000

100%|█████████████████| 5000/5000 [00:03<00:00, 1658.70it/s]

fire_type               : {0: 92, 1: 4907}
beam_position_horizontal: 18.750    23.438    28.125
window_open_fraction    : 0.000     0.811     1.000
fire_combustion_efficien: 0.800     0.900     1.000
fire_hrr_density        : 0.249     0.250     0.251
fire_load_density       : 10.000    420.059   1500.000
fire_nft_limit          : 992.583   1312.528  1473.150
fire_spread_speed       : 0.004     0.011     0.019
phi_teq                 : 1.000     1.000     1.000
timber_fire_load        : 0.000     0.000     0.000

CASE                    : Standard Case 2 (with teq_phi)
NO. OF THREADS          : 4
NO. OF SIMULATIONS      : 5000

100%|█████████████████| 5000/5000 [00:03<00:00, 1657.37it/s]

fire_type               : {1: 4896, 0: 100}
beam_position_horizontal: 18.750    23.438    28.125
window_open_fraction    : 0.000     0.812     1.000
fire_combustion_efficien: 0.800     0.900     1.000
fire_hrr_density        : 0.249     0.250     0.251
fire_load_density       : 10.000    419.923   1500.000
fire_nft_limit          : 992.583   1312.503  1473.150
fire_spread_speed       : 0.004     0.011     0.019
phi_teq                 : 0.000     1.000     3.000
timber_fire_load        : 0.000     0.000     0.000

CASE                    : Standard Case 3 (with timber)
NO. OF THREADS          : 4
NO. OF SIMULATIONS      : 5000

100%|██████████████████| 5000/5000 [00:25<00:00, 199.23it/s]

fire_type               : {0: 232, 1: 3257}
beam_position_horizontal: 18.750    23.249    28.123
window_open_fraction    : 0.000     0.818     1.000
fire_combustion_efficien: 0.800     0.898     1.000
fire_hrr_density        : 0.249     0.250     0.251
fire_load_density       : 228.671   663.885   1607.459
fire_nft_limit          : 992.583   1308.322  1472.988
fire_spread_speed       : 0.004     0.012     0.019
phi_teq                 : 1.000     1.000     1.000
timber_fire_load        : 50634.284 125581.062 309776.473

All results

results = my_mcs.mcs_out
results

	case_name	n_simulations	index	beam_cross_section_area	beam_position_vertical	beam_position_horizontal	beam_rho	fire_time_duration	fire_time_step	fire_combustion_efficiency	...	solver_steel_temperature_solved	solver_time_critical_temp_solved	solver_protection_thickness	solver_iter_count	solver_time_equivalence_solved	timber_exposed_duration	timber_solver_iter_count	timber_fire_load	timber_charred_mass	timber_charred_volume
2664	Standard Case 1	5000.0	2664	0.017	3.2	22.356346	7850.0	18000.0	10.0	0.983557	...	618.737664	3990.0	-inf	0.0	-inf	0.0	0	0.000000	0.000000	0.000000
1685	Standard Case 1	5000.0	1685	0.017	3.2	25.411332	7850.0	18000.0	10.0	0.813443	...	822.647299	810.0	0.006103	7.0	1678.944357	0.0	0	0.000000	0.000000	0.000000
2986	Standard Case 1	5000.0	2986	0.017	3.2	24.059187	7850.0	18000.0	10.0	0.968234	...	822.307370	800.0	0.006141	8.0	1687.044809	0.0	0	0.000000	0.000000	0.000000
3621	Standard Case 1	5000.0	3621	0.017	3.2	19.824590	7850.0	18000.0	10.0	0.942348	...	823.591300	810.0	0.006141	5.0	1687.169770	0.0	0	0.000000	0.000000	0.000000
2568	Standard Case 1	5000.0	2568	0.017	3.2	23.982296	7850.0	18000.0	10.0	0.838248	...	823.744601	820.0	0.006246	8.0	1709.515038	0.0	0	0.000000	0.000000	0.000000
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
4991	Standard Case 3 (with timber)	5000.0	4991	0.017	3.2	27.888703	7850.0	18000.0	10.0	0.890338	...	NaN	NaN	NaN	NaN	NaN	NaN	20	164870.112663	12490.160050	31.225400
4992	Standard Case 3 (with timber)	5000.0	4992	0.017	3.2	22.701415	7850.0	18000.0	10.0	0.864813	...	NaN	NaN	NaN	NaN	NaN	NaN	20	136270.215830	10323.501199	25.808753
4993	Standard Case 3 (with timber)	5000.0	4993	0.017	3.2	20.420959	7850.0	18000.0	10.0	0.907702	...	NaN	NaN	NaN	NaN	NaN	NaN	20	111248.304573	8427.901862	21.069755
4995	Standard Case 3 (with timber)	5000.0	4995	0.017	3.2	21.148605	7850.0	18000.0	10.0	0.868454	...	NaN	NaN	NaN	NaN	NaN	NaN	20	197247.050619	14942.958380	37.357396
4997	Standard Case 3 (with timber)	5000.0	4997	0.017	3.2	26.688463	7850.0	18000.0	10.0	0.961832	...	NaN	NaN	NaN	NaN	NaN	NaN	20	128566.224690	9739.865507	24.349664

15000 rows × 60 columns

Time equivalence results

# Get time equivalence results to a dictionary
import numpy as np

dict_teq = dict()
for case in inputs.columns:
    teq = results.loc[results['case_name']==case]['solver_time_equivalence_solved']
    teq = teq[~np.isnan(teq)]  # Remove np.nan values
    teq[teq==-np.inf] = np.amin(teq[teq>-np.inf])  # Replace -np.inf values
    teq[teq==np.inf] = np.amax(teq[teq<np.inf])  # Replace np.inf values
    dict_teq[case] = teq / 60.  # seconds to minutes

# Make time equivalence CDF plots
import matplotlib.pyplot as plt
%matplotlib inline

fig, ax = plt.subplots(dpi=100)
for case in inputs.columns:
    n, bins, patches = ax.hist(dict_teq[case], 1000, density=True, histtype='step', cumulative=True, label=case)

ax.set_xlim(0, 120); ax.set_ylim(0, 1);
ax.set_xlabel('Equivlance of time exposure [$min$]'); ax.set_ylabel('CDF');
ax.legend(loc=4, fontsize='small').set_visible(True)

MCS0: Monte Carlo Simulation Method 0 MCS2: Monte Carlo Simulation Method 2