From cf1804c797bb49837d38c6dd47fedc43be895e6a Mon Sep 17 00:00:00 2001 From: fuyans Date: Fri, 23 Oct 2020 16:03:21 +0100 Subject: [PATCH] 0.8.1 (#17) This release removed quite a few modules and migrated them to another open source project `fsetools`. The purpose of `sfeprapy` is to aid the sophisticated performance based approach of probability reliability assessment (PRA) in fire safety & structural engineering. However, many helper tools which are not so relevant to PRA had been built as parts of `sfeprapy` during research and development. These helper tools lost their purpose as part of `sfeprapy` when features are improved or new technics are discovered and implemented. Thus, to provide accommodation to these helper tools as they may still be useful other than PRA, `fsetools` has been created to serve the purpose of collecting all fire safety engineering related open source codes. * New: - New `sfeprapy.func.mcs` module, a base class for all Monte Carlo Simulation classes. - New `sfeprapy.mcs2` module behaves like `sfeprapy.mcs0` but takes different input parameters. * Improved: * Removed `probability_weight`. Only `p1`, `p2`, `p3`, `p4` and `general_room_floor_area` are required. * Multiprocessing instances in `sfeprapy.func.mcs` are being reused to improve efficiency. * One way data to hard drive jobs in `sfeprapy.func.mcs` are carried out in non-blocking threads to improve its performance. * Depreciated: - Removed all GUI components. - Removed unnecessary output from `sfeprapy.mcs0`, e.g. `fire_time`, `fire_temperature` etc. - Removed figure plotting features, i.e. `all` features relevant to `plotly` and `matplotlib`. --- .travis.yml | 2 + README.md | 9 +- profiler/mcs0.py | 19 +- requirements-dev.txt | 1 + requirements.txt | 10 +- setup.py | 13 +- sfeprapy/__init__.py | 59 +- sfeprapy/cli/__main__.py | 58 +- sfeprapy/dat/ec_3_1_2_kyT.py | 146 --- sfeprapy/dat/ec_3_1_2_kyT_probability.py | 153 --- sfeprapy/dat/intumescent.py | 34 - sfeprapy/dat/steel_carbon.py | 40 - sfeprapy/func/fire_parametric_ec.py | 271 ----- sfeprapy/func/fire_parametric_ec_din.py | 375 ------ sfeprapy/func/fire_t_equare.py | 29 - sfeprapy/func/fire_travelling.py | 334 ------ sfeprapy/func/fire_travelling_flux.py | 207 ---- .../func/heat_transfer_protected_steel_ec.py | 313 ----- .../heat_transfer_unprotected_steel_ec.py | 67 -- sfeprapy/func/mcs.py | 88 +- sfeprapy/func/mcs_gen.py | 20 +- sfeprapy/func/mcs_obj.py | 70 +- sfeprapy/func/pd_6688_1_2_2007.py | 58 +- sfeprapy/func/pp.py | 182 +++ sfeprapy/gui/__init__.py | 0 sfeprapy/gui/__main__.py | 22 - sfeprapy/gui/layout/__init__.py | 0 sfeprapy/gui/layout/main.py | 90 -- sfeprapy/gui/layout/ui/main.ui | 163 --- sfeprapy/gui/layout/ui2py.py | 21 - sfeprapy/gui/logic/__init__.py | 0 sfeprapy/gui/logic/main.py | 103 -- sfeprapy/mcs0/__init__.py | 91 +- sfeprapy/mcs0/__main__.py | 170 +-- sfeprapy/mcs0/mcs0_calc.py | 1002 +++++++---------- sfeprapy/mcs1/mcs1_func_main.py | 234 ++-- sfeprapy/mcs2/__init__.py | 83 ++ sfeprapy/mcs2/__main__.py | 25 + sfeprapy/mcs2/mcs2_calc.py | 152 +++ test/mc_figures.py | 114 +- test/test_ec_3_1_2_kyT_probability.py | 3 - test/test_func_fire_parametric_ec.py | 6 - test/test_func_travelling.py | 17 - test/test_mcs0.py | 6 +- 44 files changed, 1198 insertions(+), 3662 deletions(-) create mode 100644 requirements-dev.txt delete mode 100644 sfeprapy/dat/ec_3_1_2_kyT.py delete mode 100644 sfeprapy/dat/ec_3_1_2_kyT_probability.py delete mode 100644 sfeprapy/dat/intumescent.py delete mode 100644 sfeprapy/dat/steel_carbon.py delete mode 100644 sfeprapy/func/fire_parametric_ec.py delete mode 100644 sfeprapy/func/fire_parametric_ec_din.py delete mode 100644 sfeprapy/func/fire_t_equare.py delete mode 100644 sfeprapy/func/fire_travelling.py delete mode 100644 sfeprapy/func/fire_travelling_flux.py delete mode 100644 sfeprapy/func/heat_transfer_protected_steel_ec.py delete mode 100644 sfeprapy/func/heat_transfer_unprotected_steel_ec.py create mode 100644 sfeprapy/func/pp.py delete mode 100644 sfeprapy/gui/__init__.py delete mode 100644 sfeprapy/gui/__main__.py delete mode 100644 sfeprapy/gui/layout/__init__.py delete mode 100644 sfeprapy/gui/layout/main.py delete mode 100644 sfeprapy/gui/layout/ui/main.ui delete mode 100644 sfeprapy/gui/layout/ui2py.py delete mode 100644 sfeprapy/gui/logic/__init__.py delete mode 100644 sfeprapy/gui/logic/main.py create mode 100644 sfeprapy/mcs2/__init__.py create mode 100644 sfeprapy/mcs2/__main__.py create mode 100644 sfeprapy/mcs2/mcs2_calc.py delete mode 100644 test/test_ec_3_1_2_kyT_probability.py delete mode 100644 test/test_func_fire_parametric_ec.py delete mode 100644 test/test_func_travelling.py diff --git a/.travis.yml b/.travis.yml index 3bea263..cbc9153 100644 --- a/.travis.yml +++ b/.travis.yml @@ -2,6 +2,8 @@ language: python python: - "3.7" install: + - pip install -r requirements-dev.txt + - pip install --upgrade "git+https://github.com/fsepy/fsetools.git@dev" - pip install . - pip install codecov - pip install pytest-cov diff --git a/README.md b/README.md index 13c9f7d..7448ece 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,12 @@ # SfePrapy [![GitHub version](https://badge.fury.io/gh/fsepy%2Fsfeprapy.svg)](https://github.com/fsepy/SfePrapy) -[![Updates](https://pyup.io/repos/github/fsepy/SfePrapy/shield.svg)](https://pyup.io/repos/github/fsepy/SfePrapy/) [![Build Status](https://img.shields.io/travis/fsepy/SfePrapy.svg?branch=master&label=build%20(master)&style=flat)](https://travis-ci.org/fsepy/SfePrapy) [![Build Status](https://img.shields.io/travis/fsepy/SfePrapy.svg?branch=dev&label=build%20(dev)&style=flat)](https://travis-ci.org/fsepy/SfePrapy) [![codecov](https://codecov.io/gh/fsepy/SfePrapy/branch/dev/graph/badge.svg)](https://codecov.io/gh/fsepy/SfePrapy) Structural fire engineering (Sfe) probabilistic reliability assessment (Pra) Python (py) is a probabilistic analysis tool. It calculates equivalent of time exposure to ISO 834 standard fire and this can be used to assess the appropriate fire resistance rating for structural elements using reliability based methods. -`sfeprapy` is evolving and actively used in research and real engineering design problems. +`sfeprapy` is under continuous development and actively used in research and real engineering design problems. A publication summarising the capabilities can be found [here](https://www.researchgate.net/publication/333202825_APPLICATION_OF_PYTHON_PROGRAMMING_LANGUAGE_IN_STRUCTURAL_FIRE_ENGINEERING_-_MONTE_CARLO_SIMULATION). @@ -64,12 +63,6 @@ sfeprapy mcs0 -p 4 example_input.csv `sfeprapy.mcs0` uses the [multiprocessing](https://docs.python.org/3.4/library/multiprocessing.html#module-multiprocessing) library to utilise full potential performance of multi-core CPUs. The `-p 4` defines 4 threads will be used in running the simulation, 1 is the default value. -#### To produce a figure (once a `sfeprapy.mcs0` simulation is complete) - -```sh -sfeprapy mcs0 figure mcs.out.csv -``` - ## Authors **Ian Fu** - *ian.fu@ofrconsultants.com* diff --git a/profiler/mcs0.py b/profiler/mcs0.py index 43525bf..b814d3a 100644 --- a/profiler/mcs0.py +++ b/profiler/mcs0.py @@ -4,10 +4,8 @@ def profile_standard_case(): import copy - from sfeprapy.func.mcs_obj import MCS from sfeprapy.mcs0 import EXAMPLE_INPUT_DICT, EXAMPLE_CONFIG_DICT - from sfeprapy.mcs0.mcs0_calc import teq_main, teq_main_wrapper, mcs_out_post - from sfeprapy.func.mcs_gen import main as gen + from sfeprapy.mcs0.mcs0_calc import MCS0 # increase the number of simulations so it gives sensible results mcs_input = copy.deepcopy(EXAMPLE_INPUT_DICT) @@ -18,19 +16,14 @@ def profile_standard_case(): mcs_input[k]["timber_exposed_area"] = 0 mcs_input[k].pop("beam_position_horizontal") mcs_input[k]["beam_position_horizontal:dist"] = "uniform_" - mcs_input[k]["beam_position_horizontal:ubound"] = ( - mcs_input[k]["room_depth"] * 0.9 - ) - mcs_input[k]["beam_position_horizontal:lbound"] = ( - mcs_input[k]["room_depth"] * 0.6 - ) + mcs_input[k]["beam_position_horizontal:ubound"] = mcs_input[k]["room_depth"] * 0.9 + mcs_input[k]["beam_position_horizontal:lbound"] = mcs_input[k]["room_depth"] * 0.6 # increase the number of threads so it runs faster mcs_config["n_threads"] = 3 - mcs = MCS() - mcs.define_problem(data=mcs_input, config=mcs_config) - mcs.define_stochastic_parameter_generator(gen) - mcs.define_calculation_routine(teq_main, teq_main_wrapper, mcs_out_post) + mcs = MCS0() + mcs.mcs_inputs = mcs_input + mcs.mcs_config = mcs_config mcs.run_mcs() diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 0000000..584909d --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1 @@ +cython>=0.29.21 diff --git a/requirements.txt b/requirements.txt index d9e6802..cd4fa4f 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,10 +1,8 @@ numpy>=1.17.1 -pandas>=0.25.1 -scipy>=1.3.1 -matplotlib>=3.1.1 -tqdm>=4.36 +pandas>=1.1.1 +scipy>=1.5.2 +matplotlib>=3.2.2 +tqdm>=4.50.2 xlrd>=1.2.0 -plotly>=4.2.1 docopt>=0.6 -pyside2>=5.12 openpyxl>=3.0.3 diff --git a/setup.py b/setup.py index 981e1b2..ff67b79 100644 --- a/setup.py +++ b/setup.py @@ -10,8 +10,12 @@ with open(os.path.join(os.path.dirname(os.path.realpath(__file__)), "README.md")) as f: long_description = f.read() -with open("requirements.txt") as f: - requirements = f.read().splitlines() +try: + with open("requirements.txt") as f: + requirements = f.read().splitlines() +except FileNotFoundError: + with open(os.path.join(os.path.dirname(os.path.realpath(__file__)), "requirements.txt")) as f: + requirements = f.read().splitlines() setuptools.setup( name="sfeprapy", @@ -38,12 +42,9 @@ "sfeprapy.dat", "sfeprapy.mcs0", "sfeprapy.mcs1", + "sfeprapy.mcs2", "sfeprapy.func", - "sfeprapy.dist_fit", "sfeprapy.cli", - "sfeprapy.gui", - "sfeprapy.gui.layout", - "sfeprapy.gui.logic", ], install_requires=requirements, include_package_data=True, diff --git a/sfeprapy/__init__.py b/sfeprapy/__init__.py index 2199384..eab9439 100644 --- a/sfeprapy/__init__.py +++ b/sfeprapy/__init__.py @@ -1,3 +1,47 @@ +import logging +import os + + +# setup logger +def get_logger(f_handler_fp: str = None, f_handler_level=logging.WARNING, c_handler_level=logging.INFO): + logger_ = logging.getLogger('sfeprapy') + + if f_handler_fp: + f_handler = logging.FileHandler(os.path.realpath(f_handler_fp)) + else: + f_handler = logging.FileHandler(os.path.join(os.path.expanduser('~'), 'fsetoolsgui.log')) + f_handler.setLevel(f_handler_level) + f_handler.setFormatter(logging.Formatter('%(asctime)s %(levelname)-8s [%(filename)s:%(lineno)d] %(message)s')) + logger_.addHandler(f_handler) + + c_handler = logging.StreamHandler() + c_handler.setLevel(c_handler_level) + c_handler.setFormatter(logging.Formatter('%(asctime)s %(levelname)-8s [%(filename)s:%(lineno)d] %(message)s')) + logger_.addHandler(c_handler) + + logger_.setLevel(logging.DEBUG) + + return logger_ + + +logger = get_logger() + +# make root directory of this app which will be used 1. when running the app; 2. pyinstaller at compiling the app. +if os.path.exists(os.path.dirname(__file__)): + # this path should be used when running the app as a Python package (non compiled) and/or pyinstaller at compiling + # stage. + __root_dir__ = os.path.realpath(os.path.dirname(__file__)) +elif os.path.exists(os.path.dirname(os.path.dirname(__file__))): + # the path will become invalid when the app run after compiled as the dirname `fsetoolsGUI` will disappear. + # instead, the parent folder of the project dir will be used. + __root_dir__ = os.path.realpath(os.path.dirname(os.path.dirname(__file__))) +else: + raise IsADirectoryError( + f'Project root directory undefined: ' + f'{os.path.dirname(__file__)} nor ' + f'{os.path.dirname(os.path.dirname(__file__))}' + ) + """ VERSION IDENTIFICATION RULES DOCUMENTED IN PEP 440 ARE FOLLOWED. @@ -34,11 +78,10 @@ """ -__version__ = "0.7.1.post6" +__version__ = "0.8.1" def check_pip_upgrade(): - # Parse the latest version string import subprocess from subprocess import STDOUT, check_output @@ -70,13 +113,15 @@ def check_pip_upgrade(): if __name__ == "__main__": import re + def is_canonical(version): return ( - re.match( - r"^([1-9][0-9]*!)?(0|[1-9][0-9]*)(\.(0|[1-9][0-9]*))*((a|b|rc)(0|[1-9][0-9]*))?(\.post(0|[1-9][0-9]*))?(\.dev(0|[1-9][0-9]*))?$", - version, - ) - is not None + re.match( + r"^([1-9][0-9]*!)?(0|[1-9][0-9]*)(\.(0|[1-9][0-9]*))*((a|b|rc)(0|[1-9][0-9]*))?(\.post(0|[1-9][0-9]*))?(\.dev(0|[1-9][0-9]*))?$", + version, + ) + is not None ) + assert is_canonical(__version__) diff --git a/sfeprapy/cli/__main__.py b/sfeprapy/cli/__main__.py index 7d99724..99c053f 100644 --- a/sfeprapy/cli/__main__.py +++ b/sfeprapy/cli/__main__.py @@ -1,11 +1,16 @@ """SfePrapy CLI Help. Usage: sfeprapy - sfeprapy mcs0 gui sfeprapy mcs0 run [-p=] - sfeprapy mcs0 figure sfeprapy mcs0 template - sfeprapy distfit [--data_t=] [--dist_g=] + sfeprapy mcs2 run [-p=] + sfeprapy mcs2 template + +Examples: + sfeprapy mcs0 template inputs.csv + sfeprapy mcs0 template inputs.xlsx + sfeprapy mcs0 run -p 2 inputs.csv + sfeprapy mcs0 figure mcs.out.csv Options: --data_t= an integer indicating data type: @@ -15,12 +20,10 @@ --dist_g= an integer indicating what distribution group to be used for fitting the data: 0 fit to all available distributions. 1 (default) fit to common distribution types. - -p= to define number of processes for MCS, positive integer only. + -p= to define number of processes for MCS, positive integer only. 2 by default. -h --help to show this message. Commands: - Will run `mcs0 gui` if no command is supplied. - mcs0 gui Experimental GUI version of sfeprapy. mcs0 run Monte Carlo Simulation to solve equivalent time exposure in ISO 834 fire, method 0. mcs0 figure produce figure from the output file . mcs0 template save example input file to . @@ -29,10 +32,12 @@ from docopt import docopt from sfeprapy.func.stats_dist_fit import auto_fit_2 -from sfeprapy.gui.__main__ import main as main_gui -from sfeprapy.mcs0 import EXAMPLE_INPUT_CSV +from sfeprapy.mcs0 import EXAMPLE_INPUT_CSV as EXAMPLE_INPUT_CSV_MCS0 +from sfeprapy.mcs0 import EXAMPLE_INPUT_DF as EXAMPLE_INPUT_DF_MCS0 from sfeprapy.mcs0.__main__ import main as mcs0 -from sfeprapy.mcs0.__main__ import save_figure as mcs0_figure +from sfeprapy.mcs2 import EXAMPLE_INPUT_CSV as EXAMPLE_INPUT_CSV_MCS2 +from sfeprapy.mcs2 import EXAMPLE_INPUT_DF as EXAMPLE_INPUT_DF_MCS2 +from sfeprapy.mcs2.__main__ import main as mcs2 def main(): @@ -44,42 +49,37 @@ def main(): arguments[""] = os.path.realpath(arguments[""]) if arguments["mcs0"]: - - if arguments["gui"]: - main_gui() - return 0 - - elif arguments["figure"]: - - mcs0_figure(fp_mcs0_out=arguments[""]) - - elif arguments["template"]: - from sfeprapy.mcs0 import EXAMPLE_INPUT_CSV, EXAMPLE_INPUT_DF - + if arguments["template"]: if arguments[""].endswith('.xlsx'): - EXAMPLE_INPUT_DF.to_excel(arguments[""]) + EXAMPLE_INPUT_DF_MCS0.to_excel(arguments[""]) else: with open(arguments[""], "w+", encoding='utf-8') as f: - f.write(EXAMPLE_INPUT_CSV) + f.write(EXAMPLE_INPUT_CSV_MCS0) else: - fp_mcs_in = arguments[""] n_threads = arguments["-p"] or 2 mcs0(fp_mcs_in=fp_mcs_in, n_threads=int(n_threads)) - elif arguments["distfit"]: + elif arguments['mcs2']: + if arguments['template']: + if arguments[""].endswith('.xlsx'): + EXAMPLE_INPUT_DF_MCS2.to_excel(arguments[""]) + else: + with open(arguments[""], "w+", encoding='utf-8') as f: + f.write(EXAMPLE_INPUT_CSV_MCS2) + else: + fp_mcs_in = arguments[""] + n_threads = arguments["-p"] or 2 + mcs2(fp_mcs_in=fp_mcs_in, n_threads=int(n_threads)) + elif arguments["distfit"]: # Default values data_type = arguments["--data_t"] or 2 distribution_list = arguments["--dist_g"] or 1 - # Main auto_fit_2( data_type=int(data_type), distribution_list=int(distribution_list), data=arguments[""], ) - - else: - main_gui() diff --git a/sfeprapy/dat/ec_3_1_2_kyT.py b/sfeprapy/dat/ec_3_1_2_kyT.py deleted file mode 100644 index c61e97f..0000000 --- a/sfeprapy/dat/ec_3_1_2_kyT.py +++ /dev/null @@ -1,146 +0,0 @@ -import numpy as np -import scipy.stats as stats - - -def ky2T(T): - ky = 0 - if T < 673.15: - ky = 1 - elif T <= 773.15: - ky = 1.00 + (0.78 - 1.00) / 100 * (T - 673.15) - elif T <= 873.15: - ky = 0.78 + (0.47 - 0.78) / 100 * (T - 773.15) - elif T <= 973.15: - ky = 0.47 + (0.23 - 0.47) / 100 * (T - 873.15) - elif T <= 1073.15: - ky = 0.23 + (0.11 - 0.23) / 100 * (T - 973.15) - elif T <= 1173.15: - ky = 0.11 + (0.06 - 0.11) / 100 * (T - 1073.15) - elif T <= 1273.15: - ky = 0.06 + (0.04 - 0.06) / 100 * (T - 1173.15) - elif T <= 1373.15: - ky = 0.04 + (0.02 - 0.04) / 100 * (T - 1273.15) - elif T <= 1473.15: - ky = 0.02 + (0.00 - 0.02) / 100 * (T - 1373.15) - - return ky - - -def ky2T_vectorised(T: np.ndarray): - - ky = np.where(T <= 673.15, 1, 0) - ky = np.where( - (673.15 < T) & (T <= 773.15), 1.00 + (0.78 - 1.00) / 100 * (T - 673.15), ky - ) - ky = np.where( - (773.15 < T) & (T <= 873.15), 0.78 + (0.47 - 0.78) / 100 * (T - 773.15), ky - ) - ky = np.where( - (873.15 < T) & (T <= 973.15), 0.47 + (0.23 - 0.47) / 100 * (T - 873.15), ky - ) - ky = np.where( - (973.15 < T) & (T <= 1073.15), 0.23 + (0.11 - 0.23) / 100 * (T - 973.15), ky - ) - ky = np.where( - (1073.15 < T) & (T <= 1173.15), 0.11 + (0.06 - 0.11) / 100 * (T - 1073.15), ky - ) - ky = np.where( - (1173.15 < T) & (T <= 1273.15), 0.06 + (0.04 - 0.06) / 100 * (T - 1173.15), ky - ) - ky = np.where( - (1273.15 < T) & (T <= 1373.15), 0.04 + (0.02 - 0.04) / 100 * (T - 1273.15), ky - ) - ky = np.where( - (1373.15 < T) & (T <= 1473.15), 0.02 + (0.00 - 0.02) / 100 * (T - 1373.15), ky - ) - ky = np.where(1473.15 < T, 0, ky) - - return ky - - -def ky2T_probabilistic_vectorised(T: np.ndarray, epsilon_q: np.ndarray): - - k_y_2_T_bar = ky2T_vectorised(T) - - T = T - 273.15 - - k_y_2_T_star = (k_y_2_T_bar + 1e-6) / 1.7 - - epsilon = stats.norm.ppf(epsilon_q) - - # ky2T_probabilistic_vectorised = (a1 * exp(b1 * b2 * b3 * b4 * b5)) / (exp(b1 + b2 + b3 + b4 + b5) + c1) - # ky2T_probabilistic_vectorised = (a1 * b6) / (b6 + c1) - - b1 = np.log(k_y_2_T_star / (1 - k_y_2_T_star)) - b2 = 0.412 - b3 = -0.81e-3 * T - b4 = 0.58e-6 * (T ** 1.9) - b5 = 0.43 * epsilon - b6 = np.exp(b1 + b2 + b3 + b4 + b5) - - k_y_2_T_ = (1.7 * b6) / (b6 + 1) - - return k_y_2_T_ - - -def func_test(): - - T = np.linspace(0, 1600, 100000) + 273.15 - - ky = [] - for i in T: - ky.append(ky2T(i)) - - return T, ky - - -def func_vector_test(): - - T = np.linspace(0, 1600, 100000) + 273.15 - ky = ky2T_vectorised(T) - - return T, ky - - -def func_prob_vector_test(): - - T = np.linspace(273.15 + 0, 273.15 + 1500, 5000) - q = np.random.random_sample(len(T)) - - ky = ky2T_probabilistic_vectorised(T, q) - - return T, ky - - -if __name__ == "__main__": - from timeit import default_timer as timer - import matplotlib.pyplot as plt - - t1 = timer() - x1, y1 = func_test() - t2 = timer() - x2, y2 = func_vector_test() - t3 = timer() - x3, y3 = func_prob_vector_test() - - print(t2 - t1) - print(t3 - t2) - - fig, ax3 = plt.subplots(figsize=(3.94 * 1.2, 2.76 * 1.2)) - - ax3.scatter(x3 - 273.15, y3, c="grey", s=1, label="Random Sampled Points") - ax3.plot( - x3 - 273.15, - ky2T_probabilistic_vectorised(x3, 0.5), - "--k", - label="$\epsilon$ Percentile 0.05, 0.5, 0.95", - ) - ax3.plot(x3 - 273.15, ky2T_probabilistic_vectorised(x3, 0.05), "--k") - ax3.plot(x3 - 273.15, ky2T_probabilistic_vectorised(x3, 0.95), "--k") - ax3.plot(x3 - 273.15, ky2T_vectorised(x3), "k", label="Eurocode $k_{y,\Theta}$") - ax3.set_xlabel("Temperature [$^\circ C$]") - ax3.set_ylabel("$k_{y,ach}$") - - plt.legend(loc=0, fontsize=9) - plt.tight_layout() - plt.show() diff --git a/sfeprapy/dat/ec_3_1_2_kyT_probability.py b/sfeprapy/dat/ec_3_1_2_kyT_probability.py deleted file mode 100644 index 4c8a712..0000000 --- a/sfeprapy/dat/ec_3_1_2_kyT_probability.py +++ /dev/null @@ -1,153 +0,0 @@ -import numpy as np -import scipy.stats as stats - - -def ky2T(T): - ky = 0 - if T < 673.15: - ky = 1 - elif T <= 773.15: - ky = 1.00 + (0.78 - 1.00) / 100 * (T - 673.15) - elif T <= 873.15: - ky = 0.78 + (0.47 - 0.78) / 100 * (T - 773.15) - elif T <= 973.15: - ky = 0.47 + (0.23 - 0.47) / 100 * (T - 873.15) - elif T <= 1073.15: - ky = 0.23 + (0.11 - 0.23) / 100 * (T - 973.15) - elif T <= 1173.15: - ky = 0.11 + (0.06 - 0.11) / 100 * (T - 1073.15) - elif T <= 1273.15: - ky = 0.06 + (0.04 - 0.06) / 100 * (T - 1173.15) - elif T <= 1373.15: - ky = 0.04 + (0.02 - 0.04) / 100 * (T - 1273.15) - elif T <= 1473.15: - ky = 0.02 + (0.00 - 0.02) / 100 * (T - 1373.15) - - return ky - - -def ky2T_vectorised(T: np.ndarray): - - ky = np.where(T <= 673.15, 1, 0) - ky = np.where( - (673.15 < T) & (T <= 773.15), 1.00 + (0.78 - 1.00) / 100 * (T - 673.15), ky - ) - ky = np.where( - (773.15 < T) & (T <= 873.15), 0.78 + (0.47 - 0.78) / 100 * (T - 773.15), ky - ) - ky = np.where( - (873.15 < T) & (T <= 973.15), 0.47 + (0.23 - 0.47) / 100 * (T - 873.15), ky - ) - ky = np.where( - (973.15 < T) & (T <= 1073.15), 0.23 + (0.11 - 0.23) / 100 * (T - 973.15), ky - ) - ky = np.where( - (1073.15 < T) & (T <= 1173.15), 0.11 + (0.06 - 0.11) / 100 * (T - 1073.15), ky - ) - ky = np.where( - (1173.15 < T) & (T <= 1273.15), 0.06 + (0.04 - 0.06) / 100 * (T - 1173.15), ky - ) - ky = np.where( - (1273.15 < T) & (T <= 1373.15), 0.04 + (0.02 - 0.04) / 100 * (T - 1273.15), ky - ) - ky = np.where( - (1373.15 < T) & (T <= 1473.15), 0.02 + (0.00 - 0.02) / 100 * (T - 1373.15), ky - ) - ky = np.where(1473.15 < T, 0, ky) - - return ky - - -def ky2T_probabilistic_vectorised(T: np.ndarray, epsilon_q: np.ndarray): - - k_y_2_T_bar = ky2T_vectorised(T) - - k_y_2_T_star = (k_y_2_T_bar + 1e-6) / 1.7 - - epsilon = stats.norm.ppf(epsilon_q) - - # ky2T_probabilistic_vectorised = (a1 * exp(b1 * b2 * b3 * b4 * b5)) / (exp(b1 + b2 + b3 + b4 + b5) + c1) - # ky2T_probabilistic_vectorised = (a1 * b6) / (b6 + c1) - - b1 = np.log(k_y_2_T_star / (1 - k_y_2_T_star)) - b2 = 0.412 - b3 = -0.81e-3 * T - b4 = 0.58e-6 * (T ** 1.9) - b5 = 0.43 * epsilon - b6 = np.exp(b1 + b2 + b3 + b4 + b5) - - k_y_2_T_ = (1.7 * b6) / (b6 + 1) - - return k_y_2_T_ - - -def func_test(): - - T = np.linspace(0, 1600, 100000) + 273.15 - - ky = [] - for i in T: - ky.append(ky2T(i)) - - return T, ky - - -def func_vector_test(): - - T = np.linspace(0, 1600, 100000) + 273.15 - ky = ky2T_vectorised(T) - - return T, ky - - -def func_prob_vector_test(): - - T = np.linspace(273.15 + 0, 273.15 + 1500, 5000) - q = np.random.random_sample(len(T)) - - ky = ky2T_probabilistic_vectorised(T, q) - - return T, ky - - -def _test_probabilistic(): - assert abs(ky2T_probabilistic_vectorised(0, 0.5) - 1.161499) <= 0.00001 - assert abs(ky2T_probabilistic_vectorised(673.15, 0.5) - 1.001560) <= 0.0001 - - -if __name__ == "__main__": - _test_probabilistic() - from timeit import default_timer as timer - import matplotlib.pyplot as plt - - t1 = timer() - x1, y1 = func_test() - t2 = timer() - x2, y2 = func_vector_test() - t3 = timer() - x3, y3 = func_prob_vector_test() - - print(t2 - t1) - print(t3 - t2) - - fig, ax3 = plt.subplots(figsize=(3.94 * 1.2, 2.76 * 1.2)) - - ax3.scatter(x3 - 273.15, y3, c="grey", s=1, label="Random Sampled Points") - ax3.plot( - x3 - 273.15, - ky2T_probabilistic_vectorised(x3, 0.5), - "--k", - label="$\epsilon$ Percentile 0.05, 0.5, 0.95", - ) - - print(ky2T_probabilistic_vectorised(673.15, 0.5)) - - ax3.plot(x3 - 273.15, ky2T_probabilistic_vectorised(x3, 0.05), "--k") - ax3.plot(x3 - 273.15, ky2T_probabilistic_vectorised(x3, 0.95), "--k") - ax3.plot(x3 - 273.15, ky2T_vectorised(x3), "k", label="Eurocode $k_{y,\Theta}$") - ax3.set_xlabel("Temperature [$^\circ C$]") - ax3.set_ylabel("$k_{y,ach}$") - - plt.legend(loc=0, fontsize=9) - plt.tight_layout() - plt.show() diff --git a/sfeprapy/dat/intumescent.py b/sfeprapy/dat/intumescent.py deleted file mode 100644 index 63e23ad..0000000 --- a/sfeprapy/dat/intumescent.py +++ /dev/null @@ -1,34 +0,0 @@ -# -*- coding: utf-8 -*- -import numpy as np -import pandas as pd -from scipy.interpolate import interp1d -import os - - -def thermal(property_name): - dir_package = os.path.dirname(os.path.abspath(__file__)) - - dir_files = { - "thermal conductivity (thomas2002)": "k_1_T_intumescent_thomas2002.csv" - } - - # read file - data = pd.read_csv( - "/".join([dir_package, dir_files[property_name]]), delimiter=",", dtype=float - ) - x, y = data.values[:, 0], data.values[:, 1] - - return interp1d(x, y) - - -# if __name__ == "__main__": -# temperature = np.arange(20.,1200.+0.5,0.5) + 273.15 -# prop = temperature * 0 -# for i,v in enumerate(temperature): -# prop[i] = steel_specific_heat_carbon_steel(v) -# df = pd.DataFrame( -# { -# "Temperature [K]": temperature, -# "Specific Heat Capacity [J/kg/K]": prop -# } -# ) diff --git a/sfeprapy/dat/steel_carbon.py b/sfeprapy/dat/steel_carbon.py deleted file mode 100644 index ac05866..0000000 --- a/sfeprapy/dat/steel_carbon.py +++ /dev/null @@ -1,40 +0,0 @@ -# -*- coding: utf-8 -*- -import numpy as np - - -def steel_specific_heat_carbon_steel(temperature: float) -> float: - """Returns steel specific heat in accordance with EN 1993-1-2. SI units. - - :param temperature: [K] - :return: [J/kg/K] - """ - - temperature -= 273.15 # unit conversion, K -> deg. C - - if 20 <= temperature < 600: - return ( - 425 - + 0.773 * temperature - - 1.69e-3 * np.power(temperature, 2) - + 2.22e-6 * np.power(temperature, 3) - ) - elif 600 <= temperature < 735: - return 666 + 13002 / (738 - temperature) - elif 735 <= temperature < 900: - return 545 + 17820 / (temperature - 731) - elif 900 <= temperature <= 1200: - return 650 - else: - return 0 - - -if __name__ == "__main__": - - import timeit - - test_number = 100000 - - def aa(): - steel_specific_heat_carbon_steel(np.random.rand() * 1180 + 20 + 273.15) - - print(timeit.timeit(aa, number=test_number)) diff --git a/sfeprapy/func/fire_parametric_ec.py b/sfeprapy/func/fire_parametric_ec.py deleted file mode 100644 index b625f7f..0000000 --- a/sfeprapy/func/fire_parametric_ec.py +++ /dev/null @@ -1,271 +0,0 @@ -import numpy as np - - -def fire( - t, A_t, A_f, A_v, h_eq, q_fd, lambda_, rho, c, t_lim, temperature_initial=293.15 -): - """Function Description: (SI UNITS ONLY) - This function calculates the time-temperature curve according to Eurocode 1 part 1-2, Appendix A. - :param t: numpy.ndarray, [s], time evolution. - :param A_t: float, [m2], total surface area (including openings). - :param A_f: float, [m2], floor area. - :param A_v: float, [m2], opening area. - :param h_eq: float, [m2], opening height. - :param q_fd: float, [J/m2], fuel density. - :param lambda_: float, [K/kg/m], lining thermal conductivity. - :param rho: float, [kg/m3], lining density. - :param c: float, [J/K/kg], lining thermal capacity. - :param t_lim: float, [s], limiting time for the fire. - :return T_g: numpy.ndarray, [s], temperature evolution. - """ - # Reference: Eurocode 1991-1-2; Jean-Marc Franssen, Paulo Vila Real (2010) - Fire Design of Steel Structures - - # UNITS: SI -> Equations - q_fd /= 1e6 # [J/m2] -> [MJ/m2] - t_lim /= 3600 # [s] -> [hr] - t = t / 3600 # [s] -> [hr] - temperature_initial -= 273.15 # [K] -> [C] - - # ACQUIRING REQUIRED VARIABLES - # t = np.arange(time_start, time_end, time_step, dtype=float) - - b = (lambda_ * rho * c) ** 0.5 - O = A_v * h_eq ** 0.5 / A_t - q_td = q_fd * A_f / A_t - Gamma = ((O / 0.04) / (b / 1160)) ** 2 - - t_max = 0.0002 * q_td / O - - # check criteria - # todo: raise warnings to higher level - # if not 50 <= q_td <= 1000: - # if is_cap_q_td: - # msg = "q_td ({:4.1f}) EXCEEDED [50, 1000] AND IS CAPPED.".format(q_td, is_cap_q_td) - # else: - # msg = "q_td ({:4.1f}) EXCEEDED [50, 1000] AND IS UNCAPPED.".format(q_td) - # warnings.warn(msg) - - # CALCULATION - def _temperature_heating(t_star, temperature_initial): - # eq. 3.12 - T_g = 1325 * ( - 1 - - 0.324 * np.exp(-0.2 * t_star) - - 0.204 * np.exp(-1.7 * t_star) - - 0.472 * np.exp(-19 * t_star) - ) - T_g += temperature_initial - return T_g - - def _temperature_cooling_vent(t_star_max, T_max, t_star): # ventilation controlled - # eq. 3.16 - if t_star_max <= 0.5: - T_g = T_max - 625 * (t_star - t_star_max) - elif 0.5 < t_star_max < 2.0: - T_g = T_max - 250 * (3 - t_star_max) * (t_star - t_star_max) - elif 2.0 <= t_star_max: - T_g = T_max - 250 * (t_star - t_star_max) - else: - T_g = np.nan - return T_g - - def _temperature_cooling_fuel( - t_star_max, T_max, t_star, Gamma, t_lim - ): # fuel controlled - # eq. 3.22 - if t_star_max <= 0.5: - T_g = T_max - 625 * (t_star - Gamma * t_lim) - elif 0.5 < t_star_max < 2.0: - T_g = T_max - 250 * (3 - t_star_max) * (t_star - Gamma * t_lim) - elif 2.0 <= t_star_max: - T_g = T_max - 250 * (t_star - Gamma * t_lim) - else: - T_g = np.nan - return T_g - - def _variables(t, Gamma, t_max): - t_star = Gamma * t - t_star_max = Gamma * t_max - return t_star, t_star_max - - def _variables_2(t, t_lim, q_td, b, O): - O_lim = 0.0001 * q_td / t_lim - Gamma_lim = ((O_lim / 0.04) / (b / 1160)) ** 2 - - if O > 0.04 and q_td < 75 and b < 1160: - k = 1 + ((O - 0.04) / (0.04)) * ((q_td - 75) / (75)) * ((1160 - b) / (1160)) - Gamma_lim *= k - - t_star_ = Gamma_lim * t - t_star_max_ = Gamma_lim * t_lim - return t_star_, t_star_max_ - - t_star, t_star_max = _variables(t, Gamma, t_max) - - if t_max >= t_lim: # ventilation controlled fire - T_max = _temperature_heating(t_star_max, temperature_initial) - T_heating_g = _temperature_heating(Gamma * t, temperature_initial) - T_cooling_g = _temperature_cooling_vent(t_star_max, T_max, t_star) - fire_type = "ventilation controlled" - else: # fuel controlled fire - t_star_, t_star_max_ = _variables_2(t, t_lim, q_td, b, O) - T_max = _temperature_heating(t_star_max_, temperature_initial) - T_heating_g = _temperature_heating(t_star_, temperature_initial) - T_cooling_g = _temperature_cooling_fuel(t_star_max, T_max, t_star, Gamma, t_lim) - fire_type = "fuel controlled" - - T_g = np.minimum(T_heating_g, T_cooling_g) - T_g[T_g < temperature_initial] = temperature_initial - - # UNITS: Eq. -> SI - t *= 3600 - T_g += 273.15 - - return T_g - - -def example_plot_interflam(): - import matplotlib.pyplot as plt - import seaborn as sns - - sns.set_style("ticks") - sns.set_context(context="paper") - - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel(r"Temperature [$^{\circ}C$]") - - # define geometry - w = 16 - l = 32 - h = 3 - - h_v = 2 - - oo = (0.02, 0.04, 0.06, 0.10, 0.14, 0.20) # desired opening factor - - def opening_factor_2_window_width(of, w, l, h, h_v): - # opening factor = Av * sqrt(hv) / At - # Av = h_v * w_v - # of * At / sqrt(h_v) = h_v * w_v - # w_v = of * At / sqrt(h_v) / h_v - return of * (2 * (w * l + l * h + h * w)) / h_v ** 0.5 / h_v - - w_v_ = [opening_factor_2_window_width(o, w, l, h, h_v) for o in oo] - - print(w_v_) - - # Calculate fire curves - for w_v in w_v_: - x = np.arange(0, 5 * 60 * 60, 1) - y = fire( - t=x, - A_v=h_v * w_v, - A_t=2 * (w * h + h * l + l * w), - A_f=w * l, - h_eq=h_v, - q_fd=600e6, - lambda_=1, - rho=1, - c=720 ** 2, - t_lim=20 * 60, - temperature_initial=293.15, - ) - ax.plot( - x / 60, - y - 273.15, - label="Opening Factor {:.2f}".format( - (h_v * w_v) * h_v ** 0.5 / (2 * (w * h + h * l + l * w)) - ), - ) - - ax.legend().set_visible(True) - ax.grid(color="k", linestyle="--") - ax.set_ylim((0, 1400)) - ax.set_xlim((0, 300)) - plt.tight_layout() - - plt.savefig(fname="fire-ec.png", dpi=300) - plt.show() - - -def _test_fire(): - """This is a test function to `fire` within this module `fire_parametric_ec`, it compares the function against - Figure 7 in Holicky, M. et al [1]. - yan fu, 1 oct 2018 - - REFERENCES - [1] Holicky, M., Meterna, A., Sedlacek, G. and Schleich, J.B., 2005. Implementation of eurocodes, handbook 5, design - of buildings for the fire situation. Leonardo da Vinci Pilot Project: Luxembourg.""" - - import copy - import pandas - - # LOAD VERIFICATION DATA - # data are from the referenced document - from io import StringIO - - verification_data = StringIO( - """time_1,temperature_1,time_2,temperature_2,time_3,temperature_3,time_4,temperature_4,time_5,temperature_5,time_6,temperature_6,time_7,temperature_7 - 20,689.1109391,20,689.1109391,20,689.1109391,20,843.9616854,20,741.3497451,20,601.424372,20,269.3348197 - 30,782.7282474,30,782.7317977,30,782.7317977,30,946.9144525,30,827.5114674,30,726.7651987,30,430.1234077 - 40,167.4009827,40,402.4756096,40,508.8188932,40,781.2034026,40,885.6845648,40,777.4756096,40,533.0726245 - 50,20,50,20,50,236.7716603,50,540.8761379,50,818.8612125,50,813.2641976,50,604.3089737 - 60,20,60,20,60,20,60,304.2731159,60,692.3328174,60,843.4522203,60,651.2880412 - 70,20,70,20,70,20,70,69.53931579,70,565.8079725,70,780.3602113,70,685.2109576 - 80,20,80,20,80,20,80,20,80,441.145249,80,719.1303236,80,709.8019654""" - ) - verification_data = pandas.read_csv( - verification_data, skip_blank_lines=True, skipinitialspace=True - ) - - # CALCULATE TIME TEMPERATURE CURVE BASED ON THE VERIFICATION DATA INPUTS - - # prepare inputs - kws = dict( - A_t=360, - A_f=100, - h_eq=1, - q_fd=600e6, - lambda_=1, - rho=1, - c=2250000, - t_lim=20 * 60, - t=np.arange(0, 2 * 60 * 60 + 1, 1), - temperature_initial=293.15, - ) - - # define opening area - A_v_list = [72, 50.4, 36.000000001, 32.4, 21.6, 14.4, 7.2] - - # calculate fire curves - x1_list = [] # calculated time array - y1_list = [] # calculated temperature array - for i in A_v_list: - y = fire(A_v=i, **copy.copy(kws)) - x = np.arange(0, 2 * 60 * 60 + 1, 1) + 10 * 60 - x1_list.append(x / 60) - y1_list.append(y - 273.15) - - # COMPARE CALCULATED AND THE VERIFICATION DATA - - def r_square(x1, y1, x2, y2): - from scipy.interpolate import interp1d - from scipy import stats - - slope, intercept, r_value, p_value, std_err = stats.linregress( - y2, interp1d(x1, y1)(x2) - ) - return r_value ** 2 - - for i in range(len(A_v_list)): - x1 = x1_list[i] - y1 = y1_list[i] - x2 = verification_data["time_{}".format(int(i + 1))] - y2 = verification_data["temperature_{}".format(int(i + 1))] - assert r_square(x1, y1, x2, y2) > 0.99 - - -if __name__ == "__main__": - # exmaple_plot_interflam() - _test_fire() diff --git a/sfeprapy/func/fire_parametric_ec_din.py b/sfeprapy/func/fire_parametric_ec_din.py deleted file mode 100644 index b50061b..0000000 --- a/sfeprapy/func/fire_parametric_ec_din.py +++ /dev/null @@ -1,375 +0,0 @@ -# -*- coding: utf-8 -*- -import warnings - -import numpy as np - - -def fire( - t_array_s, - A_w_m2, - h_w_m2, - A_t_m2, - A_f_m2, - t_alpha_s, - b_Jm2s05K, - q_x_d_MJm2, - gamma_fi_Q=1.0, - q_ref=1300, - alpha=0.0117, - hrrpua_MWm2=0.25, -): - """This piece of code calculates a time dependent temperature array in accordance of Appendix AA in German Annex " - Simplified natural fire model for fully developed room fires" to Eurocode 1991-1-2 (DIN EN 1991-1-2/NA:2010-12). - Limitations are detailed in the Section AA.2 and they are: - - minimum fuel load 100 MJ/sq.m - - maximum fuel load 1,300 MJ/sq.m - - maximum floor area 400 sq.m - - maximum ceiling height 5 m - - minimum vertical vent opening to floor area 12.5 % - - maximum vertical vent opening to floor area 50 % - - PARAMETERS: - :param t_array_s: numpy.array, in [s], time - :param A_w_m2: float, [m2], is window opening area - :param h_w_m2: float, [m2], is weighted window opening height - :param A_t_m2: float, [m2], is total enclosure internal surface area, including openings - :param A_f_m2: float, [m2], is total floor area - :param t_alpha_s: float, [s], is the fire growth factor, Table BB.2, for residential/office t_alpha = 300 [s] - :param b_Jm2s05K: float, [J/m2/s0.5/K], is the weighted heat storage capacity, see equation AA.31 and Table AA.1 - :param q_x_d_MJm2: float, [MJ/m2], is the design value for fire load density, same as q_f_d - :param gamma_fi_Q: float, is the partial factor according to BB.5.3 - :param q_ref: float, [MJ/m2], is the reference upper bound heat release rate, 1300 [MJ/m2] for offices and residential - :param alpha: float, [kW/s2], is the growth constant for t-square fire - :param hrrpua_MWm2: float, [MW], is the heat release rate per unit area - :return T: numpy.ndarray, in [K], is calculated gas temperature within fire enclosure - - SUPPLEMENTAL DATA: - - DIN EN 1991-1-2/NA:1010-12 (English translation) - - Table AA.1 - Influence groups as a function of heat storage capacity b - | Line | Influence group | Heat storage capacity b [J/m2/s0.5/K] | - |------|-----------------|---------------------------------------| - | 1 | 1 | 2500 | - | 2 | 2 | 1500 | - | 3 | 3 | 750 | - - Table BB.2 - | Line | Occupancy | Fire propagation | t_alpha [s] | RHR_f [MW/m2] | - |------|--------------------------------------------|------------------|-------------|---------------| - | 1 | Residential building | medium | 300 | 0.25 | - | 2 | Office building | medium | 300 | 0.25 | - | 3 | Hospital (room) | medium | 300 | 0.25 | - | 4 | Hotel (room) | medium | 450 | 0.25 | - | 5 | Library | medium | 300 | 0.25 | - | 6 | School (classroom) | medium | 300 | 0.15 | - | 7 | Shop, shopping centre | fast | 150 | 0.25 | - | 8 | Place of public assembly (theatre, cinema) | fast | 150 | 0.50 | - | 9 | Public transport | slow | 600 | 0.25 | - """ - - # AA.1 - Q_max_v_k = 1.21 * A_w_m2 * np.sqrt(h_w_m2) # [MW] AA.1 - - # AA.2 - Q_max_f_k = hrrpua_MWm2 * A_f_m2 # [MW] AA.2 - - # AA.3, Characteristic value of the maximum HRR, is the lower value of the Q_max_f_k and Q_max_v_k - Q_max_k = min(Q_max_f_k, Q_max_v_k) # [MW] - - # AA.5 - Q_max_v_d = gamma_fi_Q * Q_max_v_k # [MW] AA.5 - Q_max_f_d = gamma_fi_Q * Q_max_f_k # [MW] AA.6, gamma_fi_Q see BB.5.3 - - # AA.6 - Q_max_d = gamma_fi_Q * Q_max_k - - # Work out fire type - if Q_max_v_k == Q_max_k: - fire_type = 0 - elif Q_max_f_k == Q_max_k: - fire_type = 1 - else: - fire_type = np.nan - - # print(fire_type) - - # AA.7 - AA.19: Calculate location of t and Q - O = A_w_m2 * h_w_m2 ** 0.5 / A_t_m2 - Q_d = q_ref * A_f_m2 # [MJ], total fire load in the compartment - if fire_type == 0: # ventilation controlled fire - - # AA.7 - t_1 = t_alpha_s * np.sqrt(Q_max_v_d) # [s] AA.7 - - # AA.8 - T_1_v = -8.75 / O - 0.1 * b_Jm2s05K + 1175 # [deg.C] - - Q_1 = t_1 ** 3 / (3 * t_alpha_s ** 2) - - # AA.9 - Q_2 = 0.7 * Q_d - Q_1 # AA.9(a) - t_2 = t_1 + Q_2 / Q_max_v_d # [s] AA.9(b) - - # AA.10 - T_2_v = min( - 1134, (0.004 * b_Jm2s05K - 17) / O - 0.4 * b_Jm2s05K + 2175 - ) # [deg.C] AA.10 - - # AA.11 - Q_3 = 0.3 * Q_d - t_3 = t_2 + (2 * Q_3) - - # AA.12 - T_3_v = -5.0 / O - 0.16 * b_Jm2s05K + 1060 # [deg.C] - - T_1, T_2, T_3 = T_1_v, T_2_v, T_3_v - - elif fire_type == 1: # fuel controlled fire - - # AA.19 - k = ( - (Q_max_f_d ** 2) / (A_w_m2 * h_w_m2 ** 0.5 * (A_t_m2 - A_w_m2) * b_Jm2s05K) - ) ** (1 / 3) - # print(h_w_m2, A_t_m2, A_w_m2, A_w_m2 * h_w_m2**0.5 * (A_t_m2-A_w_m2) * b_Jm2s05K) - - # AA.13 - t_1 = t_alpha_s * Q_max_f_d ** 0.5 - Q_1 = 1 / 3 * alpha * t_alpha_s ** 3 - Q_1 *= 1e-3 - - # AA.14 - T_1_f = min(980, 24000 * k + 20) # [deg.C] - - # AA.15 - Q_2 = 0.7 * Q_d - Q_1 - t_2 = t_1 + Q_2 / Q_max_f_d - - # AA.16 - T_2_f = min(1340, 33000 * k + 20) # [deg.C] - - # AA.17 - Q_3 = 0.3 * Q_d - t_3 = t_2 + (2 * Q_3) / Q_max_f_d - - # AA.18 - T_3_f = min(660, 16000 * k + 20) # [deg.C] - - # AA.19 - # See above - - T_1, T_2, T_3 = T_1_f, T_2_f, T_3_f - - else: - wmsg = "WTH, I do not know what type of fire ({}) this is, bugs?".format( - fire_type - ) - warnings.warn(wmsg) - k = t_1 = T_1 = Q_2 = t_2 = Q_2_f = Q_3 = t_3 = Q_3_f = T_2 = T_3 = np.nan - - # Prerequisite for AA.20 and AA.21 - Q_1 = t_1 ** 3 / (3 * t_alpha_s ** 2) # [MW] - Q_x_d = q_x_d_MJm2 * A_f_m2 - - if Q_1 < 0.7 * Q_x_d: - # AA.20 - t_2_x = t_1 + (0.7 * Q_x_d - t_1 ** 3 / (3 * t_alpha_s ** 2)) / Q_max_d # [s] - - # AA.21 - T_2_x = (T_2 - T_1) * ((t_2_x - t_1) / (t_2 - t_1)) ** 0.5 + T_1 # [deg.C] - elif Q_1 >= 0.7: - # AA.22 - t_1_x = (0.7 * Q_x_d * 3 * t_alpha_s ** 2) ** (1 / 3) # [s] - t_2_x = (0.7 * Q_x_d * 3 * t_alpha_s ** 2) ** (1 / 3) # [s] - - # AA.23 - T_2_x = (T_1 - 20) / (t_1 ** 2) * t_1_x ** 2 + 20 # [deg.C] - else: - warnings.warn("Q_1 out of bound for AA.20 to AA.23.") - t_1_x = t_2_x = T_2_x = np.nan - - # AA.25 - t_3_x = 0.6 * Q_x_d / Q_max_d + t_2_x # [s] - - # AA.24 - T_3_x = T_3 * np.log10(t_3_x / 60 + 1) / np.log10(t_3 / 60 + 1) # [deg.C] - - T = T_t( - t_array_s, - t_1, - t_2, - t_2_x, - t_3_x, - T_1, - T_2_x, - T_3_x, - A_t_m2, - A_w_m2, - h_w_m2, - t_alpha_s, - ) - - # CONVERT UNITS TO SI - T += 273.15 - - return T - - -# AA.26 - AA.28 -def T_t( - t, t_1, t_2, t_2_x, t_3_x, T_1, T_2_x, T_3_x, A_t, A_w, h_w, t_alpha, T_initial=20 -): - - # Initialise container for return value - T = np.zeros(len(t)) - - # Check flash-over - # AA.30 - Q_fo = 0.0078 * A_t + 0.378 * A_w * h_w ** 0.5 # [MW] - # AA.29 - t_1_fo = (t_alpha ** 2 * Q_fo) ** 0.5 # [s] - - t_1 = min(t_1, t_1_fo) - - # AA.26 - t_1_ = np.logical_and(0 <= t, t <= t_1) - - T_1_ = (T_1 - 20) / t_1 ** 2 * t[t_1_] ** 2 + 20 - T[t_1_] = T_1_ # [deg.C] - - # AA.27 - # t_2_ = np.logical_and(t_1 <= t, t <= t_2) - # T_2_ = (T_2_x - T_1) * ((t[t_2_] - t_1) / (t_2_x - t_1)) ** 0.5 + T_1 - # T[t_2_] = T_2_ # [deg.C] - - # AA.27 MOD - t_2_ = np.logical_and(t_1 <= t, t <= t_2_x) - T_2_ = (T_2_x - T_1) * ((t[t_2_] - t_1) / (t_2_x - t_1)) ** 0.5 + T_1 - T[t_2_] = T_2_ # [deg.C] - - # AA.28 - # t_3_ = t > t_2 - # T_3_ = (T_3_x - T_2_x) * ((t[t_3_] - t_2) / (t_3_x - t_2_x)) ** 0.5 + T_2_x - # T[t_3_] = T_3_ # [deg.C] - - # AA.28 MOD - t_3_ = t > t_2_x - T_3_ = (T_3_x - T_2_x) * ((t[t_3_] - t_2_x) / (t_3_x - t_2_x)) ** 0.5 + T_2_x - T[t_3_] = T_3_ # [deg.C] - - # No temperate below T_initial - T[T < T_initial] = T_initial - - return T - - -def example_plot_interflam(): - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel("Temperature [$^{℃}C$]") - - # define geometry - w = 16 - l = 32 - h = 3 - - h_v = 2 - - oo = (0.02, 0.04, 0.06, 0.10, 0.14, 0.20) # desired opening factor - - def of2wv(of, w, l, h, h_v): - # opening factor = Av * sqrt(hv) / At - # Av = h_v * w_v - # of * At / sqrt(h_v) = h_v * w_v - # w_v = of * At / sqrt(h_v) / h_v - return of * (2 * (w * l + l * h + h * w)) / h_v ** 0.5 / h_v - - w_v_ = [of2wv(o, w, l, h, h_v) for o in oo] - - print(w_v_) - - for w_v in w_v_: - q_ref = 1300 - x = np.arange(0, 5 * 60 * 60 + 1, 1) - y = fire( - t_array_s=x, - A_w_m2=h_v * w_v, - h_w_m2=h_v, - A_t_m2=2 * (l * w + w * h + h * l), - A_f_m2=w * l, - t_alpha_s=300, - b_Jm2s05K=720, - q_x_d_MJm2=600, - gamma_fi_Q=1.0, - q_ref=q_ref, - ) - ax.plot( - x / 60, - y - 273.15, - label="Opening Factor {:.2f}".format( - (h_v * w_v) * h_v ** 0.5 / (2 * (w * h + h * l + l * w)) - ), - ) - - ax.legend().set_visible(True) - ax.grid(color="k", linestyle="--") - ax.set_ylim((0, 1400)) - ax.set_xlim((0, 300)) - plt.tight_layout() - plt.savefig(fname="fire-ec_din.png", dpi=300) - - -def plot_variable_qfd(): - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel("Temperature [$^{℃}C$]") - - # define geometry - - w = 16 - l = 32 - h = 3 - - h_v = 2 - w_v = 5 - - q_fd_ = [150, 300, 450, 600] - - for q_fd in q_fd_: - q_ref = 1300 - x = np.arange(0, 5 * 60 * 60 + 1, 1) - y = fire( - t_array_s=x, - A_w_m2=h_v * w_v, - h_w_m2=h_v, - A_t_m2=2 * (w * l + l * h + h * w), # 80, - A_f_m2=w * l, - t_alpha_s=300, - b_Jm2s05K=750, - q_x_d_MJm2=q_fd, - gamma_fi_Q=1.0, - q_ref=q_ref, - ) - ax.plot( - x / 60, - y - 273.15, - label="Fuel load density {q_fd:4.0f} $MJ/m^2$".format(q_fd=q_fd), - ) - - ax.legend(loc=4).set_visible(True) - ax.grid(color="k", linestyle="--") - ax.set_ylim((-10, 1310)) - plt.tight_layout() - plt.savefig(fname="fire-ec_din.png", dpi=300) - plt.show() - - -if __name__ == "__main__": - # plot_variable_qfd() - example_plot_interflam() diff --git a/sfeprapy/func/fire_t_equare.py b/sfeprapy/func/fire_t_equare.py deleted file mode 100644 index d78269f..0000000 --- a/sfeprapy/func/fire_t_equare.py +++ /dev/null @@ -1,29 +0,0 @@ -def fire(time, growth_factor, cap_hrr=0, cap_hrr_to_time=0): - if isinstance(growth_factor, str): - growth_factor = str.lower(growth_factor) - growth_dict = { - "slow": 0.0029, - "medium": 0.0117, - "fast": 0.0469, - "ultra-fast": 0.1876, - } - try: - growth_factor = growth_dict[growth_factor] - except KeyError: - err_msg = "{} should be one of the following if not a number: {}" - err_msg = err_msg.format( - "growth_factor", ", ".join(list(growth_dict.keys())) - ) - raise ValueError(err_msg) - - heat_release_rate = growth_factor * time ** 2 * 1000 - - # cap hrr - if cap_hrr: - heat_release_rate[heat_release_rate > cap_hrr] = cap_hrr - - if cap_hrr_to_time: - heat_release_rate[time > cap_hrr_to_time] = -1 - heat_release_rate[heat_release_rate == -1] = max(heat_release_rate) - - return heat_release_rate diff --git a/sfeprapy/func/fire_travelling.py b/sfeprapy/func/fire_travelling.py deleted file mode 100644 index f74b4f0..0000000 --- a/sfeprapy/func/fire_travelling.py +++ /dev/null @@ -1,334 +0,0 @@ -# -*- coding: utf-8 -*- -import copy -from typing import Union - -import numpy as np - - -def fire( - t: np.array, - fire_load_density_MJm2: float, - fire_hrr_density_MWm2: float, - room_length_m: float, - room_width_m: float, - fire_spread_rate_ms: float, - beam_location_height_m: float, - beam_location_length_m: Union[float, list, np.ndarray], - fire_nft_limit_c: float, - opening_fraction: float = 0, - opening_width_m: float = 0, - opening_height_m: float = 0, -): - """ - This function calculates and returns a temperature array representing travelling fire. This function is NOT in SI. - :param t: in s, is the time array - :param fire_load_density_MJm2: in MJ/m2, is the fuel density on the floor - :param fire_hrr_density_MWm2: in MW/m2, is the heat release rate density - :param room_length_m: in m, is the room length - :param room_width_m: in m, is the room width - :param fire_spread_rate_ms: in m/s, is fire spread speed - :param beam_location_height_m: in m, is the beam lateral distance to fire origin - :param beam_location_length_m: in m, is the beam height above the floor - :param fire_nft_limit_c: in deg.C, is the maximum near field temperature - :param opening_fraction: in -, is the ventilation opening proportion between 0 to 1 - :param opening_width_m: in m, is ventilation opening width - :param opening_height_m: in m, is ventilation opening height - :return T_g: in deg.C, is calculated gas temperature - """ - - # re-assign variable names for equation readability - q_fd = fire_load_density_MJm2 - HRRPUA = fire_hrr_density_MWm2 - s = fire_spread_rate_ms - h_s = beam_location_height_m - l_s = beam_location_length_m - l = room_length_m - w = room_width_m - if l < w: - l += w - w = l - w - l -= w - - # work out ventilation conditions - - # a_v = opening_height_m * opening_width_m * opening_fraction - # Qv = 1.75 * a_v * np.sqrt(opening_height_m) - - # workout burning time etc. - t_burn = max([q_fd / HRRPUA, 900.0]) - t_decay = max([t_burn, l / s]) - t_lim = min([t_burn, l / s]) - - # reduce resolution to fit time step for t_burn, t_decay, t_lim - time_interval_s = t[1] - t[0] - t_decay_ = round(t_decay / time_interval_s, 0) * time_interval_s - t_lim_ = round(t_lim / time_interval_s, 0) * time_interval_s - if t_decay_ == t_lim_: - t_lim_ -= time_interval_s - - # workout the heat release rate ARRAY (corrected with time) - Q_growth = (HRRPUA * w * s * t) * (t < t_lim_) - Q_peak = ( - min([HRRPUA * w * s * t_burn, HRRPUA * w * l]) * (t >= t_lim_) * (t <= t_decay_) - ) - Q_decay = (max(Q_peak) - (t - t_decay_) * w * s * HRRPUA) * (t > t_decay_) - Q_decay[Q_decay < 0] = 0 - Q = (Q_growth + Q_peak + Q_decay) * 1000.0 - - # workout the distance between fire median to the structural element r - l_fire_front = s * t - l_fire_front[l_fire_front < 0] = 0 - l_fire_front[l_fire_front > l] = l - l_fire_end = s * (t - t_lim) - l_fire_end[l_fire_end < 0] = 0.0 - l_fire_end[l_fire_end > l] = l - l_fire_median = (l_fire_front + l_fire_end) / 2.0 - - # workout the far field temperature of gas T_g - if isinstance(l_s, float) or isinstance(l_s, int): - r = np.absolute(l_s - l_fire_median) - T_g = np.where((r / h_s) > 0.8, (5.38 * np.power(Q / r, 2 / 3) / h_s) + 20.0, 0) - T_g = np.where( - (r / h_s) <= 0.8, - (16.9 * np.power(Q, 2 / 3) / np.power(h_s, 5 / 3)) + 20.0, - T_g, - ) - T_g[T_g >= fire_nft_limit_c] = fire_nft_limit_c - return T_g - elif isinstance(l_s, np.ndarray) or isinstance(l_s, list): - l_s_list = copy.copy(l_s) - T_g_list = list() - for l_s in l_s_list: - r = np.absolute(l_s - l_fire_median) - T_g = np.where( - (r / h_s) > 0.8, (5.38 * np.power(Q / r, 2 / 3) / h_s) + 20.0, 0 - ) - T_g = np.where( - (r / h_s) <= 0.8, - (16.9 * np.power(Q, 2 / 3) / np.power(h_s, 5 / 3)) + 20.0, - T_g, - ) - T_g[T_g >= fire_nft_limit_c] = fire_nft_limit_c - T_g_list.append(T_g) - return T_g_list - else: - raise TypeError('Unknown type of parameter "l_s": {}'.format(type(l_s))) - - -def fire_backup( - t: np.ndarray, - T_0: float, - q_fd: float, - hrrpua: float, - l: float, - w: float, - s: float, - e_h: float, - e_l: float, - T_max: float = 1323.15, -): - """ - :param t: ndarray, [s] An array representing time incorporating 'temperature'. - :param T_0: float, [K] ,Initial temperature. - :param q_fd: float, [J/m2], Fire load density. - :param hrrpua: float, [W/m2], Heat release rate density. - :param l: float, [m], Compartment length. - :param w: float, [m], Compartment width. - :param s: float, [m/s], Fire spread speed. - :param e_h: float, [m], Vertical distance between element to fuel bed. - :param e_l: float, [m], Horizontal distance between element to fire front. - :return temperature: [K] An array representing temperature incorporating 'time'. - """ - - # UNIT CONVERSION TO FIT EQUATIONS - T_0 -= 273.15 - q_fd /= 1e6 - hrrpua /= 1e6 - T_max -= 273.15 - - # workout time step - time_step = t[1] - t[0] - - # workout burning time etc. - t_burn = max([q_fd / hrrpua, 900.0]) - t_decay = max([t_burn, l / s]) - t_lim = min([t_burn, l / s]) - - # reduce resolution to fit time step for t_burn, t_decay, t_lim - t_decay_ = round(t_decay / time_step, 0) * time_step - t_lim_ = round(t_lim / time_step, 0) * time_step - if t_decay_ == t_lim_: - t_lim_ -= time_step - - # workout the heat release rate ARRAY (corrected with time) - Q_growth = (hrrpua * w * s * t) * (t < t_lim_) - Q_peak = ( - min([hrrpua * w * s * t_burn, hrrpua * w * l]) * (t >= t_lim_) * (t <= t_decay_) - ) - Q_decay = (max(Q_peak) - (t - t_decay_) * w * s * hrrpua) * (t > t_decay_) - Q_decay[Q_decay < 0] = 0 - Q = (Q_growth + Q_peak + Q_decay) * 1000.0 - - # workout the distance between fire_curve median to the structural element r - l_fire_front = s * t - l_fire_front[l_fire_front < 0] = 0.0 - l_fire_front[l_fire_front > l] = l - l_fire_end = s * (t - t_lim) - l_fire_end[l_fire_end < 0] = 0.0 - l_fire_end[l_fire_end > l] = l - l_fire_median = (l_fire_front + l_fire_end) / 2.0 - r = np.absolute(e_l - l_fire_median) - r[r == 0] = 0.001 # will cause crash if r = 0 - - # workout the far field temperature of gas T_g - T_g1 = (5.38 * np.power(Q / r, 2 / 3) / e_h) * ((r / e_h) > 0.18) - T_g2 = (16.9 * np.power(Q, 2 / 3) / np.power(e_h, 5 / 3)) * ((r / e_h) <= 0.18) - T_g = T_g1 + T_g2 + T_0 - - T_g[T_g >= T_max] = T_max - - # UNIT CONVERSION TO FIT OUTPUT (SI) - T_g = T_g + 273.15 # C -> K - Q *= 10e6 # MJ -> J - - return T_g - - -def example_plot_interflam(): - time = np.arange(0, 210 * 60, 30) - list_l = [25, 50, 100, 150] - - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel("Temperature [$℃$]") - - for length in list_l: - temperature = fire( - t=time, - fire_load_density_MJm2=600, - fire_hrr_density_MWm2=0.25, - room_length_m=length, - room_width_m=16, - fire_spread_rate_ms=0.012, - beam_location_height_m=3, - beam_location_length_m=length / 2, - fire_nft_limit_c=1050, - ) - - ax.plot(time / 60, temperature, label="Room length {:4.0f} m".format(length)) - - ax.legend(loc=4).set_visible(True) - ax.set_xlim((0, 180)) - ax.set_ylim((0, 1400)) - ax.grid(color="k", linestyle="--") - plt.tight_layout() - plt.savefig(fname="fire-travelling.png", dpi=300) - - -def _test_fire_backup(): - import numpy as np - - time = np.arange(0, 22080, 30) - list_l = [50, 100, 150] - - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel("Temperature [$℃$]") - - for l in list_l: - temperature = fire_backup( - t=time, - T_0=293.15, - q_fd=900e6, - hrrpua=0.15e6, - l=l, - w=17.4, - s=0.012, - e_h=3.5, - e_l=l / 2, - ) - ax.plot(time / 60, temperature - 273.15) - - ax.legend().set_visible(True) - ax.set_xlim((0, 120)) - ax.grid(color="k", linestyle="--") - plt.tight_layout() - # plt.show() - - -def _test_fire(): - time = np.arange(0, 210 * 60, 30) - list_l = [25, 50, 100, 150] - - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel(u"Temperature [$℃$]") - - for length in list_l: - temperature = fire( - t=time, - fire_load_density_MJm2=600, - fire_hrr_density_MWm2=0.25, - room_length_m=length, - room_width_m=16, - fire_spread_rate_ms=0.012, - beam_location_height_m=3, - beam_location_length_m=length / 2, - fire_nft_limit_c=1050, - ) - - ax.plot(time / 60, temperature, label="Room length {:4.0f} m".format(length)) - - ax.legend(loc=4).set_visible(True) - ax.set_xlim((-10, 190)) - ax.grid(color="k", linestyle="--") - plt.tight_layout() - # plt.show() - - -def _test_fire_multiple_beam_location(): - time = np.arange(0, 210 * 60, 30) - length = 100 - - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel("Temperature [$℃$]") - - temperature_list = fire( - t=time, - fire_load_density_MJm2=600, - fire_hrr_density_MWm2=0.25, - room_length_m=length, - room_width_m=16, - fire_spread_rate_ms=0.012, - beam_location_height_m=3, - beam_location_length_m=np.linspace(0, length, 12)[1:-1], - fire_nft_limit_c=1050, - ) - - for temperature in temperature_list: - ax.plot(time / 60, temperature, label="Room length {:4.0f} m".format(length)) - - ax.legend(loc=4).set_visible(True) - ax.set_xlim((-10, 190)) - ax.grid(color="k", linestyle="--") - plt.tight_layout() - # plt.show() - - -if __name__ == "__main__": - _test_fire() - _test_fire_multiple_beam_location() - # example_plot_interflam() diff --git a/sfeprapy/func/fire_travelling_flux.py b/sfeprapy/func/fire_travelling_flux.py deleted file mode 100644 index 110b822..0000000 --- a/sfeprapy/func/fire_travelling_flux.py +++ /dev/null @@ -1,207 +0,0 @@ -# -*- coding: utf-8 -*- -from typing import Union - -import numpy as np -import pandas as pd - - -def Q_star_H( - HRR: float, - height: float -): - # Calcualtes non-dimensional HRR. HRR must be in [W] and height in [m] - q_star = HRR / (1.11e06 * height **(5/2)) - return q_star - -def Q_star_D( - HRR: float, - diameter: float -): - # Calcualtes non-dimensional HRR. HRR must be in [W] and diameter in [m] - q_star = HRR / (1.11e06 * diameter **(5/2)) - return q_star - -def flame_ext_length( - q_star: float, - height: float -): - # Calculates flame extension under the ceiling according to Annex C of EN 1991-1-2 - # q_star_H is dimensionless - # height in [m] - - lh = (2.9 * height * (q_star**0.33)) - height - return lh - -def y_param( - q_star_D: float, - diameter: float, - lh: float, - radial_dist: float, - height: float, -): - # Calculates dimensionless parameter y for incident heat flux computation - # Correlations per Annex C of EN 1991-1-2 - # q_star_D is dimensionless - # remaining parameters in [m] - - # Calculate z' - - z = np.where(q_star_D < 1, 2.4 * diameter * ((q_star_D**(2/5))-(q_star_D**(2/3))),0) - z = np.where(q_star_D >= 1,2.4 * diameter * (1-(q_star_D**(2/5))),z) - - # Calculate y - - y = (radial_dist + height + z) / (lh + height + z) - - return y - -def fire( - t: np.array, - fire_load_density_MJm2: float, - fire_hrr_density_MWm2: float, - room_length_m: float, - room_width_m: float, - fire_spread_rate_ms: float, - beam_location_height_m: float, - beam_location_length_m: Union[float, list, np.ndarray], - fire_nff_limit_kW: float -): - """ - This function calculates and returns a temperature array representing travelling fire. This function is NOT in SI. - :param t: in s, is the time array - :param fire_load_density_MJm2: in MJ/m2, is the fuel density on the floor - :param fire_hrr_density_MWm2: in MW/m2, is the heat release rate density - :param room_length_m: in m, is the room length - :param room_width_m: in m, is the room width - :param fire_spread_rate_ms: in m/s, is fire spread speed - :param beam_location_height_m: in m, is the beam lateral distance to fire origin - :param beam_location_length_m: in m, is the beam height above the floor - :param fire_nff_limit_kW: in kW, is the maximum near field heat flux - :param opening_fraction: in -, is the ventilation opening proportion between 0 to 1 - :param opening_width_m: in m, is ventilation opening width - :param opening_height_m: in m, is ventilation opening height - :return q_inc: in kW, is calculated incident heat flux - """ - - # re-assign variable names for equation readability - q_fd = fire_load_density_MJm2 - HRRPUA = fire_hrr_density_MWm2 - s = fire_spread_rate_ms - h_s = beam_location_height_m - l_s = beam_location_length_m - l = room_length_m - w = room_width_m - if l < w: - l += w - w = l - w - l -= w - - # work out ventilation conditions - - # a_v = opening_height_m * opening_width_m * opening_fraction - # Qv = 1.75 * a_v * np.sqrt(opening_height_m) - - # workout burning time etc. - t_burn = max([q_fd / HRRPUA, 900.0]) - t_decay = max([t_burn, l / s]) - t_lim = min([t_burn, l / s]) - - # reduce resolution to fit time step for t_burn, t_decay, t_lim - time_interval_s = t[1] - t[0] - t_decay_ = round(t_decay / time_interval_s, 0) * time_interval_s - t_lim_ = round(t_lim / time_interval_s, 0) * time_interval_s - if t_decay_ == t_lim_: - t_lim_ -= time_interval_s - - # workout the heat release rate ARRAY (corrected with time) - Q_growth = (HRRPUA * w * s * t) * (t < t_lim_) - Q_peak = ( - min([HRRPUA * w * s * t_burn, HRRPUA * w * l]) * (t >= t_lim_) * (t <= t_decay_) - ) - Q_decay = (max(Q_peak) - (t - t_decay_) * w * s * HRRPUA) * (t > t_decay_) - Q_decay[Q_decay < 0] = 0 - Q = (Q_growth + Q_peak + Q_decay) * 1000.0 - - # workout the distance between fire median to the structural element r - l_fire_front = s * t - l_fire_front[l_fire_front < 0] = 0 - l_fire_front[l_fire_front > l] = l - l_fire_end = s * (t - t_lim) - l_fire_end[l_fire_end < 0] = 0.0 - l_fire_end[l_fire_end > l] = l - l_fire_median = (l_fire_front + l_fire_end) / 2.0 - - # fTFM parameters - fire_area = Q / (HRRPUA*1000) - fire_dia = ((fire_area / np.pi)**0.5)*2 - q_star_H = Q_star_H(Q*1000, h_s) - q_star_D = Q_star_D(Q * 1000, fire_dia) - lh = flame_ext_length(q_star_H, h_s) - lh[lh<0] = 0 - - # Distance related parameters for proximity of ceiling point to flame - r = np.absolute(l_s - l_fire_median) - y = y_param(q_star_D,fire_dia,lh,r,h_s) - - # Heat flux computation for near field - q_inc = np.where(y>0.5, 682 * np.exp(-3.4 * y),0) - q_inc = np.where(y <= 0.5, fire_nff_limit_kW, q_inc) - q_inc[q_inc > fire_nff_limit_kW] = fire_nff_limit_kW - - #Calculate far field heat flux - q_inc_ff = 682 * np.exp(-3.4 * 1) - q_inc = np.where(y > 1, q_inc_ff * (y**-3.7), q_inc) - - # Where HRR = 0 set Q_inc = 0 - q_inc = np.where(Q<=0,0,q_inc) - - return q_inc - -def _test_fire(): - time = np.arange(0, 120 * 60, 10) - list_l = [5, 10, 15, 20, 25, 30, 35] - - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_xlim(0,90) - ax.set_ylim(0,130) - ax.set_ylabel('Incident heat flux [kW/m$^2$]') - - dict_data = { - 'time': time, - } - - for beam_location_length_m in list_l: - heat_flux = fire( - t=time, - fire_load_density_MJm2=760*0.8, - fire_hrr_density_MWm2=0.5, - room_length_m=40, - room_width_m=16, - fire_spread_rate_ms=0.018, - beam_location_height_m=3, - beam_location_length_m= beam_location_length_m, - fire_nff_limit_kW=120, - ) - - dict_data['Ceiling location = '+ str(beam_location_length_m) + ' [m]'] = heat_flux - - ax.plot(time / 60, heat_flux, label="Ceiling position {:4.0f} m".format(beam_location_length_m)) - - ax.legend(loc=4).set_visible(True) - ax.grid(color="k", linestyle="--") - plt.tight_layout() - plt.show() - - df_data = pd.DataFrame.from_dict(dict_data) - - print(df_data) - - df_data.to_csv('test.csv') - -if __name__ == "__main__": - _test_fire() - diff --git a/sfeprapy/func/heat_transfer_protected_steel_ec.py b/sfeprapy/func/heat_transfer_protected_steel_ec.py deleted file mode 100644 index 9458b2e..0000000 --- a/sfeprapy/func/heat_transfer_protected_steel_ec.py +++ /dev/null @@ -1,313 +0,0 @@ -# -*- coding: utf-8 -*- -import numpy as np - - -def protected_steel_eurocode( - time, - temperature_ambient, - rho_steel, - area_steel_section, - k_protection, - rho_protection, - c_protection, - thickness_protection, - perimeter_protected, - terminate_when_cooling=False, - terminate_max_temperature=np.inf, -): - """ - SI UNITS! - This function calculate the temperature curve of protected steel section based on BS EN 1993-1-2:2005, Section 4 - . Ambient (fire) time-temperature data must be given, as well as the parameters specified below. - :param time: ndarray, [s], time evolution. - :param temperature_ambient: ndarray, [K], imposed temperature evolution. - :param rho_steel: float, [kg/m3], steel density. - :param area_steel_section: float, [m2], steel cross section area. - :param k_protection: float, [K/kg/m], protection thermal conductivity. - :param rho_protection: float, [kg/m3], protection density. - :param c_protection: float, [J/K/kg], protection thermal capacity. - :param thickness_protection: float, [m], protection thickness. - :param perimeter_protected: float, [m], protected perimeter. - :param terminate_when_cooling: bool, [-], if True then terminate and return values when first peak steel - temperature is observed. - :return temperature_steel: ndarray, [K], is calculated steel temperature. - """ - # todo: 4.2.5.2 (2) - thermal properties for the insulation material - # todo: revise BS EN 1993-1-2:2005, Clauses 4.2.5.2 - - # BS EN 1993-1-2:2005, 3.4.1.2 - # return self.__make_property("c") - - def c_steel_T(temperature): - - temperature -= 273.15 - if temperature < 20: - # warnings.warn('Temperature ({:.1f} °C) is below 20 °C'.format(temperature)) - return ( - 425 + 0.773 * 20 - 1.69e-3 * np.power(20, 2) + 2.22e-6 * np.power(20, 3) - ) - if 20 <= temperature < 600: - return ( - 425 - + 0.773 * temperature - - 1.69e-3 * np.power(temperature, 2) - + 2.22e-6 * np.power(temperature, 3) - ) - elif 600 <= temperature < 735: - return 666 + 13002 / (738 - temperature) - elif 735 <= temperature < 900: - return 545 + 17820 / (temperature - 731) - elif 900 <= temperature <= 1200: - return 650 - elif temperature > 1200: - # warnings.warn('Temperature ({:.1f} °C) is greater than 1200 °C'.format(temperature)) - return 650 - else: - # warnings.warn('Temperature ({:.1f} °C) is outside bound.'.format(temperature)) - return 0 - - V = area_steel_section - rho_a = rho_steel - lambda_p = k_protection - rho_p = rho_protection - d_p = thickness_protection - A_p = perimeter_protected - c_p = c_protection - - temperature_steel = time * 0.0 - temperature_rate_steel = time * 0.0 - specific_heat_steel = time * 0.0 - - # Check time step <= 30 seconds. [BS EN 1993-1-2:2005, Clauses 4.2.5.2 (3)] - # time_change = gradient(time) - # if np.max(time_change) > 30.: - # raise ValueError("Time step needs to be less than 30s: {0}".format(np.max(time))) - - flag_heating_started = False - - temperature_steel[0] = temperature_ambient[ - 0 - ] # initially, steel temperature is equal to ambient - temperature_ambient_ = iter(temperature_ambient) # skip the first item - next(temperature_ambient_) # skip the first item - for i, T_g in enumerate(temperature_ambient_): - i += 1 # actual index since the first item had been skipped. - try: - specific_heat_steel[i] = c_steel_T(temperature_steel[i - 1]) - except ValueError: - specific_heat_steel[i] = specific_heat_steel[i - 1] - # print(temperature_steel[i-1]) - - # Steel temperature equations are from [BS EN 1993-1-2:2005, Clauses 4.2.5.2, Eq. 4.27] - phi = (c_p * rho_p / specific_heat_steel[i] / rho_a) * d_p * A_p / V - - a = (lambda_p * A_p / V) / (d_p * specific_heat_steel[i] * rho_a) - b = (T_g - temperature_steel[i - 1]) / (1.0 + phi / 3.0) - c = (np.exp(phi / 10.0) - 1.0) * (T_g - temperature_ambient[i - 1]) - d = time[i] - time[i - 1] - - temperature_rate_steel[i] = ( - a * b * d - c - ) / d # deviated from e4.27, converted to rate [s-1] - - temperature_steel[i] = temperature_steel[i - 1] + temperature_rate_steel[i] * d - - if (temperature_rate_steel[i] > 0 and flag_heating_started is False) and time[ - i - ] > 1800: - flag_heating_started = True - - # Terminate steel temperature calculation if necessary - if ( - terminate_when_cooling - and flag_heating_started - and temperature_rate_steel[i] < 0 - ): - break - elif flag_heating_started and terminate_max_temperature < temperature_steel[i]: - break - - # NOTE: Steel temperature can be in cooling phase at the beginning of calculation, even the ambient temperature - # (fire) is hot. This is - # due to the factor 'phi' which intends to address the energy locked within the protection layer. - # The steel temperature is forced to be increased or remain as previous when ambient temperature and - # its previous temperature are all higher than the current calculated temperature. - # A better implementation is perhaps to use a 1-D heat transfer model. - - # DEPRECIATED 26 MAR 2019 - # if temperature_steel[i] < temperature_steel[i-1] or temperature_steel[i] < temperature_ambient[i]: - # temperature_rate_steel[i] = 0 - # temperature_steel[i] = temperature_steel[i-1] - - return temperature_steel - - -def protected_steel_eurocode_max_temperature( - time, - temperature_ambient, - rho_steel, - area_steel_section, - k_protection, - rho_protection, - c_protection, - thickness_protection, - perimeter_protected, - terminate_check_wait_time=3600, - terminate_max_temperature=np.inf, -): - """ - LIMITATIONS: - Constant time interval throughout - Only one maxima - - SI UNITS! - This function calculate the temperature curve of protected steel section based on BS EN 1993-1-2:2005, Section 4 - . Ambient (fire) time-temperature data must be given, as well as the parameters specified below. - - :param time: {ndarray} [s] - :param temperature_ambient: {ndarray} [K] - :param rho_steel: {float} [kg/m3] - :param area_steel_section: {float} [m2] - :param k_protection: {float} [K/kg/m] - :param rho_protection: {float} [kg/m3] - :param c_protection: {float} [J/K/kg] - :param thickness_protection: {float} [m] - :param perimeter_protected: {float} [m] - temperature is observed. - :return time: {ndarray, float} [s] - :return temperature_steel: {ndarray, float} [K] - :return data_all: {Dict} [-] - """ - # todo: 4.2.5.2 (2) - thermal properties for the insulation material - # todo: revise BS EN 1993-1-2:2005, Clauses 4.2.5.2 - - def c_steel_T(temperature): - - temperature -= 273.15 - if temperature < 20: - # warnings.warn('Temperature ({:.1f} °C) is below 20 °C'.format(temperature)) - return ( - 425 + 0.773 * 20 - 1.69e-3 * np.power(20, 2) + 2.22e-6 * np.power(20, 3) - ) - if 20 <= temperature < 600: - return ( - 425 - + 0.773 * temperature - - 1.69e-3 * np.power(temperature, 2) - + 2.22e-6 * np.power(temperature, 3) - ) - elif 600 <= temperature < 735: - return 666 + 13002 / (738 - temperature) - elif 735 <= temperature < 900: - return 545 + 17820 / (temperature - 731) - elif 900 <= temperature <= 1200: - return 650 - elif temperature > 1200: - # warnings.warn('Temperature ({:.1f} °C) is greater than 1200 °C'.format(temperature)) - return 650 - else: - # warnings.warn('Temperature ({:.1f} °C) is outside bound.'.format(temperature)) - return 0 - - V = area_steel_section - rho_a = rho_steel - lambda_p = k_protection - rho_p = rho_protection - d_p = thickness_protection - A_p = perimeter_protected - c_p = c_protection - - # temperature_steel = time * 0. - # temperature_rate_steel = time * 0. - # specific_heat_steel = time * 0. - - # Check time step <= 30 seconds. [BS EN 1993-1-2:2005, Clauses 4.2.5.2 (3)] - # time_change = gradient(time) - # if np.max(time_change) > 30.: - # raise ValueError("Time step needs to be less than 30s: {0}".format(np.max(time))) - - flag_heating_started = False - - # T_ini = temperature_ambient[0] # temperature at beginning - T = temperature_ambient[0] # current steel temperature - # T_max = -1 # maximum steel temperature - # c = c_steel_T(293.15) - d = time[1] - time[0] - - for i in range(1, len(temperature_ambient)): - - T_g = temperature_ambient[i] - - c_s = c_steel_T(T) - - # Steel temperature equations are from [BS EN 1993-1-2:2005, Clauses 4.2.5.2, Eq. 4.27] - phi = (c_p * rho_p / c_s / rho_a) * d_p * A_p / V - - a = (lambda_p * A_p / V) / (d_p * c_s * rho_a) - b = (T_g - T) / (1.0 + phi / 3.0) - c = (np.exp(phi / 10.0) - 1.0) * (T_g - temperature_ambient[i - 1]) - - dT = (a * b * d - c) / d # deviated from e4.27, converted to rate [s-1] - - T += dT * d - - if not (T >= 0 or T <= 0): - print("d") - - if not flag_heating_started: - if time[i] >= terminate_check_wait_time: - if dT > 0: - flag_heating_started = True - - # Terminate early if maximum temperature is reached - elif flag_heating_started: - if T > terminate_max_temperature: - break - if dT < 0: - T -= dT * d - break - # print(T) - return T - - -if __name__ == "__main__": - from sfeprapy.func.fire_iso834 import fire - import numpy as np - - import matplotlib.pyplot as plt - - plt.style.use("seaborn-paper") - fig, ax = plt.subplots(figsize=(3.94, 2.76)) - ax.set_xlabel("Time [minute]") - ax.set_ylabel("Temperature [$^{℃}C$]") - ax.legend().set_visible(True) - ax.grid(color="k", linestyle="--") - - rho = 7850 - t = np.arange(0, 700, 0.1) - T = fire(t, 20 + 273.15) - ax.plot(t / 60, T - 273.15, "k", label="ISO 834 fire") - - list_dp = np.arange(0.0001, 0.01 + 0.002, 0.002) - - for d_p in list_dp: - T_s = protected_steel_eurocode( - time=t, - temperature_ambient=T, - rho_steel=rho, - area_steel_section=0.017, - k_protection=0.2, - rho_protection=800, - c_protection=1700, - thickness_protection=d_p, - perimeter_protected=2.14, - ) - plt.plot( - t / 60, - T_s - 273.15, - label="Protection thickness {:2.0f} mm".format(d_p * 1000), - ) - - plt.legend(loc=4) - plt.tight_layout() - plt.show() diff --git a/sfeprapy/func/heat_transfer_unprotected_steel_ec.py b/sfeprapy/func/heat_transfer_unprotected_steel_ec.py deleted file mode 100644 index 628671a..0000000 --- a/sfeprapy/func/heat_transfer_unprotected_steel_ec.py +++ /dev/null @@ -1,67 +0,0 @@ -# -*- coding: utf-8 -*- -import numpy as np - - -def unprotected_steel_eurocode( - time, - temperature_ambient, - perimeter_section, - area_section, - perimeter_box, - density_steel, - c_steel_T, - h_conv, - emissivity_resultant, -): - """ - SI UNITS FOR INPUTS AND OUTPUTS. - :param time: - :param temperature_ambient: - :param perimeter_section: - :param area_section: - :param perimeter_box: - :param density_steel: - :param c_steel_T: - :param h_conv: - :param emissivity_resultant: - :return: - """ - - # Create steel temperature change array s - temperature_rate_steel = time * 0.0 - temperature_steel = time * 0.0 - heat_flux_net = time * 0.0 - c_s = time * 0.0 - - k_sh = ( - 0.9 * (perimeter_box / area_section) / (perimeter_section / area_section) - ) # BS EN 1993-1-2:2005 (e4.26a) - F = perimeter_section - V = area_section - rho_s = density_steel - h_c = h_conv - sigma = 56.7e-9 # todo: what is this? - epsilon = emissivity_resultant - - time_, temperature_steel[0], c_s[0] = iter(time), temperature_ambient[0], 0.0 - next(time_) - for i, v in enumerate(time_): - i += 1 - T_f, T_s_ = ( - temperature_ambient[i], - temperature_steel[i - 1], - ) # todo: steel specific heat - c_s[i] = c_steel_T(temperature_steel[i - 1] + 273.15) - - # BS EN 1993-1-2:2005 (e4.25) - a = h_c * (T_f - T_s_) - b = sigma * epsilon * (np.power(T_f, 4) - np.power(T_s_, 4)) - c = k_sh * F / V / rho_s / c_s[i] - d = time[i] - time[i - 1] - - heat_flux_net[i] = a + b - - temperature_rate_steel[i] = c * (a + b) * d - temperature_steel[i] = temperature_steel[i - 1] + temperature_rate_steel[i] - - return temperature_steel, temperature_rate_steel, heat_flux_net, c_s diff --git a/sfeprapy/func/mcs.py b/sfeprapy/func/mcs.py index 75890f2..c4e2195 100644 --- a/sfeprapy/func/mcs.py +++ b/sfeprapy/func/mcs.py @@ -1,6 +1,8 @@ import copy +import multiprocessing as mp import os import time +from abc import ABC, abstractmethod from typing import Union, Callable import pandas as pd @@ -9,7 +11,7 @@ from sfeprapy.func.mcs_gen import main as mcs_gen_main -class MCS: +class MCS(ABC): """ Monte Carlo Simulation (MCS) object defines the framework of a MCS process. MCS is designed to work as parent class and certain methods defined in this class are placeholders. @@ -31,7 +33,7 @@ class MCS: `MCS.mcs_deterministic_calc` A method to carry out deterministic calculation. NOTE! This method needs to be re-defined in a child class. - `MCS.mcs_post` + `MCS.mcs_post_per_case` A method to post processing results. NOTE! This method needs to be re-defined in a child class. """ @@ -45,7 +47,7 @@ def __init__(self): # ------------------------------ # instantiate internal variables # ------------------------------ - self.__cwd: str = None # work folder path + self.cwd: str = None # work folder path self.__mcs_inputs: dict = None # input parameters self.__mcs_config: dict = None # configuration parameters self.__mcs_sampler: Callable = mcs_gen_main # stochastic variable generator function @@ -57,17 +59,19 @@ def __init__(self): # assign default properties self.func_mcs_gen = mcs_gen_main + @abstractmethod def mcs_deterministic_calc(self, *args, **kwargs) -> dict: """Placeholder method. The Monte Carlo Simulation deterministic calculation routine. :return: """ - pass + raise NotImplementedError('This method should be overridden by a child class') + @abstractmethod def mcs_deterministic_calc_mp(self, *args, **kwargs) -> dict: """Placeholder method. The Monte Carlo Simulation deterministic calculation routine. :return: """ - pass + raise NotImplementedError('This method should be overridden by a child class') @property def mcs_inputs(self) -> dict: @@ -81,7 +85,8 @@ def mcs_inputs(self, fp_df_dict: Union[str, pd.DataFrame, dict]): """ if isinstance(fp_df_dict, str): fp = fp_df_dict - self.path_wd = os.path.dirname(fp) + if self.cwd is None: + self.cwd = os.path.dirname(fp) if fp.endswith(".xlsx") or fp.endswith(".xls"): self.__mcs_inputs = pd.read_excel(fp).set_index("PARAMETERS").to_dict() elif fp.endswith(".csv"): @@ -107,9 +112,9 @@ def mcs_config(self) -> dict: def mcs_config(self, config): self.__mcs_config = config if 'cwd' in config: - self.__cwd = config['cwd'] + self.cwd = config['cwd'] - def run_mcs(self): + def run_mcs(self, qt_prog_signal_0=None, qt_prog_signal_1=None): # ---------------------------- # Prepare mcs parameter inputs # ---------------------------- @@ -161,49 +166,58 @@ def run_mcs(self): # 'case_1': df1, # 'case_2': df2 # } + + m, p = mp.Manager(), mp.Pool(self.mcs_config["n_threads"], maxtasksperchild=1000) for k, v in x2.items(): + if qt_prog_signal_0: + qt_prog_signal_0.emit(f'{len(x3) + 1}/{len(x1)} {k}') x3_ = self.__mcs_mp( self.mcs_deterministic_calc, self.mcs_deterministic_calc_mp, x=v, n_threads=self.mcs_config["n_threads"], + m=m, + p=p, + qt_prog_signal_1=qt_prog_signal_1 ) - if self.mcs_post: - mcs_post = self.mcs_post - x3_ = mcs_post(x3_) + + # Post process output upon completion per case + if self.mcs_post_per_case: + self.mcs_post_per_case(df=x3_) x3[k] = copy.copy(x3_) - if self.__cwd: - if not os.path.exists(os.path.join(self.__cwd, self.DEFAULT_TEMP_FOLDER_NAME)): - os.makedirs(os.path.join(self.__cwd, self.DEFAULT_TEMP_FOLDER_NAME)) - x3_.to_csv( - os.path.join( - os.path.join(self.__cwd, self.DEFAULT_TEMP_FOLDER_NAME), - f"{k}.csv", - ) - ) + p.close() + p.join() # ------------ # Pack results # ------------ self.__mcs_out = pd.concat([v for v in x3.values()]) - if self.__cwd: - self.mcs_out.to_csv( - os.path.join(self.__cwd, self.DEFAULT_MCS_OUTPUT_FILE_NAME), - index=False, - ) + # Post process output upon completion of all cases + self.mcs_post_all_cases(self.__mcs_out) + + @abstractmethod + def mcs_post_per_case(self, *arg, **kwargs): + """ + This method will be called upon completion of each case + """ + raise NotImplementedError('This method should be overridden by a child class') - def mcs_post(self, *arg, **kwargs): - pass + @abstractmethod + def mcs_post_all_cases(self, *args, **kwargs): + """ + This method will be called upon completion of all cases + """ + raise NotImplementedError('This method should be overridden by a child class') @property def mcs_out(self): return self.__mcs_out @staticmethod - def __mcs_mp(func, func_mp, x: pd.DataFrame, n_threads: int) -> pd.DataFrame: + def __mcs_mp(func, func_mp, x: pd.DataFrame, n_threads: int, m, p, qt_prog_signal_1=None) -> pd.DataFrame: list_mcs_in = x.to_dict(orient="records") time.sleep(0.5) # to avoid clashes between the prints and progress bar @@ -212,31 +226,35 @@ def __mcs_mp(func, func_mp, x: pd.DataFrame, n_threads: int) -> pd.DataFrame: print("{:<24.24}: {}".format("NO. OF SIMULATIONS", len(x.index))) time.sleep(0.5) # to avoid clashes between the prints and progress bar + n_simulations = len(list_mcs_in) if n_threads == 1 or func_mp is None: mcs_out = list() + j = 0 for i in tqdm(list_mcs_in, ncols=60): mcs_out.append(func(**i)) + j += 1 + if qt_prog_signal_1: + qt_prog_signal_1.emit(int(j / n_simulations * 100)) else: - import multiprocessing as mp - m, p = mp.Manager(), mp.Pool(n_threads, maxtasksperchild=1000) q = m.Queue() jobs = p.map_async(func_mp, [(dict_, q) for dict_ in list_mcs_in]) - n_simulations = len(list_mcs_in) with tqdm(total=n_simulations, ncols=60) as pbar: while True: if jobs.ready(): if n_simulations > pbar.n: pbar.update(n_simulations - pbar.n) + if qt_prog_signal_1: + qt_prog_signal_1.emit(100) break else: if q.qsize() - pbar.n > 0: pbar.update(q.qsize() - pbar.n) + if qt_prog_signal_1: + qt_prog_signal_1.emit(int(q.qsize() / n_simulations * 100)) time.sleep(1) - p.close() - p.join() - mcs_out = jobs.get() - time.sleep(0.5) + mcs_out = jobs.get() + time.sleep(0.5) # clean and convert results to dataframe and return df_mcs_out = pd.DataFrame(mcs_out) diff --git a/sfeprapy/func/mcs_gen.py b/sfeprapy/func/mcs_gen.py index df695b7..91e86b1 100644 --- a/sfeprapy/func/mcs_gen.py +++ b/sfeprapy/func/mcs_gen.py @@ -230,7 +230,10 @@ def main(x: dict, num_samples: int) -> pd.DataFrame: ) t_ = d_.iloc[:, 0] v_ = d_.iloc[:, 1] - f_interp = interp1d(t_, v_, bounds_error=False, fill_value=0) + if all(v_ == v_[0]): + f_interp = v_[0] + else: + f_interp = interp1d(t_, v_, bounds_error=False, fill_value=0) dict_out[k] = np.full((num_samples,), f_interp) else: raise ValueError("Unknown input data type for {}.".format(k)) @@ -281,21 +284,6 @@ def _test_random_variable_generator(): assert abs(np.min(y["v"].values) - 50) <= 1 assert abs(np.mean(y["v"].values) - 420) <= 1 - # todo: lognormal distribution need to be verified. - # x = dict( - # v=dict( - # dist='lognorm_', - # ubound=2.5, - # lbound=0.00001, - # mean=0.00001, - # sd=0.25, - # ), - # ) - # y = main(x, 1000) - # assert len(y['v'].values) == 1000 - # print(np.mean(y['v'].values)) - # assert abs(np.mean(y['v'].values) - 0.5) <= 0.001 - if __name__ == "__main__": _test_random_variable_generator() diff --git a/sfeprapy/func/mcs_obj.py b/sfeprapy/func/mcs_obj.py index 9b395ef..70fbe91 100644 --- a/sfeprapy/func/mcs_obj.py +++ b/sfeprapy/func/mcs_obj.py @@ -2,7 +2,6 @@ import json import os import time -from tkinter import filedialog, Tk, StringVar from typing import Union, Callable import pandas as pd @@ -94,38 +93,30 @@ def mcs_out(self, df_out: pd.DataFrame): self._df_mcs_out = df_out def define_problem( - self, - data: Union[str, pd.DataFrame, dict] = None, - config: dict = None, - path_wd: str = None, + self, + data: Union[str, pd.DataFrame, dict], + config: dict = None, + path_wd: str = None, ): # to get problem definition: try to parse from csv/xls/xlsx - if data is None: - fp = os.path.realpath(self._get_file_path_gui()) - self.path_wd = os.path.dirname(fp) - if fp.endswith(".xlsx") or fp.endswith(".xls"): - self.input = pd.read_excel(fp).set_index("PARAMETERS").to_dict() - elif fp.endswith(".csv"): - self.input = pd.read_csv(fp).set_index("PARAMETERS").to_dict() - else: - raise ValueError("Unknown input file format.") - elif isinstance(data, str): + if isinstance(data, str): fp = data self.path_wd = os.path.dirname(fp) - if fp.endswith(".xlsx") or fp.endswith(".xls"): - self.input = pd.read_excel(fp).set_index("PARAMETERS").to_dict() - elif fp.endswith(".csv"): - self.input = pd.read_csv(fp).set_index("PARAMETERS").to_dict() - else: - raise ValueError("Unknown input file format.") + try: + self.input = pd.read_excel(fp, index_col=0).to_dict() + except Exception: + try: + self.input = pd.read_csv(fp, index_col=0).to_dict() + except Exception as e: + raise ValueError(f"Unknown input file format, {e}") elif isinstance(data, pd.DataFrame): self.input = data.to_dict() elif isinstance(data, dict): self.input = data else: - raise TypeError("Unknown input data type.") + raise TypeError("Unknown input data type, input can only be: None, a file path, DataFrame or dict") # to get configuration: try to parse from cwd if there is any, otherwise chose default values @@ -134,7 +125,7 @@ def define_problem( else: # otherwise try: with open( - os.path.join(self.path_wd, self.DEFAULT_CONFIG_FILE_NAME), "r" + os.path.join(self.path_wd, self.DEFAULT_CONFIG_FILE_NAME), "r" ) as f: self.config = json.load(f) except (TypeError, FileNotFoundError): @@ -144,7 +135,7 @@ def define_stochastic_parameter_generator(self, func): self.func_mcs_gen = func def define_calculation_routine( - self, func_mcs, func_mcs_mp=None, func_mcs_out_post=None + self, func_mcs, func_mcs_mp=None, func_mcs_out_post=None ): self.func_mcs_calc = func_mcs self.func_mcs_calc_mp = func_mcs_mp @@ -210,7 +201,7 @@ def run_mcs(self): if self.path_wd: if not os.path.exists( - os.path.join(self.path_wd, self.DEFAULT_TEMP_FOLDER_NAME) + os.path.join(self.path_wd, self.DEFAULT_TEMP_FOLDER_NAME) ): os.makedirs( os.path.join(self.path_wd, self.DEFAULT_TEMP_FOLDER_NAME) @@ -224,11 +215,12 @@ def run_mcs(self): self.mcs_out = pd.concat([v for v in x3.values()]) - if self.path_wd: - self.mcs_out.to_csv( - os.path.join(self.path_wd, self.DEFAULT_MCS_OUTPUT_FILE_NAME), - index=False, - ) + # DEPRECIATED ON 28/08/2020 + # if self.path_wd: + # self.mcs_out.to_csv( + # os.path.join(self.path_wd, self.DEFAULT_MCS_OUTPUT_FILE_NAME), + # index=False, + # ) @staticmethod def _mcs_mp(func, func_mp, x: pd.DataFrame, n_threads: int) -> pd.DataFrame: @@ -272,21 +264,3 @@ def _mcs_mp(func, func_mp, x: pd.DataFrame, n_threads: int) -> pd.DataFrame: "solver_time_equivalence_solved", inplace=True ) # sort base on time equivalence return df_mcs_out - - @staticmethod - def _get_file_path_gui(): - root = Tk() - root.withdraw() - folder_path = StringVar() - - path_input_file_csv = filedialog.askopenfile( - title="Select Input File", filetypes=[("MCS IN", [".csv", ".xlsx"])] - ) - folder_path.set(path_input_file_csv) - root.update() - - try: - path_input_file_csv = os.path.realpath(path_input_file_csv.name) - return path_input_file_csv - except AttributeError: - raise FileNotFoundError("File not found.") diff --git a/sfeprapy/func/pd_6688_1_2_2007.py b/sfeprapy/func/pd_6688_1_2_2007.py index 79d7233..1904cae 100644 --- a/sfeprapy/func/pd_6688_1_2_2007.py +++ b/sfeprapy/func/pd_6688_1_2_2007.py @@ -21,46 +21,28 @@ def annex_b_equivalent_time_of_fire_exposure( SUPPLEMENT INFO: Table A.4 - Design fire growth rates (from PD 7974-1:2003, Table 3) - ╔══════════════════════════════════════════════════════════════════╦════════════╗ - ║ BUILDING ║ USE ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Picture gallery ║ Slow ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Passenger stations and termini for air, rail, road or sea travel ║ Slow ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Classroom (school) ║ Medium ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Dwelling ║ Medium ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Office ║ Medium ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Hotel reception ║ Medium ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Hotel bedroom ║ Medium ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Hospital room ║ Medium ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Library ║ Fast ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Theatre (cinema) ║ Fast ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Shop ║ Fast ║ - ╠══════════════════════════════════════════════════════════════════╬════════════╣ - ║ Industrial storage or plant room ║ Ultra-fast ║ - ╚══════════════════════════════════════════════════════════════════╩════════════╝ + | BUILDING | USE | + |------------------------------------------------------------------|------------| + | Picture gallery | Slow | + | Passenger stations and termini for air_ rail_ road or sea travel | Slow | + | Classroom (school) | Medium | + | Dwelling | Medium | + | Office | Medium | + | Hotel reception | Medium | + | Hotel bedroom | Medium | + | Hospital room | Medium | + | Library | Fast | + | Theatre (cinema) | Fast | + | Shop | Fast | + | Industrial storage or plant room | Ultra-fast | Table A.5 - Heat release rate per unit area of fire for different occupancies (from PD 7974-1:2003) - ╔═════════════╦══════════════════════════════════════════════╗ - ║ OCCUPANCY ║ HEAT RELEASE RATE PER UNIT AREA, Q'' [kW/m2] ║ - ╠═════════════╬══════════════════════════════════════════════╣ - ║ Shops ║ 550 ║ - ╠═════════════╬══════════════════════════════════════════════╣ - ║ Offices ║ 290 ║ - ╠═════════════╬══════════════════════════════════════════════╣ - ║ Hotel rooms ║ 249 ║ - ╠═════════════╬══════════════════════════════════════════════╣ - ║ Industrial ║ 86-620 ║ - ╚═════════════╩══════════════════════════════════════════════╝ + | OCCUPANCY | HEAT RELEASE RATE PER UNIT AREA | $Q''$ [kW/m2] | | + |-------------|---------------------------------|---------------|---| + | Shops | 550 | | | + | Offices | 290 | | | + | Hotel rooms | 249 | | | + | Industrial | 86-620 | | | EXAMPLE: diff --git a/sfeprapy/func/pp.py b/sfeprapy/func/pp.py new file mode 100644 index 0000000..04454c7 --- /dev/null +++ b/sfeprapy/func/pp.py @@ -0,0 +1,182 @@ +import itertools +from typing import Tuple, Union + +import matplotlib.pyplot as plt +import numpy as np +import seaborn as sns + + +def lineplot( + x: Union[list, np.ndarray], + y: Union[list, np.ndarray], + legend_labels: Union[list, np.ndarray], + acceptable_reliability: float = None, + figsize: tuple = (3.5, 3.5), + fig=None, + ax=None, + fp_figure: str = None, + xlim: Tuple[float, float] = (15, 180), + n_legend_col: int = 1, + legend_no_overlap: bool = True, + xlabel: str = None, + ylabel: str = None, + plot_in_set=False, # @param {type:"boolean"} + plot_in_set_zoom=20, # @param {type:"number"} + plot_in_set_x=75, # @param {type:"number"} + plot_in_set_x_tol=1, # @param {type:"number"} + plot_in_set_y=0, # acceptable_reliability + plot_in_set_y_tol=0.01, # @param {type:"number"} +): + # Calculate the maximum time equivalence value (i.e. x-axis limit) + + # Generate x-axis array, based upon the calculated x-axis upper limit + # x = (edges[1:] + edges[:-1]) / 2 + + # Calculate the PDF and CDF of time equivalence, based upon the x-axis array + + # Generate plot + if fig is None and ax is None: + fig, ax = plt.subplots(1, 1, figsize=figsize, sharex=True) + + for i in range(len(y)): + sns.lineplot(x[i], y[i], label=legend_labels[i], ax=ax, palette=sns.color_palette("husl", n_colors=20)) + + if isinstance(acceptable_reliability, (float, int)): + ax.axhline(acceptable_reliability, ls='--', c='k') + ax.text(15, acceptable_reliability, f'{acceptable_reliability:.3f}', ha='left', va='bottom') + + # Update plot settings + if legend_no_overlap: + ax.legend( + shadow=False, edgecolor='k', fancybox=False, ncol=n_legend_col, fontsize='small',bbox_to_anchor=(1.04, 1), + loc="upper left" + ).set_visible(True) + else: + ax.legend(shadow=False, edgecolor='k', fancybox=False, ncol=n_legend_col, fontsize='small').set_visible(True) + ax.set_xticks(np.arange(xlim[0], xlim[1] + 1, 15)) + ax.set_xlim(*xlim) + if xlabel: + ax.set_xlabel(xlabel) + if ylabel: + ax.set_ylabel(ylabel) + + # plot inset + if plot_in_set: + from mpl_toolkits.axes_grid1.inset_locator import zoomed_inset_axes, mark_inset + + ax_ins = zoomed_inset_axes(parent_axes=ax, zoom=plot_in_set_zoom, + loc=1) # zoom-factor: 2.5, location: upper-left + + for k, v in y.items(): + pass + for i in range(len(y)): + sns.lineplot(x, y[i], ax=ax_ins, palette=sns.color_palette("husl", n_colors=20)) + ax_ins.axhline(acceptable_reliability, ls='--', c='k') + + ax_ins.set_xlim(plot_in_set_x - plot_in_set_x_tol, plot_in_set_x + plot_in_set_x_tol) + ax_ins.set_ylim(plot_in_set_y - plot_in_set_y_tol, plot_in_set_y + plot_in_set_y_tol) + ax_ins.set_xticks([plot_in_set_x, ]) + ax_ins.set_yticks([]) + # ax_ins.xaxis.set_visible(False) + ax_ins.yaxis.set_visible(False) + ax_ins.legend().set_visible(False) + + mark_inset(parent_axes=ax, inset_axes=ax_ins, loc1=2, loc2=3, ec='k', fc='w', zorder=3) + + # Save plot + if fp_figure and fig: + fig.savefig(fp_figure, dpi=300, bbox_inches='tight', transparent=True) + + return fig, ax + + +def lineplot_matrix( + dict_teq: dict, + n_cols: int = 9, + fp_figure: str = None, + xlim: Tuple[float, float] = (0, 120), + bin_width: float = 2.5, + figsize: Tuple[float, float] = (1.2, 1.2) +): + n_rows = int(np.ceil(len(dict_teq.keys()) / n_cols)) + figsize = (n_cols * figsize[0], n_rows * figsize[1]) + fig, axes1 = plt.subplots(n_rows, n_cols, figsize=figsize, sharex=True, sharey=True) + + try: + axes1 = list(itertools.chain(*axes1)) + except: + pass + + i = 0 + for k, v in dict_teq.items(): + v[v == np.inf] = np.max(v[v != np.inf]) + v[v == -np.inf] = np.min(v[v != -np.inf]) + + edges = np.arange(np.min(v) - bin_width, np.max(v) + bin_width * 2, bin_width) + centres = (edges[1:] + edges[:-1]) / 2 + hist, _ = np.histogram(v, edges, weights=np.full_like(v, 1)) + + ax1 = axes1[i] + ax2 = ax1.twinx() + + sns.distplot(v, label='PDF', ax=ax2, bins=edges, kde=False, hist=True, + hist_kws=dict(cumulative=False, density=True, color='grey')) + sns.lineplot(centres, np.cumsum(hist) / len(v), label='CDF', ax=ax1, color='k') + + ax1.set_ylim((0, 1)) + ax1.set_yticks([0, 1]) + ax1.legend().set_visible(False) + ax1.text(0.95, 0.5, k, transform=ax2.transAxes, ha='right', va='center') + ax1.tick_params(axis='y', direction='in') + + ax2.set_xlabel('') + ax2.set_yticks([]) + ax2.set_ylim(0, 0.2) + ax2.legend().set_visible(False) + ax2.tick_params(axis='both', direction='in') + + i = i + 1 + + ax1.set_xticks(np.arange(xlim[0], xlim[1] + 1, 30)) + ax1.set_xlim(*xlim) + fig.text(-0.01, 0.5, 'CDF [-]', ha='center', va='top', rotation=90) + fig.text(0.5, 0, 'Equivalent of time exposure [$min$]', ha='center', va='top') + + fig.tight_layout(pad=0.1) + + if fp_figure: + fig.savefig(fp_figure, dpi=300, bbox_inches='tight') + + return fig, axes1 + + +def influence_factors(): + # calcualte CCDF fractile for each curve + fire_rating_target = 120 # @param {type:"number"} + list_teq_for_fr_target = list() + for i, case_name_ in enumerate(list_case_name): + teq_cdf = dict_teq_cdf[case_name_] + teq_cdf_120 = np.max(teq_cdf[x <= fire_rating_target]) + list_teq_for_fr_target.append(teq_cdf_120) + + list_teq_for_fr_target = np.asarray(list_teq_for_fr_target) + print(list_teq_for_fr_target) + + # build all data into a dict before saving + # influence_factor = (np.max(list_teq_for_fr_target)-list_teq_for_fr_target) * np.array(list_weight) + influence_factor = (1 - list_teq_for_fr_target) * np.array(list_weight) + influence_factor = influence_factor / np.linalg.norm(influence_factor) + + dict_data_influence_factor = { + 'case name': list_case_name, + 'probability weight': list_weight, + f'fractile at {fire_rating_target:.0f} minutes fire resistance rating': list_teq_for_fr_target, + 'influence factor': influence_factor + } + + print(influence_factor) + + df_influence = pd.DataFrame(data=dict_data_influence_factor) + + # save influence factor numerical data + df_influence.to_csv(os.path.join(path_work_directory, 'mcs.out.influence.csv')) diff --git a/sfeprapy/gui/__init__.py b/sfeprapy/gui/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/sfeprapy/gui/__main__.py b/sfeprapy/gui/__main__.py deleted file mode 100644 index 8106faf..0000000 --- a/sfeprapy/gui/__main__.py +++ /dev/null @@ -1,22 +0,0 @@ - -import sys - -from PySide2 import QtWidgets - -from sfeprapy.gui.logic.main import MainWindow - - -def main(app_window=None): - """main function to fire up the application cycle""" - - app = QtWidgets.QApplication(sys.argv) - if app_window is None: - window = MainWindow() - else: - window = app_window() - window.show() - sys.exit(app.exec_()) - - -if __name__ == "__main__": - main() diff --git a/sfeprapy/gui/layout/__init__.py b/sfeprapy/gui/layout/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/sfeprapy/gui/layout/main.py b/sfeprapy/gui/layout/main.py deleted file mode 100644 index 7bff0e0..0000000 --- a/sfeprapy/gui/layout/main.py +++ /dev/null @@ -1,90 +0,0 @@ -# -*- coding: utf-8 -*- - -################################################################################ -## Form generated from reading UI file 'main.ui' -## -## Created by: Qt User Interface Compiler version 5.14.1 -## -## WARNING! All changes made in this file will be lost when recompiling UI file! -################################################################################ - -from PySide2.QtCore import (QCoreApplication, QMetaObject, QObject, QPoint, - QRect, QSize, QUrl, Qt) -from PySide2.QtGui import (QBrush, QColor, QConicalGradient, QCursor, QFont, - QFontDatabase, QIcon, QLinearGradient, QPalette, QPainter, QPixmap, - QRadialGradient) -from PySide2.QtWidgets import * - - -class Ui_MainWindow(object): - def setupUi(self, MainWindow): - if MainWindow.objectName(): - MainWindow.setObjectName(u"MainWindow") - MainWindow.resize(446, 186) - sizePolicy = QSizePolicy(QSizePolicy.Minimum, QSizePolicy.Minimum) - sizePolicy.setHorizontalStretch(0) - sizePolicy.setVerticalStretch(0) - sizePolicy.setHeightForWidth(MainWindow.sizePolicy().hasHeightForWidth()) - MainWindow.setSizePolicy(sizePolicy) - MainWindow.setMinimumSize(QSize(426, 178)) - self.centralwidget = QWidget(MainWindow) - self.centralwidget.setObjectName(u"centralwidget") - self.pushButton_run = QPushButton(self.centralwidget) - self.pushButton_run.setObjectName(u"pushButton_run") - self.pushButton_run.setGeometry(QRect(335, 130, 96, 26)) - font = QFont() - font.setPointSize(9) - self.pushButton_run.setFont(font) - self.pushButton_test = QPushButton(self.centralwidget) - self.pushButton_test.setObjectName(u"pushButton_test") - self.pushButton_test.setGeometry(QRect(230, 130, 96, 26)) - self.pushButton_test.setFont(font) - self.groupBox = QGroupBox(self.centralwidget) - self.groupBox.setObjectName(u"groupBox") - self.groupBox.setGeometry(QRect(15, 15, 416, 101)) - self.groupBox.setFont(font) - self.pushButton_input_file = QPushButton(self.groupBox) - self.pushButton_input_file.setObjectName(u"pushButton_input_file") - self.pushButton_input_file.setGeometry(QRect(300, 25, 96, 26)) - self.label = QLabel(self.groupBox) - self.label.setObjectName(u"label") - self.label.setGeometry(QRect(10, 25, 101, 26)) - self.lineEdit_input_file = QLineEdit(self.groupBox) - self.lineEdit_input_file.setObjectName(u"lineEdit_input_file") - self.lineEdit_input_file.setGeometry(QRect(132, 25, 156, 26)) - self.lineEdit_n_proc = QLineEdit(self.groupBox) - self.lineEdit_n_proc.setObjectName(u"lineEdit_n_proc") - self.lineEdit_n_proc.setGeometry(QRect(132, 55, 156, 26)) - self.label_2 = QLabel(self.groupBox) - self.label_2.setObjectName(u"label_2") - self.label_2.setGeometry(QRect(10, 55, 101, 26)) - MainWindow.setCentralWidget(self.centralwidget) - self.menubar = QMenuBar(MainWindow) - self.menubar.setObjectName(u"menubar") - self.menubar.setGeometry(QRect(0, 0, 446, 22)) - MainWindow.setMenuBar(self.menubar) - self.statusbar = QStatusBar(MainWindow) - self.statusbar.setObjectName(u"statusbar") - MainWindow.setStatusBar(self.statusbar) - - self.retranslateUi(MainWindow) - - QMetaObject.connectSlotsByName(MainWindow) - # setupUi - - def retranslateUi(self, MainWindow): - MainWindow.setWindowTitle(QCoreApplication.translate("MainWindow", u"MainWindow", None)) -#if QT_CONFIG(statustip) - self.pushButton_run.setStatusTip(QCoreApplication.translate("MainWindow", u"Run simulation", None)) -#endif // QT_CONFIG(statustip) - self.pushButton_run.setText(QCoreApplication.translate("MainWindow", u"Run", None)) -#if QT_CONFIG(statustip) - self.pushButton_test.setStatusTip(QCoreApplication.translate("MainWindow", u"Make and save a template input file", None)) -#endif // QT_CONFIG(statustip) - self.pushButton_test.setText(QCoreApplication.translate("MainWindow", u"Make a Temp.", None)) - self.groupBox.setTitle(QCoreApplication.translate("MainWindow", u"Inputs", None)) - self.pushButton_input_file.setText(QCoreApplication.translate("MainWindow", u"Select", None)) - self.label.setText(QCoreApplication.translate("MainWindow", u"Input File", None)) - self.label_2.setText(QCoreApplication.translate("MainWindow", u"No. of Processes", None)) - # retranslateUi - diff --git a/sfeprapy/gui/layout/ui/main.ui b/sfeprapy/gui/layout/ui/main.ui deleted file mode 100644 index 0d618a9..0000000 --- a/sfeprapy/gui/layout/ui/main.ui +++ /dev/null @@ -1,163 +0,0 @@ - - - MainWindow - - - - 0 - 0 - 446 - 186 - - - - - 0 - 0 - - - - - 426 - 178 - - - - MainWindow - - - - - - 335 - 130 - 96 - 26 - - - - - 9 - - - - Run simulation - - - Run - - - - - - 230 - 130 - 96 - 26 - - - - - 9 - - - - Make and save a template input file - - - Make a Temp. - - - - - - 15 - 15 - 416 - 101 - - - - - 9 - - - - Inputs - - - - - 300 - 25 - 96 - 26 - - - - Select - - - - - - 10 - 25 - 101 - 26 - - - - Input File - - - - - - 132 - 25 - 156 - 26 - - - - - - - 132 - 55 - 156 - 26 - - - - - - - 10 - 55 - 101 - 26 - - - - No. of Processes - - - - - - - - 0 - 0 - 446 - 22 - - - - - - - - diff --git a/sfeprapy/gui/layout/ui2py.py b/sfeprapy/gui/layout/ui2py.py deleted file mode 100644 index 6dff975..0000000 --- a/sfeprapy/gui/layout/ui2py.py +++ /dev/null @@ -1,21 +0,0 @@ - - -def ui2py(): - import os - import subprocess - - list_ui_file_names = [ - 'main.ui', - ] - - cwd = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'ui') - destination_dir = os.path.dirname(os.path.realpath(__file__)) - - for ui_file_name in list_ui_file_names: - cmd = f'pyside2-uic {os.path.join(cwd, ui_file_name)} > {os.path.join(destination_dir, ui_file_name.replace(".ui", ".py"))}' - print(cmd) - subprocess.call(cmd, shell=True) - - -if __name__ == '__main__': - ui2py() diff --git a/sfeprapy/gui/logic/__init__.py b/sfeprapy/gui/logic/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/sfeprapy/gui/logic/main.py b/sfeprapy/gui/logic/main.py deleted file mode 100644 index b890648..0000000 --- a/sfeprapy/gui/logic/main.py +++ /dev/null @@ -1,103 +0,0 @@ -import sys - -from PySide2 import QtWidgets, QtGui, QtCore - -from sfeprapy.gui.layout.main import Ui_MainWindow -from sfeprapy.mcs0 import EXAMPLE_INPUT_CSV -from sfeprapy.mcs0.__main__ import main as mcs0 - - -class MainWindow(QtWidgets.QMainWindow): - - def __init__(self): - super().__init__() # inherit base class - - self.ui = Ui_MainWindow() # instantiate ui - self.ui.setupUi(self) # set up ui - - self.setWindowTitle('SfePrapy TEQ MCS Alpha') - - # set ui elements initial status - self.ui.lineEdit_input_file.setReadOnly(True) - self.ui.lineEdit_input_file.setText('Click the Select button.') - - # validators - self.ui.lineEdit_n_proc.setValidator(QtGui.QRegExpValidator(QtCore.QRegExp(r'^[0-9]+$'))) - - # set ui elements connections - self.ui.pushButton_input_file.clicked.connect(self.select_input_file) - self.ui.pushButton_test.clicked.connect(self.make_input_template) - self.ui.pushButton_run.clicked.connect(self.run) - - def select_input_file(self): - """select input file and copy its path to ui object""" - - # dialog to select file - path_to_file, _ = QtWidgets.QFileDialog.getOpenFileName( - self, - "Select File", - "~/", - "Spreadsheet (*.csv *.xls *.xlsx)") - - # paste the select file path to the ui object - self.ui.lineEdit_input_file.setText(path_to_file) - - def make_input_template(self): - """save an input template file to a user specified location""" - - # get the user desired file path - path_to_file, _ = QtWidgets.QFileDialog.getSaveFileName( - parent=self, - caption='Save Template Input File', - dir='mcs0.csv' - ) - - # save file - with open(path_to_file, "w+") as f: - f.write(EXAMPLE_INPUT_CSV) - - def run(self): - """to fetch input and run simulation""" - - # get input file path - try: - fp_mcs_in = self.ui.lineEdit_input_file.text() - if len(fp_mcs_in) == 0: # raise error if path length is zero. - raise - except: - self.ui.statusbar.showMessage('Run failed, unable to fetch input file path.') - self.repaint() - return -1 - - # get number of processes - try: - n_threads = int(self.ui.lineEdit_n_proc.text()) - except ValueError: - self.ui.statusbar.showMessage('Run failed, unable to fetch number of processes.') - self.repaint() - return -1 - - # run simulation - mcs0(fp_mcs_in=fp_mcs_in, n_threads=int(n_threads)) - - def update_status_bar_text(self, text: str): - self.ui.statusbar.showMessage(text=text) - - -def main(app_window=None): - """main function to fire up the application cycle""" - - app = QtWidgets.QApplication(sys.argv) - if app_window is None: - window = MainWindow() - else: - window = app_window() - window.show() - sys.exit(app.exec_()) - - -if __name__ == "__main__": - """update ui files and run the application. this is used only for testing""" - from sfeprapy.gui.layout.ui2py import ui2py - ui2py() - main() diff --git a/sfeprapy/mcs0/__init__.py b/sfeprapy/mcs0/__init__.py index 0f875f9..e6efc3e 100644 --- a/sfeprapy/mcs0/__init__.py +++ b/sfeprapy/mcs0/__init__.py @@ -12,27 +12,18 @@ def __example_input_dict(): y = { "Standard Case 1": dict( case_name="Standard Case 1", - n_simulations=500, - probability_weight=1 / 3, - fire_time_step=30, + n_simulations=1000, + fire_time_step=10, fire_time_duration=18000, - fire_hrr_density=dict( - dist="uniform_", lbound=0.25 - 0.001, ubound=0.25 + 0.001 - ), - fire_load_density=dict( - dist="gumbel_r_", lbound=10, ubound=1500, mean=420, sd=126 - ), + fire_hrr_density=dict(dist="uniform_", lbound=0.25 - 0.001, ubound=0.25 + 0.001), + fire_load_density=dict(dist="gumbel_r_", lbound=10, ubound=1500, mean=420, sd=126), fire_spread_speed=dict(dist="uniform_", lbound=0.0035, ubound=0.0190), - fire_nft_limit=dict( - dist="norm_", lbound=623.15, ubound=2023.15, mean=1323.15, sd=93 - ), + fire_nft_limit=dict(dist="norm_", lbound=623.15, ubound=2023.15, mean=1323.15, sd=93), fire_combustion_efficiency=dict(dist="uniform_", lbound=0.8, ubound=1.0), - window_open_fraction=dict( - dist="lognorm_mod_", ubound=0.9999, lbound=0.0001, mean=0.2, sd=0.2 - ), + window_open_fraction=dict(dist="lognorm_mod_", ubound=0.9999, lbound=0.0001, mean=0.2, sd=0.2), phi_teq=dict(dist="constant_", ubound=1, lbound=1, mean=0, sd=0), beam_cross_section_area=0.017, - beam_position_horizontal=-1, + beam_position_horizontal=dict(dist="uniform_", lbound=0.6 * 31.25, ubound=0.9 * 31.25), beam_position_vertical=3.2, beam_rho=7850, fire_mode=3, @@ -50,7 +41,7 @@ def __example_input_dict(): solver_temperature_goal=893.15, solver_max_iter=20, solver_thickness_lbound=0.0001, - solver_thickness_ubound=0.0500, + solver_thickness_ubound=0.0300, solver_tol=1.0, window_height=2.8, window_width=72, @@ -61,30 +52,26 @@ def __example_input_dict(): timber_density=400, # [kg/m3] timber_solver_ilim=20, timber_solver_tol=1, + p1=3e-7, + p2=0.1, + p3=0.25, + p4=0.09, + general_room_floor_area=500, ), "Standard Case 2 (with teq_phi)": dict( case_name="Standard Case 2 (with teq_phi)", - n_simulations=500, - probability_weight=1 / 3, - fire_time_step=30, + n_simulations=1000, + fire_time_step=10, fire_time_duration=18000, - fire_hrr_density=dict( - dist="uniform_", lbound=0.25 - 0.001, ubound=0.25 + 0.001 - ), - fire_load_density=dict( - dist="gumbel_r_", lbound=10, ubound=1500, mean=420, sd=126 - ), + fire_hrr_density=dict(dist="uniform_", lbound=0.25 - 0.001, ubound=0.25 + 0.001), + fire_load_density=dict(dist="gumbel_r_", lbound=10, ubound=1500, mean=420, sd=126), fire_spread_speed=dict(dist="uniform_", lbound=0.0035, ubound=0.0190), - fire_nft_limit=dict( - dist="norm_", lbound=623.15, ubound=2023.15, mean=1323.15, sd=93 - ), + fire_nft_limit=dict(dist="norm_", lbound=623.15, ubound=2023.15, mean=1323.15, sd=93), fire_combustion_efficiency=dict(dist="uniform_", lbound=0.8, ubound=1.0), - window_open_fraction=dict( - dist="lognorm_mod_", ubound=0.9999, lbound=0.0001, mean=0.2, sd=0.2 - ), + window_open_fraction=dict(dist="lognorm_mod_", ubound=0.9999, lbound=0.0001, mean=0.2, sd=0.2), phi_teq=dict(dist="lognorm_", ubound=3, lbound=0.00001, mean=1, sd=0.25), beam_cross_section_area=0.017, - beam_position_horizontal=-1, + beam_position_horizontal=dict(dist="uniform_", lbound=0.6 * 31.25, ubound=0.9 * 31.25), beam_position_vertical=3.2, beam_rho=7850, fire_mode=3, @@ -102,7 +89,7 @@ def __example_input_dict(): solver_temperature_goal=893.15, solver_max_iter=20, solver_thickness_lbound=0.0001, - solver_thickness_ubound=0.0500, + solver_thickness_ubound=0.0300, solver_tol=1.0, window_height=2.8, window_width=72, @@ -113,30 +100,26 @@ def __example_input_dict(): timber_density=400, # [kg/m3] timber_solver_ilim=20, timber_solver_tol=1, + p1=3e-7, + p2=0.1, + p3=0.25, + p4=0.09, + general_room_floor_area=500, ), "Standard Case 3 (with timber)": dict( case_name="Standard Case 3 (with timber)", - n_simulations=500, - probability_weight=1 / 3, - fire_time_step=30, + n_simulations=1000, + fire_time_step=10, fire_time_duration=18000, - fire_hrr_density=dict( - dist="uniform_", lbound=0.25 - 0.001, ubound=0.25 + 0.001 - ), - fire_load_density=dict( - dist="gumbel_r_", lbound=10, ubound=1500, mean=420, sd=126 - ), + fire_hrr_density=dict(dist="uniform_", lbound=0.25 - 0.001, ubound=0.25 + 0.001), + fire_load_density=dict(dist="gumbel_r_", lbound=10, ubound=1500, mean=420, sd=126), fire_spread_speed=dict(dist="uniform_", lbound=0.0035, ubound=0.0190), - fire_nft_limit=dict( - dist="norm_", lbound=623.15, ubound=2023.15, mean=1323.15, sd=93 - ), + fire_nft_limit=dict(dist="norm_", lbound=623.15, ubound=2023.15, mean=1323.15, sd=93), fire_combustion_efficiency=dict(dist="uniform_", lbound=0.8, ubound=1.0), - window_open_fraction=dict( - dist="lognorm_mod_", ubound=0.9999, lbound=0.0001, mean=0.2, sd=0.2 - ), + window_open_fraction=dict(dist="lognorm_mod_", ubound=0.9999, lbound=0.0001, mean=0.2, sd=0.2), phi_teq=dict(dist="constant_", ubound=1, lbound=1, mean=0, sd=0), beam_cross_section_area=0.017, - beam_position_horizontal=-1, + beam_position_horizontal=dict(dist="uniform_", lbound=0.6 * 31.25, ubound=0.9 * 31.25), beam_position_vertical=3.2, beam_rho=7850, fire_mode=3, @@ -154,7 +137,7 @@ def __example_input_dict(): solver_temperature_goal=893.15, solver_max_iter=20, solver_thickness_lbound=0.0001, - solver_thickness_ubound=0.0500, + solver_thickness_ubound=0.0300, solver_tol=1.0, window_height=2.8, window_width=72, @@ -165,6 +148,11 @@ def __example_input_dict(): timber_density=400, # [kg/m3] timber_solver_ilim=20, timber_solver_tol=1, + p1=3e-7, + p2=0.1, + p3=0.25, + p4=0.09, + general_room_floor_area=500, ), } return y @@ -190,6 +178,7 @@ def __example_input_df(x: dict) -> pd.DataFrame: EXAMPLE_INPUT_CSV = __example_input_csv(__example_input_dict()) EXAMPLE_INPUT_DF = __example_input_df(__example_input_dict()) + if __name__ == "__main__": print(EXAMPLE_CONFIG_DICT, "\n") print(EXAMPLE_INPUT_DICT, "\n") diff --git a/sfeprapy/mcs0/__main__.py b/sfeprapy/mcs0/__main__.py index 5839f0b..83b2a6d 100644 --- a/sfeprapy/mcs0/__main__.py +++ b/sfeprapy/mcs0/__main__.py @@ -1,172 +1,26 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- +import os +import warnings +from sfeprapy.mcs0.mcs0_calc import MCS0 -def save_figure(fp_mcs0_out: str): - import os - import numpy as np - import pandas as pd - import plotly.io - import plotly.graph_objects as go - - mcs_out = pd.read_csv(fp_mcs0_out) - - # A helper function to produce (x, y) line plot based upon samples - def cdf_xy(v, xlim, bin_width=0.1, weights=None): - edges = np.arange(*xlim, bin_width) - x = (edges[1:] + edges[:-1]) / 2 - y_pdf = np.histogram(v, bins=edges, weights=weights)[0] / len(v) - y_cdf = np.cumsum(y_pdf) - return x, y_cdf - - # Settings - # fig_size = np.array([3.5, 3.5]) * 1 # in inch - fig_x_limit = (0, 180) - fig_y_limit = (0, 1) - bin_width = 0.5 - - # Process time equivalence value, obtain `probability_weight` and `n_simulations` - list_case_name = sorted(list(set(mcs_out["case_name"].values))) - list_t_eq = list() - list_weight = list() - list_n_simulation = list() - for case_name in list_case_name: - teq = np.asarray( - mcs_out[mcs_out["case_name"] == case_name][ - "solver_time_equivalence_solved" - ].values, - float, - ) - teq[teq == np.inf] = np.max(teq[teq != np.inf]) - teq[teq == -np.inf] = np.min(teq[teq != -np.inf]) - teq = teq[~np.isnan(teq)] - teq[ - teq > 18000.0 - ] = 18000.0 # limit maximum time equivalence plot value to 5 hours - list_t_eq.append(teq / 60.0) - list_weight.append( - np.average( - mcs_out[mcs_out["case_name"] == case_name]["probability_weight"].values - ) - ) - list_n_simulation.append( - np.average( - mcs_out[mcs_out["case_name"] == case_name]["n_simulations"].values - ) - ) - - # Time equivalence samples -> x, y of cumulative density function - xlim = (0, np.max([np.max(v) for v in list_t_eq]) + bin_width) - x = np.arange(*xlim, bin_width) - list_cdf_t_eq_y = [ - cdf_xy(t_eq, xlim, bin_width, weights=None)[1] for t_eq in list_t_eq - ] - - # Combined time equivalence cdf - cdf_t_eq_y_combined = np.array( - [list_weight[i] * v for i, v in enumerate(list_cdf_t_eq_y)] - ) - cdf_t_eq_y_combined = np.sum(cdf_t_eq_y_combined, axis=0) - - # Plot figure - fig = go.Figure() - - # add combined time equivalence to the plot - if len(list_case_name) > 1 and np.max(cdf_t_eq_y_combined) > 0.9: - fig.add_trace( - go.Scatter(x=x, y=cdf_t_eq_y_combined, mode="lines", name="Combined") - ) - - # add individual time equivalence to the plot - for i, t_eq in enumerate(list_cdf_t_eq_y): - fig.add_trace(go.Scatter(x=x, y=t_eq, mode="lines", name=list_case_name[i])) - - fig.update_layout( - autosize=True, - paper_bgcolor="White", - plot_bgcolor="White", - xaxis=dict( - title="Equivalent of time exposure [minute]", - dtick=30, - range=fig_x_limit, - showline=True, - linewidth=1, - linecolor="black", - mirror=True, - visible=True, - showgrid=True, - gridcolor="Black", - gridwidth=1, - ticks="outside", - zeroline=False, - ), - yaxis=dict( - title="CDF", - dtick=0.2, - range=fig_y_limit, - showline=True, - linewidth=1, - linecolor="black", - mirror=True, - visible=True, - showgrid=True, - gridcolor="Black", - gridwidth=1, - ticks="outside", - zeroline=False, - ), - legend=dict( - x=0.98, - xanchor="right", - y=0.02, - yanchor="bottom", - traceorder="normal", - font=dict(family="sans-serif", color="black"), - bgcolor="White", - bordercolor="Black", - borderwidth=1, - ), - ) - - config = { - "scrollZoom": False, - "displayModeBar": True, - "editable": True, - "showLink": False, - "displaylogo": False, - } - - plotly.io.write_html( - fig, - file=os.path.join(os.path.dirname(fp_mcs0_out), "mcs.out.html"), - auto_open=True, - config=config, - ) +warnings.filterwarnings("ignore") def main(fp_mcs_in: str, n_threads: int = None): - - import os - import warnings - from sfeprapy.func.mcs_obj import MCS - from sfeprapy.mcs0.mcs0_calc import teq_main, teq_main_wrapper, mcs_out_post - from sfeprapy.func.mcs_gen import main as gen - - warnings.filterwarnings("ignore") - fp_mcs_in = os.path.realpath(fp_mcs_in) - mcs = MCS() - mcs.define_problem(fp_mcs_in) - mcs.define_stochastic_parameter_generator(gen) - mcs.define_calculation_routine(teq_main, teq_main_wrapper, mcs_out_post) + mcs = MCS0() + try: + mcs.mcs_inputs = fp_mcs_in + except Exception as e: + raise e try: if n_threads: - mcs.config = ( + mcs.mcs_config = ( dict(n_threads=n_threads) - if mcs.config - else mcs.config.update(dict(n_threads=n_threads)) + if mcs.mcs_config + else mcs.mcs_config.update(dict(n_threads=n_threads)) ) except KeyError: pass diff --git a/sfeprapy/mcs0/mcs0_calc.py b/sfeprapy/mcs0/mcs0_calc.py index 457e226..f8a424e 100644 --- a/sfeprapy/mcs0/mcs0_calc.py +++ b/sfeprapy/mcs0/mcs0_calc.py @@ -1,47 +1,47 @@ # -*- coding: utf-8 -*- import copy +import os +import threading import warnings -from typing import Union +from typing import Union, Callable import numpy as np import pandas as pd +from fsetools.lib.fse_bs_en_1991_1_2_parametric_fire import temperature as _fire_param +from fsetools.lib.fse_bs_en_1993_1_2_heat_transfer_c import protection_thickness as _protection_thickness +from fsetools.lib.fse_bs_en_1993_1_2_heat_transfer_c import temperature as _steel_temperature +from fsetools.lib.fse_bs_en_1993_1_2_heat_transfer_c import temperature_max as _steel_temperature_max +from fsetools.lib.fse_din_en_1991_1_2_parametric_fire import temperature as _fire_param_ger +from fsetools.lib.fse_travelling_fire import temperature as fire_travelling from scipy.interpolate import interp1d from sfeprapy.func.asciiplot import AsciiPlot -from sfeprapy.func.fire_parametric_ec import fire as _fire_param -from sfeprapy.func.fire_parametric_ec_din import fire as _fire_param_ger -from sfeprapy.func.fire_travelling import fire as fire_travelling -from sfeprapy.func.heat_transfer_protected_steel_ec import ( - protected_steel_eurocode as _steel_temperature, -) -from sfeprapy.func.heat_transfer_protected_steel_ec import ( - protected_steel_eurocode_max_temperature as _steel_temperature_max, -) from sfeprapy.func.mcs import MCS def _fire_travelling(**kwargs): if isinstance(kwargs["beam_location_length_m"], list) or isinstance( - kwargs["beam_location_length_m"], np.ndarray + kwargs["beam_location_length_m"], np.ndarray ): kwarg_ht_ec = dict( - time=kwargs["t"], - rho_steel=7850, - area_steel_section=0.017, - k_protection=0.2, - rho_protection=800, - c_protection=1700, - thickness_protection=0.005, - perimeter_protected=2.14, + fire_time=kwargs["t"], + beam_rho=7850, + beam_cross_section_area=0.017, + protection_k=0.2, + protection_rho=800, + protection_c=1700, + protection_thickness=0.005, + protection_protected_perimeter=2.14, ) temperature_steel_list = list() temperature_gas_list = fire_travelling(**kwargs) for temperature in temperature_gas_list: - kwarg_ht_ec["temperature_ambient"] = temperature + 273.15 - temperature_steel_list.append(_steel_temperature_max(**kwarg_ht_ec)) + kwarg_ht_ec["fire_temperature"] = temperature + 273.15 + T_a_max, t = _steel_temperature_max(**kwarg_ht_ec) + temperature_steel_list.append(T_a_max) return ( temperature_gas_list[np.argmax(temperature_steel_list)] + 273.15, @@ -49,25 +49,26 @@ def _fire_travelling(**kwargs): ) elif isinstance(kwargs["beam_location_length_m"], float) or isinstance( - kwargs["beam_location_length_m"], int + kwargs["beam_location_length_m"], int ): return fire_travelling(**kwargs) + 273.15, kwargs["beam_location_length_m"] def decide_fire( - window_height: float, - window_width: float, - window_open_fraction: float, - room_breadth: float, - room_depth: float, - room_height: float, - fire_mode: int, - fire_load_density: float, - fire_combustion_efficiency: float, - fire_hrr_density: float, - fire_spread_speed: float, - **_, + window_height: float, + window_width: float, + window_open_fraction: float, + room_breadth: float, + room_depth: float, + room_height: float, + fire_mode: int, + fire_load_density: float, + fire_combustion_efficiency: float, + fire_hrr_density: float, + fire_spread_speed: float, + *_, + **__, ) -> dict: """Calculates equivalent time exposure for a protected steel element member in more realistic fire environment opposing to the standard fire curve ISO 834. @@ -99,13 +100,11 @@ def decide_fire( room_floor_area = room_breadth * room_depth # Room internal surface area, total, including window openings - room_total_area = (2 * room_floor_area) + ( - (room_breadth + room_depth) * 2 * room_height - ) + room_total_area = (2 * room_floor_area) + ((room_breadth + room_depth) * 2 * room_height) # Fire load density related to the total surface area A_t fire_load_density_total = ( - fire_load_density_deducted * room_floor_area / room_total_area + fire_load_density_deducted * room_floor_area / room_total_area ) # Opening factor @@ -115,15 +114,13 @@ def decide_fire( fire_spread_entire_room_time = room_depth / fire_spread_speed burn_out_time = max([fire_load_density_deducted / fire_hrr_density, 900.0]) - if ( - fire_mode == 0 or fire_mode == 1 or fire_mode == 2 - ): # enforced to selected fire, i.e. 0 is ec parametric; 1 is travelling; and 2 is din ec parametric + if fire_mode == 0 or fire_mode == 1 or fire_mode == 2: # enforced to selected fire, i.e. 0 is ec parametric; 1 is travelling; and 2 is din ec parametric fire_type = fire_mode elif fire_mode == 3: # enforced to ec parametric + travelling if ( - fire_spread_entire_room_time < burn_out_time - and 0.01 < opening_factor <= 0.2 - and 50 <= fire_load_density_total <= 1000 + fire_spread_entire_room_time < burn_out_time + and 0.01 < opening_factor <= 0.2 + and 50 <= fire_load_density_total <= 1000 ): fire_type = 0 # parametric fire else: # Otherwise, it is a travelling fire @@ -131,8 +128,8 @@ def decide_fire( elif fire_mode == 4: # enforced to german parametric + travelling # If fire spreads throughout compartment and ventilation is within EC limits = Parametric fire if ( - fire_spread_entire_room_time < burn_out_time - and 0.125 <= (window_area / room_floor_area) <= 0.5 + fire_spread_entire_room_time < burn_out_time + and 0.125 <= (window_area / room_floor_area) <= 0.5 ): fire_type = 2 # german parametric else: # Otherwise, it is a travelling fire @@ -140,32 +137,31 @@ def decide_fire( else: raise ValueError("Unknown fire mode {fire_mode}.".format(fire_mode=fire_mode)) - results = dict(fire_type=fire_type) - - return results + return dict(fire_type=fire_type) def evaluate_fire_temperature( - window_height: float, - window_width: float, - window_open_fraction: float, - room_breadth: float, - room_depth: float, - room_height: float, - room_wall_thermal_inertia: float, - fire_tlim: float, - fire_type: float, - fire_time: Union[list, np.ndarray], - fire_nft_limit: float, - fire_load_density: float, - fire_combustion_efficiency: float, - fire_hrr_density: float, - fire_spread_speed: float, - fire_t_alpha: float, - fire_gamma_fi_q: float, - beam_position_vertical: float, - beam_position_horizontal: Union[np.ndarray, list, float] = -1.0, - **_, + window_height: float, + window_width: float, + window_open_fraction: float, + room_breadth: float, + room_depth: float, + room_height: float, + room_wall_thermal_inertia: float, + fire_tlim: float, + fire_type: float, + fire_time: Union[list, np.ndarray], + fire_nft_limit: float, + fire_load_density: float, + fire_combustion_efficiency: float, + fire_hrr_density: float, + fire_spread_speed: float, + fire_t_alpha: float, + fire_gamma_fi_q: float, + beam_position_vertical: float, + beam_position_horizontal: Union[np.ndarray, list, float] = -1.0, + *_, + **__, ) -> dict: """Calculate temperature array of pre-defined fire type `fire_type`. @@ -202,9 +198,7 @@ def evaluate_fire_temperature( room_floor_area = room_breadth * room_depth # Room internal surface area, total, including window openings - room_total_area = (2 * room_floor_area) + ( - (room_breadth + room_depth) * 2 * room_height - ) + room_total_area = 2 * room_floor_area + (room_breadth + room_depth) * 2 * room_height if fire_type == 0: kwargs_fire_0_paramec = dict( @@ -224,9 +218,7 @@ def evaluate_fire_temperature( elif fire_type == 1: if beam_position_horizontal < 0: - beam_position_horizontal = np.linspace(0.5 * room_depth, room_depth, 7)[ - 1:-1 - ] + beam_position_horizontal = np.linspace(0.5 * room_depth, room_depth, 7)[1:-1] kwargs_fire_1_travel = dict( t=fire_time, @@ -242,9 +234,7 @@ def evaluate_fire_temperature( opening_height_m=window_height, opening_fraction=window_open_fraction, ) - fire_temperature, beam_position_horizontal = _fire_travelling( - **kwargs_fire_1_travel - ) + fire_temperature, beam_position_horizontal = _fire_travelling(**kwargs_fire_1_travel) if beam_position_horizontal <= 0: raise ValueError("Beam position less or equal to 0.") @@ -266,34 +256,30 @@ def evaluate_fire_temperature( else: fire_temperature = None - results = dict( + return dict( fire_temperature=fire_temperature, beam_position_horizontal=beam_position_horizontal, ) - return results - - -def solve_time_equivalence( - fire_time: Union[list, np.ndarray], - fire_temperature: Union[list, np.ndarray], - beam_cross_section_area: float, - beam_rho: float, - protection_k: float, - protection_rho: float, - protection_c: float, - protection_protected_perimeter: float, - fire_time_iso834: Union[list, np.ndarray], - fire_temperature_iso834: Union[list, np.ndarray], - solver_temperature_goal: float, - solver_max_iter: int, - solver_thickness_ubound: float, - solver_thickness_lbound: float, - solver_tol: float, - phi_teq: float, - **_, + +def solve_time_equivalence_iso834( + beam_cross_section_area: float, + beam_rho: float, + protection_k: float, + protection_rho: float, + protection_c: float, + protection_protected_perimeter: float, + fire_time_iso834: Union[list, np.ndarray], + fire_temperature_iso834: Union[list, np.ndarray], + solver_temperature_goal: float, + solver_protection_thickness: float, + phi_teq: float, + *_, + **__, ) -> dict: - """Calculates equivalent time exposure for a protected steel element member in more realistic fire environment + """ + **WIP** + Calculates equivalent time exposure for a protected steel element member in more realistic fire environment opposing to the standard fire curve ISO 834. PARAMETERS: @@ -301,13 +287,13 @@ def solve_time_equivalence( :param fire_temperature: [K], temperature array :param beam_cross_section_area: [m2], the steel beam element cross section area :param beam_rho: [kg/m3], steel beam element density - :param solver_temperature_goal: [K], steel beam element expected failure temperature :param protection_k: steel beam element protection material thermal conductivity :param protection_rho: steel beam element protection material density :param protection_c: steel beam element protection material specific heat :param protection_protected_perimeter: [m], steel beam element protection material perimeter :param fire_time_iso834: [s], the time (array) component of ISO 834 fire curve :param fire_temperature_iso834: [K], the temperature (array) component of ISO 834 fire curve + :param solver_temperature_goal: [K], steel beam element expected failure temperature :param solver_max_iter: Maximum allowable iteration counts for seeking solution for time equivalence :param solver_thickness_ubound: [m], protection layer thickness upper bound initial condition for solving time equivalence :param solver_thickness_lbound: [m], protection layer thickness lower bound initial condition for solving time equivalence @@ -323,214 +309,127 @@ def solve_time_equivalence( # MATCH PEAK STEEL TEMPERATURE BY ADJUSTING PROTECTION LAYER THICKNESS - solver_iter_count = 0 # count how many iterations for the seeking process - solver_convergence_status = ( - False - ) # flag used to indicate when the seeking is successful + # Solve equivalent time exposure in ISO 834 + solver_d_p = solver_protection_thickness + + if -np.inf < solver_d_p < np.inf: + steel_temperature = _steel_temperature( + fire_time=fire_time_iso834, + fire_temperature=fire_temperature_iso834, + beam_rho=beam_rho, + beam_cross_section_area=beam_cross_section_area, + protection_k=protection_k, + protection_rho=protection_rho, + protection_c=protection_c, + protection_thickness=solver_d_p, + protection_protected_perimeter=protection_protected_perimeter, + ) + func_teq = interp1d(steel_temperature, fire_time_iso834, kind="linear", bounds_error=False, fill_value=-1) + solver_time_equivalence_solved = func_teq(solver_temperature_goal) + solver_time_equivalence_solved = solver_time_equivalence_solved * phi_teq - # Default values - solver_time_equivalence_solved = -1 - solver_steel_temperature_solved = -1 - solver_protection_thickness = -1 + elif solver_d_p == np.inf: + solver_time_equivalence_solved = np.inf - if solver_temperature_goal > 0: # check seeking temperature, opt out if less than 0 + elif solver_d_p == -np.inf: + solver_time_equivalence_solved = -np.inf - # Solve heat transfer using EC3 correlations - # SI UNITS FOR INPUTS! - kwarg_ht_ec = dict( - time=fire_time, - temperature_ambient=fire_temperature, - rho_steel=beam_rho, - area_steel_section=beam_cross_section_area, - k_protection=protection_k, - rho_protection=protection_rho, - c_protection=protection_c, - perimeter_protected=protection_protected_perimeter, - terminate_max_temperature=solver_temperature_goal + 2 * solver_tol, - ) + elif solver_d_p is np.nan: + solver_time_equivalence_solved = np.nan - # SOLVER START FROM HERE - # ====================== - - # Minimise f(x) - θ - # where: - # f(x) is the steel maximum temperature. - # x is the thickness, - # θ is the steel temperature goal - - time_at_max_temperature = fire_time[np.argmax(fire_temperature)] - - def f_(x, terminate_check_wait_time): - kwarg_ht_ec["thickness_protection"] = x - T_ = _steel_temperature_max( - **kwarg_ht_ec, terminate_check_wait_time=terminate_check_wait_time - ) - return T_ - - # Check whether there is a solution within predefined protection thickness boundaries - x1, x2 = solver_thickness_lbound, solver_thickness_ubound - y1, y2 = f_(x1, time_at_max_temperature), f_(x2, time_at_max_temperature) - t1, t2 = ( - solver_temperature_goal - solver_tol, - solver_temperature_goal + solver_tol, - ) + else: + raise ValueError(f'This error should not occur, solver_d_p = {solver_d_p}') + + return dict(solver_time_equivalence_solved=solver_time_equivalence_solved) + + +def solve_protection_thickness( + fire_time: Union[list, np.ndarray], + fire_temperature: Union[list, np.ndarray], + beam_cross_section_area: float, + beam_rho: float, + protection_k: float, + protection_rho: float, + protection_c: float, + protection_protected_perimeter: float, + solver_temperature_goal: float, + solver_max_iter: int, + solver_thickness_ubound: float, + solver_thickness_lbound: float, + solver_tol: float, + *_, + **__, +) -> dict: + """ + Calculates equivalent time exposure for a protected steel element member in more realistic fire environment + opposing to the standard fire curve ISO 834. - if y2 <= solver_temperature_goal <= y1: - - while True: - - solver_iter_count += 1 - - # Work out linear equation: f(x) = y = a x + b - a = (y1 - y2) / (x1 - x2) - b = y1 - a * x1 - - # work out new y based upon interpolated y - x3 = solver_protection_thickness = (solver_temperature_goal - b) / a - y3 = solver_steel_temperature_solved = f_(x3, time_at_max_temperature) - - if x1 < 0 or x2 < 0 or x3 < 0: - print("check") - - if y3 < t1: # steel temperature is too low, decrease thickness - x2 = x3 - y2 = y3 - elif y3 > t2: # steel temperature is too high, increase thickness - x1 = x3 - y1 = y3 - else: - solver_convergence_status = True - - if solver_convergence_status: - - # CALCULATE BEAM FIRE RESISTANCE PERIOD IN ISO 834 - # ================================================ - # Make steel time-temperature curve when exposed to the given ambient temperature, i.e. ISO 834. - - kwarg_ht_ec["time"] = fire_time_iso834 - kwarg_ht_ec["temperature_ambient"] = fire_temperature_iso834 - kwarg_ht_ec["thickness_protection"] = x3 - steel_temperature = _steel_temperature(**kwarg_ht_ec) - - steel_time = np.concatenate( - ( - np.array([0]), - fire_time_iso834, - np.array([fire_time_iso834[-1]]), - ) - ) - steel_temperature = np.concatenate( - (np.array([-1]), steel_temperature, np.array([1e12])) - ) - func_teq = interp1d( - steel_temperature, - steel_time, - kind="linear", - bounds_error=False, - fill_value=-1, - ) - solver_time_equivalence_solved = func_teq(solver_temperature_goal) - - break - - elif ( - solver_iter_count >= solver_max_iter - ): # Terminate if maximum solving iteration is reached - - if solver_temperature_goal > y3: - solver_time_equivalence_solved = 0. - elif solver_temperature_goal < y3: - solver_time_equivalence_solved = np.inf - else: - solver_time_equivalence_solved = ( - np.nan - ) # theoretically impossible, for syntax error only - break - - elif solver_temperature_goal - 2 <= y1 <= solver_temperature_goal + 2: - solver_protection_thickness = x1 - solver_steel_temperature_solved = y1 - solver_convergence_status = True - - # CALCULATE BEAM FIRE RESISTANCE PERIOD IN ISO 834 - # ================================================ - # Make steel time-temperature curve when exposed to the given ambient temperature, i.e. ISO 834. - - kwarg_ht_ec["time"] = fire_time_iso834 - kwarg_ht_ec["temperature_ambient"] = fire_temperature_iso834 - kwarg_ht_ec["thickness_protection"] = x1 - steel_temperature = _steel_temperature(**kwarg_ht_ec) - - steel_time = np.concatenate( - (np.array([0]), fire_time_iso834, np.array([fire_time_iso834[-1]])) - ) - steel_temperature = np.concatenate( - (np.array([-1]), steel_temperature, np.array([1e12])) - ) - func_teq = interp1d( - steel_temperature, - steel_time, - kind="linear", - bounds_error=False, - fill_value=-1, - ) - solver_time_equivalence_solved = func_teq(solver_temperature_goal) - - elif solver_temperature_goal - 2 <= y2 <= solver_temperature_goal + 2: - solver_protection_thickness = x2 - solver_steel_temperature_solved = y2 - solver_convergence_status = True - - # CALCULATE BEAM FIRE RESISTANCE PERIOD IN ISO 834 - # ================================================ - # Make steel time-temperature curve when exposed to the given ambient temperature, i.e. ISO 834. - - kwarg_ht_ec["time"] = fire_time_iso834 - kwarg_ht_ec["temperature_ambient"] = fire_temperature_iso834 - kwarg_ht_ec["thickness_protection"] = x2 - steel_temperature = _steel_temperature(**kwarg_ht_ec) - - steel_time = np.concatenate( - (np.array([0]), fire_time_iso834, np.array([fire_time_iso834[-1]])) - ) - steel_temperature = np.concatenate( - (np.array([-1]), steel_temperature, np.array([1e12])) - ) - func_teq = interp1d( - steel_temperature, - steel_time, - kind="linear", - bounds_error=False, - fill_value=-1, - ) - solver_time_equivalence_solved = func_teq(solver_temperature_goal) - - # No solution, thickness upper bound is not thick enough - elif solver_temperature_goal > y1: - solver_protection_thickness = x1 - solver_steel_temperature_solved = y1 - solver_time_equivalence_solved = -np.inf - - # No solution, thickness lower bound is not thin enough - elif solver_temperature_goal < y2: - solver_protection_thickness = x2 - solver_steel_temperature_solved = y2 - solver_time_equivalence_solved = np.inf - - solver_time_equivalence_solved *= phi_teq - - results = dict( - solver_convergence_status=solver_convergence_status, - solver_time_equivalence_solved=solver_time_equivalence_solved, - solver_steel_temperature_solved=solver_steel_temperature_solved, - solver_protection_thickness=solver_protection_thickness, + PARAMETERS: + :param fire_time: [s], time array + :param fire_temperature: [K], temperature array + :param beam_cross_section_area: [m2], the steel beam element cross section area + :param beam_rho: [kg/m3], steel beam element density + :param protection_k: steel beam element protection material thermal conductivity + :param protection_rho: steel beam element protection material density + :param protection_c: steel beam element protection material specific heat + :param protection_protected_perimeter: [m], steel beam element protection material perimeter + :param fire_time_iso834: [s], the time (array) component of ISO 834 fire curve + :param fire_temperature_iso834: [K], the temperature (array) component of ISO 834 fire curve + :param solver_temperature_goal: [K], steel beam element expected failure temperature + :param solver_max_iter: Maximum allowable iteration counts for seeking solution for time equivalence + :param solver_thickness_ubound: [m], protection layer thickness upper bound initial condition for solving time equivalence + :param solver_thickness_lbound: [m], protection layer thickness lower bound initial condition for solving time equivalence + :param solver_tol: [K], tolerance for solving time equivalence + :param phi_teq: [-], model uncertainty factor + :return results: dict + EXAMPLE: + """ + + # ============================================ + # GOAL SEEK TO MATCH STEEL FAILURE TEMPERATURE + # ============================================ + + # MATCH PEAK STEEL TEMPERATURE BY ADJUSTING PROTECTION LAYER THICKNESS + + # Solve protection properties for `solver_temperature_goal` + solver_d_p, solver_T_max_a, solver_t, solver_iter_count = _protection_thickness( + fire_time=fire_time, + fire_temperature=fire_temperature, + beam_rho=beam_rho, + beam_cross_section_area=beam_cross_section_area, + protection_k=protection_k, + protection_rho=protection_rho, + protection_c=protection_c, + protection_protected_perimeter=protection_protected_perimeter, + solver_temperature_goal=solver_temperature_goal, + solver_temperature_goal_tol=solver_tol, + solver_max_iter=solver_max_iter, + d_p_1=solver_thickness_lbound, + d_p_2=solver_thickness_ubound, + ) + + return dict( + solver_convergence_status=-np.inf < solver_d_p < np.inf, + solver_steel_temperature_solved=solver_T_max_a, + solver_time_solved=solver_t, + solver_protection_thickness=solver_d_p, solver_iter_count=solver_iter_count, ) - return results +def mcs_out_post_per_case(df: pd.DataFrame, fp: str) -> pd.DataFrame: + # save outputs if work direction is provided per iteration + if fp: + def _save_(fp_: str): + try: + if not os.path.exists(os.path.dirname(fp_)): + os.makedirs(os.path.dirname(fp_)) + except Exception as e: + print(e) + df.to_csv(os.path.join(fp_), index=False) + + threading.Thread(target=_save_, kwargs=dict(fp_=fp)).start() -def mcs_out_post(df: pd.DataFrame) -> pd.DataFrame: df_res = copy.copy(df) df_res = df_res.replace(to_replace=[np.inf, -np.inf], value=np.nan) df_res = df_res.dropna(axis=0, how="any") @@ -572,7 +471,7 @@ def mcs_out_post(df: pd.DataFrame) -> pd.DataFrame: y = np.linspace(0, 1, len(x), dtype=float) aplot = AsciiPlot(size=(55, 15)) aplot.plot(x=x, y=y, xlim=(20, min([180, np.amax(x)]))) - print(aplot.show()) + aplot.show() except Exception as e: print(f'Failed to plot time equivalence, {e}') @@ -589,50 +488,50 @@ def teq_main_wrapper(args): def teq_main( - case_name: str, - n_simulations: int, - probability_weight: float, - index: int, - beam_cross_section_area: float, - beam_position_vertical: float, - beam_position_horizontal: float, - beam_rho: float, - fire_time_duration: float, - fire_time_step: float, - fire_combustion_efficiency: float, - fire_gamma_fi_q: float, - fire_hrr_density: float, - fire_load_density: float, - fire_mode: int, - fire_nft_limit: float, - fire_spread_speed: float, - fire_t_alpha: float, - fire_tlim: float, - protection_c: float, - protection_k: float, - protection_protected_perimeter: float, - protection_rho: float, - room_breadth: float, - room_depth: float, - room_height: float, - room_wall_thermal_inertia: float, - solver_temperature_goal: float, - solver_max_iter: int, - solver_thickness_lbound: float, - solver_thickness_ubound: float, - solver_tol: float, - window_height: float, - window_open_fraction: float, - window_width: float, - window_open_fraction_permanent: float, - phi_teq: float = 1.0, - timber_charring_rate=None, - timber_hc: float = None, - timber_density: float = None, - timber_exposed_area: float = None, - timber_solver_tol: float = None, - timber_solver_ilim: float = None, - **_, + case_name: str, + n_simulations: int, + index: int, + beam_cross_section_area: float, + beam_position_vertical: float, + beam_position_horizontal: float, + beam_rho: float, + fire_time_duration: float, + fire_time_step: float, + fire_combustion_efficiency: float, + fire_gamma_fi_q: float, + fire_hrr_density: float, + fire_load_density: float, + fire_mode: int, + fire_nft_limit: float, + fire_spread_speed: float, + fire_t_alpha: float, + fire_tlim: float, + protection_c: float, + protection_k: float, + protection_protected_perimeter: float, + protection_rho: float, + room_breadth: float, + room_depth: float, + room_height: float, + room_wall_thermal_inertia: float, + solver_temperature_goal: float, + solver_max_iter: int, + solver_thickness_lbound: float, + solver_thickness_ubound: float, + solver_tol: float, + window_height: float, + window_open_fraction: float, + window_width: float, + window_open_fraction_permanent: float, + phi_teq: float = 1.0, + timber_charring_rate=None, + timber_hc: float = None, + timber_density: float = None, + timber_exposed_area: float = None, + timber_solver_tol: float = None, + timber_solver_ilim: float = None, + *_, + **__, ) -> dict: # Make the longest dimension between (room_depth, room_breadth) as room_depth if room_depth < room_breadth: @@ -640,196 +539,108 @@ def teq_main( room_breadth = room_depth - room_breadth room_depth -= room_breadth - window_open_fraction = ( - window_open_fraction * (1 - window_open_fraction_permanent) - + window_open_fraction_permanent - ) + window_open_fraction = (window_open_fraction * (1 - window_open_fraction_permanent) + window_open_fraction_permanent) # Fix ventilation opening size so it doesn't exceed wall area - if window_height > room_height: window_height = room_height - room_perimeter = 2 * (room_breadth + room_depth) - - if window_width > room_perimeter: - window_width = room_perimeter - # Calculate fire time, this is used for all fire curves in the calculation - fire_time = np.arange(0, fire_time_duration + fire_time_step, fire_time_step) # Calculate ISO 834 fire temperature - fire_time_iso834 = fire_time - fire_temperature_iso834 = ( - 345.0 * np.log10((fire_time / 60.0) * 8.0 + 1.0) + 20.0 - ) + 273.15 # in [K] - - # Inject results, i.e. these input values will also be in the output file - res = dict( - case_name=case_name, - n_simulations=n_simulations, - probability_weight=probability_weight, - index=index, - # beam_cross_section_area=beam_cross_section_area, - # beam_position_vertical=beam_position_vertical, - beam_position_horizontal=beam_position_horizontal, - # beam_rho=beam_rho, - # fire_time=fire_time, - fire_combustion_efficiency=fire_combustion_efficiency, - # fire_gamma_fi_q=fire_gamma_fi_q, - fire_hrr_density=fire_hrr_density, - fire_load_density=fire_load_density, - # fire_mode=fire_mode, - fire_nft_limit=fire_nft_limit, - fire_spread_speed=fire_spread_speed, - # fire_t_alpha=fire_t_alpha, - # fire_tlim=fire_tlim, - # fire_temperature_iso834=fire_temperature_iso834, - # fire_time_iso834=fire_time_iso834, - # protection_c=protection_c, - # protection_k=protection_k, - # protection_protected_perimeter=protection_protected_perimeter, - # protection_rho=protection_rho, - room_breadth=room_breadth, - room_depth=room_depth, - # room_height=room_height, - # room_wall_thermal_inertia=room_wall_thermal_inertia, - # solver_temperature_goal=solver_temperature_goal, - # solver_max_iter=solver_max_iter, - # solver_thickness_lbound=solver_thickness_lbound, - # solver_thickness_ubound=solver_thickness_ubound, - # solver_tol=solver_tol, - window_height=window_height, - window_open_fraction=window_open_fraction, - window_width=window_width, - phi_teq=phi_teq, - # timber_charring_rate=timber_charring_rate, - # timber_hc=timber_hc, - # timber_density=timber_density, - # timber_exposed_area=timber_exposed_area, - # timber_solver_tol=timber_solver_tol, - # timber_solver_ilim=timber_solver_ilim, - ) - - # initial timber exposure time - if timber_exposed_area > 0: - timber_exposed_duration = 1200 - else: - timber_exposed_duration = 0 + fire_temperature_iso834 = (345.0 * np.log10((fire_time / 60.0) * 8.0 + 1.0) + 20.0) + 273.15 # in [K] - timber_solver_iter = ( - -1 - ) # initialise solver iteration count for timber fuel contribution + inputs = copy.deepcopy(locals()) + inputs.pop('_'), inputs.pop('__') - if isinstance(timber_charring_rate, (int, float)): - timber_charring_rate_ = copy.copy(timber_charring_rate) - timber_charring_rate = lambda x: timber_charring_rate_ + # initialise solver iteration count for timber fuel contribution + timber_solver_iter_count = -1 + timber_exposed_duration = 0 # initial condition, timber exposed duration + _fire_load_density_ = inputs.pop('fire_load_density') # preserve original fire load density while True: - timber_solver_iter += 1 - timber_charring_rate_ = timber_charring_rate(timber_exposed_duration) - timber_charring_rate_ *= 1 / 1000 # [mm/min] -> [m/min] - timber_charring_rate_ *= 1 / 60 # [m/min] -> [m/s] - timber_charred_depth = timber_charring_rate_ * timber_exposed_duration + timber_solver_iter_count += 1 + if isinstance(timber_charring_rate, (float, int)): + timber_charring_rate_i = timber_charring_rate + elif isinstance(timber_charring_rate, Callable): + timber_charring_rate_i = timber_charring_rate(timber_exposed_duration) + else: + raise TypeError('`timber_charring_rate_i` is not numerical nor Callable type') + timber_charring_rate_i *= 1 / 1000 # [mm/min] -> [m/min] + timber_charring_rate_i *= 1 / 60 # [m/min] -> [m/s] + timber_charred_depth = timber_charring_rate_i * timber_exposed_duration timber_charred_volume = timber_charred_depth * timber_exposed_area timber_charred_mass = timber_density * timber_charred_volume timber_fire_load = timber_charred_mass * timber_hc timber_fire_load_density = timber_fire_load / (room_breadth * room_depth) - # To check what design fire to use + inputs['fire_load_density'] = _fire_load_density_ + timber_fire_load_density - res_decide_fire = decide_fire( - window_height, - window_width, - window_open_fraction, - room_breadth, - room_depth, - room_height, - fire_mode, - fire_load_density + timber_fire_load_density, - fire_combustion_efficiency, - fire_hrr_density, - fire_spread_speed, - ) + # To check what design fire to use + inputs.update(decide_fire(**inputs)) # To calculate design fire temperature + inputs.update(evaluate_fire_temperature(**inputs)) - res_evaluate_fire_temperature = evaluate_fire_temperature( - window_height, - window_width, - window_open_fraction, - room_breadth, - room_depth, - room_height, - room_wall_thermal_inertia, - fire_tlim, - res_decide_fire["fire_type"], - fire_time, - fire_nft_limit, - fire_load_density + timber_fire_load_density, - fire_combustion_efficiency, - fire_hrr_density, - fire_spread_speed, - fire_t_alpha, - fire_gamma_fi_q, - beam_position_vertical, - beam_position_horizontal, - ) - - # To calculate time equivalence - - res_solve_time_equivalence = solve_time_equivalence( - fire_time, - res_evaluate_fire_temperature["fire_temperature"], - beam_cross_section_area, - beam_rho, - protection_k, - protection_rho, - protection_c, - protection_protected_perimeter, - fire_time_iso834, - fire_temperature_iso834, - solver_temperature_goal, - solver_max_iter, - solver_thickness_ubound, - solver_thickness_lbound, - solver_tol, - phi_teq, - ) + # To solve protection thickness at critical temperature + inputs.update(solve_protection_thickness(**inputs)) # additional fuel contribution from timber - if timber_exposed_area <= 0 or timber_exposed_area is None: # no timber exposed + # Exit timber fuel contribution solver if: + # 1. no timber exposed + # 2. timber exposed area undefined break - elif not res_solve_time_equivalence["solver_convergence_status"]: # no time equivalence solution + elif timber_solver_iter_count >= timber_solver_ilim: + inputs['solver_convergence_status'] = np.nan + inputs['solver_steel_temperature_solved'] = np.nan + inputs['solver_time_solved'] = np.nan + inputs['solver_protection_thickness'] = np.nan + inputs['solver_iter_count'] = np.nan + timber_exposed_duration = np.nan break - elif (timber_solver_iter >= timber_solver_ilim): # over the solver iteration limit + elif not -np.inf < inputs["solver_protection_thickness"] < np.inf: + # no protection thickness solution + timber_exposed_duration = inputs['solver_protection_thickness'] break - elif (abs(timber_exposed_duration- res_solve_time_equivalence["solver_time_equivalence_solved"])<= timber_solver_tol): # convergence sought successfully + elif abs(timber_exposed_duration - inputs["solver_time_solved"]) <= timber_solver_tol: + # convergence sought successfully break else: - timber_exposed_duration = res_solve_time_equivalence[ - "solver_time_equivalence_solved" - ] - - res_timber_solver = dict( - timber_charring_rate=timber_charring_rate_, - timber_exposed_duration=timber_exposed_duration, - timber_solver_iter=timber_solver_iter, - timber_fire_load=timber_fire_load, - timber_charred_depth=timber_charred_depth, - timber_charred_mass=timber_charred_mass, - timber_charred_volume=timber_charred_volume, + timber_exposed_duration = inputs["solver_time_solved"] + + inputs.update(solve_time_equivalence_iso834(**inputs)) + + inputs.update( + dict( + timber_charring_rate=timber_charring_rate_i, + timber_exposed_duration=timber_exposed_duration, + timber_solver_iter_count=timber_solver_iter_count, + timber_fire_load=timber_fire_load, + timber_charred_depth=timber_charred_depth, + timber_charred_mass=timber_charred_mass, + timber_charred_volume=timber_charred_volume, + ) ) - res.update(res_timber_solver) - res.update(res_decide_fire) - res.update(res_evaluate_fire_temperature) - res.update(res_solve_time_equivalence) + # Prepare results to be returned, only the items in the list below will be returned + # add keys accordingly if more parameters are desired to be returned + outputs = { + i: inputs[i] for i in + ['phi_teq', 'fire_spread_speed', 'fire_nft_limit', 'fire_mode', 'fire_load_density', 'fire_hrr_density', 'fire_combustion_efficiency', 'beam_position_horizontal', + # 'beam_position_vertical', 'index', 'probability_weight', 'case_name', 'fire_type', 'solver_convergence_status', 'solver_time_equivalence_solved', + 'beam_position_vertical', 'index', 'case_name', 'fire_type', 'solver_convergence_status', 'solver_time_equivalence_solved', + 'solver_steel_temperature_solved', 'solver_protection_thickness', 'solver_iter_count', 'window_open_fraction', 'timber_solver_iter_count', 'timber_charred_depth'] + } + + return outputs + - return res +def mcs_out_post_all_cases(df: pd.DataFrame, fp: str): + if fp: + df[['case_name', 'index', 'solver_time_equivalence_solved']].to_csv(fp, index=False) class MCS0(MCS): @@ -842,8 +653,26 @@ def mcs_deterministic_calc(self, *args, **kwargs) -> dict: def mcs_deterministic_calc_mp(self, *args, **kwargs) -> dict: return teq_main_wrapper(*args, **kwargs) - def mcs_post(self, *args, **kwargs): - return mcs_out_post(*args, **kwargs) + def mcs_post_per_case(self, df: pd.DataFrame): + + case_name = df['case_name'].to_numpy() + assert (case_name == case_name[0]).all() + case_name = case_name[0] + + try: + fp = os.path.join(self.cwd, self.DEFAULT_TEMP_FOLDER_NAME, f'{case_name}.csv') + except TypeError: + fp = None + + return mcs_out_post_per_case(df=df, fp=fp) + + def mcs_post_all_cases(self, df: pd.DataFrame): + try: + fp = os.path.join(self.cwd, self.DEFAULT_MCS_OUTPUT_FILE_NAME) + except TypeError: + fp = None + + return mcs_out_post_all_cases(df=df, fp=fp) def _test_teq_phi(): @@ -857,14 +686,14 @@ def _test_teq_phi(): input_param = dict( index=0, case_name="Standard 1", - probability_weight=1, - fire_time_step=30, - fire_time_duration=5 * 60 * 60, + probability_weight=1., + fire_time_step=1., + fire_time_duration=5. * 60 * 60, n_simulations=1, beam_cross_section_area=0.017, beam_position_vertical=2.5, beam_position_horizontal=18, - beam_rho=7850, + beam_rho=7850., fire_combustion_efficiency=0.8, fire_gamma_fi_q=1, fire_hrr_density=0.25, @@ -876,19 +705,19 @@ def _test_teq_phi(): fire_tlim=0.333, fire_temperature_iso834=fire_temperature_iso834_, fire_time_iso834=fire_time_, - protection_c=1700, + protection_c=1700., protection_k=0.2, protection_protected_perimeter=2.14, - protection_rho=7850, + protection_rho=800., room_breadth=16, room_depth=31.25, room_height=3, room_wall_thermal_inertia=720, solver_temperature_goal=620 + 273.15, - solver_max_iter=20, + solver_max_iter=200, solver_thickness_lbound=0.0001, solver_thickness_ubound=0.0500, - solver_tol=1.0, + solver_tol=0.01, window_height=2, window_open_fraction=0.8, window_width=72, @@ -908,136 +737,59 @@ def _test_teq_phi(): input_param["phi_teq"] = 0.1 teq_01 = teq_main(**input_param)["solver_time_equivalence_solved"] - assert abs(teq_10 / teq_01 - 10) < 0.001 + print( + f'Time equivalence at phi_teq=0.1: {teq_01:<8.3f}\n' + f'Time equivalence at phi_teq=1.0: {teq_10:<8.3f}\n' + f'Ratio between the above: {teq_10 / teq_01:<8.3f}\n' + ) + + assert abs(teq_10 / teq_01 - 10) < 0.01 def _test_standard_case(): import copy - from sfeprapy.func.mcs_obj import MCS from sfeprapy.mcs0 import EXAMPLE_INPUT_DICT, EXAMPLE_CONFIG_DICT - # from sfeprapy.mcs0.mcs0_calc import teq_main, teq_main_wrapper, mcs_out_post - from sfeprapy.func.mcs_gen import main as gen from scipy.interpolate import interp1d - import numpy as np # increase the number of simulations so it gives sensible results mcs_input = copy.deepcopy(EXAMPLE_INPUT_DICT) mcs_config = copy.deepcopy(EXAMPLE_CONFIG_DICT) - for k in list(mcs_input.keys()): - mcs_input[k]["phi_teq"] = 1 - mcs_input[k]["n_simulations"] = 333 - mcs_input[k]["probability_weight"] = 1 / 3.0 - mcs_input[k]["fire_time_duration"] = 10000 - mcs_input[k]["timber_exposed_area"] = 0 - mcs_input[k].pop("beam_position_horizontal") - mcs_input[k]["beam_position_horizontal:dist"] = "uniform_" - mcs_input[k]["beam_position_horizontal:ubound"] = ( - mcs_input[k]["room_depth"] * 0.9 - ) - mcs_input[k]["beam_position_horizontal:lbound"] = ( - mcs_input[k]["room_depth"] * 0.6 - ) - - # increase the number of threads so it runs faster - mcs_config["n_threads"] = 1 # coverage does not support - mcs = MCS() - mcs.define_problem(data=mcs_input, config=mcs_config) - mcs.define_stochastic_parameter_generator(gen) - mcs.define_calculation_routine(teq_main, teq_main_wrapper, mcs_out_post) - mcs.run_mcs() - mcs_out = mcs.mcs_out - teq = mcs_out["solver_time_equivalence_solved"] / 60.0 - hist, edges = np.histogram(teq, bins=np.arange(0, 181, 0.5)) - x, y = (edges[:-1] + edges[1:]) / 2, np.cumsum(hist / np.sum(hist)) - teq_at_80_percentile = interp1d(y, x)(0.8) - print(teq_at_80_percentile) - target, target_tol = 60, 2 - assert target - target_tol < teq_at_80_percentile < target + target_tol - + mcs_config["n_threads"] = 1 -def _test_standard_oversized_case(): - """the same as standard case but window_height and window_width are set to over the floor limit.""" - import copy - from sfeprapy.func.mcs_obj import MCS - from sfeprapy.mcs0 import EXAMPLE_INPUT_DICT, EXAMPLE_CONFIG_DICT - # from sfeprapy.mcs0.mcs0_calc import teq_main, teq_main_wrapper, mcs_out_post - from sfeprapy.func.mcs_gen import main as gen - from scipy.interpolate import interp1d - import numpy as np + mcs = MCS0() - # increase the number of simulations so it gives sensible results - mcs_input = copy.deepcopy(EXAMPLE_INPUT_DICT) - mcs_config = copy.deepcopy(EXAMPLE_CONFIG_DICT) - for k in list(mcs_input.keys()): - mcs_input[k]["phi_teq"] = 1 - mcs_input[k]["n_simulations"] = 10 - mcs_input[k]["probability_weight"] = 1 / 3.0 - mcs_input[k]["fire_time_duration"] = 10000 - mcs_input[k]["timber_exposed_area"] = 0 - mcs_input[k].pop("beam_position_horizontal") - mcs_input[k]["beam_position_horizontal:dist"] = "uniform_" - mcs_input[k]["beam_position_horizontal:ubound"] = ( - mcs_input[k]["room_depth"] * 0.9 - ) - mcs_input[k]["beam_position_horizontal:lbound"] = ( - mcs_input[k]["room_depth"] * 0.6 - ) - mcs_input[k]['window_width'] = 500 - mcs_input[k]['window_height'] = 5 - - # increase the number of threads so it runs faster - mcs_config["n_threads"] = 1 # coverage does not support - mcs = MCS() - mcs.define_problem(data=mcs_input, config=mcs_config) - mcs.define_stochastic_parameter_generator(gen) - mcs.define_calculation_routine(teq_main, teq_main_wrapper, mcs_out_post) + mcs.mcs_inputs = mcs_input + mcs.mcs_config = mcs_config mcs.run_mcs() mcs_out = mcs.mcs_out - teq = mcs_out["solver_time_equivalence_solved"] / 60.0 - hist, edges = np.histogram(teq, bins=np.arange(0, 181, 0.5)) - x, y = (edges[:-1] + edges[1:]) / 2, np.cumsum(hist / np.sum(hist)) - teq_at_80_percentile = interp1d(y, x)(0.8) - - assert len(mcs_out['window_height']) == sum(mcs_out['window_height'].values == 3.3) - assert len(mcs_out['window_width']) == sum(mcs_out['window_width'].values == 94.5) + def get_time_equivalence(data, fractile: float): + hist, edges = np.histogram(data, bins=np.arange(0, 181, 0.5)) + x, y = (edges[:-1] + edges[1:]) / 2, np.cumsum(hist / np.sum(hist)) + return interp1d(y, x)(fractile) -def _test_standard_case_new(): - import copy - from sfeprapy.mcs0 import EXAMPLE_INPUT_DICT, EXAMPLE_CONFIG_DICT - # from sfeprapy.mcs0.mcs0_calc import teq_main, teq_main_wrapper, mcs_out_post - from scipy.interpolate import interp1d - import numpy as np + mcs_out_standard_case_1 = mcs_out.loc[mcs_out['case_name'] == 'Standard Case 1'] + teq = mcs_out_standard_case_1["solver_time_equivalence_solved"] / 60.0 + teq_at_80_percentile = get_time_equivalence(teq, 0.8) + print(f'Time equivalence at CDF 0.8 is {teq_at_80_percentile:<6.3f} min') + target, target_tol = 60, 2 + assert target - target_tol < teq_at_80_percentile < target + target_tol - # increase the number of simulations so it gives sensible results - mcs_input = copy.deepcopy(EXAMPLE_INPUT_DICT) - mcs_config = copy.deepcopy(EXAMPLE_CONFIG_DICT) - for k in list(mcs_input.keys()): - mcs_input[k]["phi_teq"] = 1 - mcs_input[k]["n_simulations"] = 333 - mcs_input[k]["probability_weight"] = 1 / 3.0 - mcs_input[k]["fire_time_duration"] = 10000 - mcs_input[k]["timber_exposed_area"] = 0 - mcs_input[k].pop("beam_position_horizontal") - mcs_input[k]["beam_position_horizontal:dist"] = "uniform_" - mcs_input[k]["beam_position_horizontal:ubound"] = ( - mcs_input[k]["room_depth"] * 0.9 - ) - mcs_input[k]["beam_position_horizontal:lbound"] = ( - mcs_input[k]["room_depth"] * 0.6 - ) + mcs_out_standard_case_2 = mcs_out.loc[mcs_out['case_name'] == 'Standard Case 2 (with teq_phi)'] + teq = mcs_out_standard_case_2["solver_time_equivalence_solved"] / 60.0 + teq_at_80_percentile = get_time_equivalence(teq, 0.8) + print(f'Time equivalence at CDF 0.8 is {teq_at_80_percentile:<6.3f} min') + target, target_tol = 64, 2 # 64 minutes based on a test run on 2nd Oct 2020 + assert target - target_tol < teq_at_80_percentile < target + target_tol - # increase the number of threads so it runs faster - mcs_config["n_threads"] = 1 # coverage does not support - mcs0 = MCS0() - mcs0.mcs_inputs = mcs_input - mcs0.mcs_config = EXAMPLE_CONFIG_DICT - mcs0.run_mcs() - mcs_out = mcs0.mcs_out - teq = mcs_out["solver_time_equivalence_solved"] / 60.0 - hist, edges = np.histogram(teq, bins=np.arange(0, 181, 0.5)) - x, y = (edges[:-1] + edges[1:]) / 2, np.cumsum(hist / np.sum(hist)) - teq_at_80_percentile = interp1d(y, x)(0.8) - print(teq_at_80_percentile) - target, target_tol = 60, 2 + mcs_out_standard_case_3 = mcs_out.loc[mcs_out['case_name'] == 'Standard Case 3 (with timber)'] + teq = mcs_out_standard_case_3["solver_time_equivalence_solved"] / 60.0 + teq_at_80_percentile = get_time_equivalence(teq, 0.8) + print(f'Time equivalence at CDF 0.8 is {teq_at_80_percentile:<6.3f} min') + target, target_tol = 90, 2 # 81 minutes based on a test run on 2nd Oct 2020 assert target - target_tol < teq_at_80_percentile < target + target_tol + + +if __name__ == '__main__': + _test_teq_phi() + _test_standard_case() diff --git a/sfeprapy/mcs1/mcs1_func_main.py b/sfeprapy/mcs1/mcs1_func_main.py index 4ee65f7..edcd66d 100644 --- a/sfeprapy/mcs1/mcs1_func_main.py +++ b/sfeprapy/mcs1/mcs1_func_main.py @@ -1,48 +1,42 @@ # -*- coding: utf-8 -*- import numpy as np import pandas as pd - -from sfeprapy.dat.ec_3_1_2_kyT import ky2T_probabilistic_vectorised -from sfeprapy.func.fire_parametric_ec import fire as _fire_param -from sfeprapy.func.fire_parametric_ec_din import fire as _fire_param_ger -from sfeprapy.func.fire_travelling import fire as _fire_travelling -from sfeprapy.func.heat_transfer_protected_steel_ec import ( - protected_steel_eurocode_max_temperature as _steel_temperature_max, -) +from fsetools.lib.fse_bs_en_1991_1_2_parametric_fire import temperature as _fire_param +from fsetools.lib.fse_bs_en_1993_1_2_heat_transfer import temperature_max as _steel_temperature_max +from fsetools.lib.fse_bs_en_1993_1_2_strength_reduction_factor import k_y_theta_prob +from fsetools.lib.fse_din_en_1991_1_2_parametric_fire import temperature as _fire_param_ger +from fsetools.lib.fse_travelling_fire import temperature as _fire_travelling def a_solve_steel_protection_thickness_in_iso_fire( - fire_iso834_time, - fire_iso834_temperature, - beam_rho, - beam_cross_section_area, - protection_k, - protection_rho, - protection_c, - protection_protected_perimeter, - solver_temperature_target, - solver_fire_duration, - solver_thickness_lbound, - solver_thickness_ubound, - solver_tolerance, - solver_iteration_limit, - **_ + fire_iso834_time, + fire_iso834_temperature, + beam_rho, + beam_cross_section_area, + protection_k, + protection_rho, + protection_c, + protection_protected_perimeter, + solver_temperature_target, + solver_fire_duration, + solver_thickness_lbound, + solver_thickness_ubound, + solver_tolerance, + solver_iteration_limit, + **_ ): # Solve heat transfer using EC3 correlations # SI UNITS FOR INPUTS! kwarg_ht_ec = dict( - time=fire_iso834_time[fire_iso834_time <= solver_fire_duration], - temperature_ambient=fire_iso834_temperature[ - fire_iso834_time <= solver_fire_duration - ], - rho_steel=beam_rho, - area_steel_section=beam_cross_section_area, - k_protection=protection_k, - rho_protection=protection_rho, - c_protection=protection_c, - perimeter_protected=protection_protected_perimeter, - terminate_max_temperature=solver_temperature_target + 5 * solver_tolerance, + fire_time=fire_iso834_time[fire_iso834_time <= solver_fire_duration], + fire_temperature=fire_iso834_temperature[fire_iso834_time <= solver_fire_duration], + beam_rho=beam_rho, + beam_cross_section_area=beam_cross_section_area, + protection_k=protection_k, + protection_rho=protection_rho, + protection_c=protection_c, + protection_protected_perimeter=protection_protected_perimeter, ) solver_iteration_count = 0 # count how many iterations for the seeking process @@ -52,7 +46,7 @@ def a_solve_steel_protection_thickness_in_iso_fire( solver_thickness_solved = -1 if ( - solver_temperature_target > 0 + solver_temperature_target > 0 ): # check seeking temperature, opt out if less than 0 # def f_(x): @@ -66,13 +60,11 @@ def a_solve_steel_protection_thickness_in_iso_fire( # y1, y2 = f_(x1), f_(x2) y1 = _steel_temperature_max( **kwarg_ht_ec, - thickness_protection=x1, - terminate_check_wait_time=solver_fire_duration - 30 + protection_thickness=x1, ) y2 = _steel_temperature_max( **kwarg_ht_ec, - thickness_protection=x2, - terminate_check_wait_time=solver_fire_duration - 30 + protection_thickness=x2, ) t1, t2 = ( @@ -81,9 +73,9 @@ def a_solve_steel_protection_thickness_in_iso_fire( ) if ( - (y2 - solver_tolerance) - <= solver_temperature_target - <= (y1 + solver_tolerance) + (y2 - solver_tolerance) + <= solver_temperature_target + <= (y1 + solver_tolerance) ): while True: @@ -99,8 +91,7 @@ def a_solve_steel_protection_thickness_in_iso_fire( # y3 = f_(x3) y3 = _steel_temperature_max( **kwarg_ht_ec, - thickness_protection=x3, - terminate_check_wait_time=solver_fire_duration - 30 + protection_thickness=x3, ) if y3 < t1: # steel temperature is too low, decrease thickness @@ -113,7 +104,7 @@ def a_solve_steel_protection_thickness_in_iso_fire( flag_solver_status = True if flag_solver_status or ( - solver_iteration_count >= solver_iteration_limit + solver_iteration_count >= solver_iteration_limit ): break @@ -129,32 +120,32 @@ def a_solve_steel_protection_thickness_in_iso_fire( def b_calc_steel_temperature_in_design_fire( - fire_time, - fire_mode, - room_depth, - room_breadth, - room_height, - room_wall_thermal_inertia, - window_height, - window_width, - window_open_fraction, - beam_loc_z, - beam_position, - beam_rho, - beam_cross_section_area, - protection_k, - protection_rho, - protection_c, - protection_protected_perimeter, - protection_thickness, - fire_load_density, - fire_hrr_density, - fire_tlim, - fire_spread_speed, - fire_nft_ubound, - fire_t_alpha, - fire_gamma_fi_q, - **_ + fire_time, + fire_mode, + room_depth, + room_breadth, + room_height, + room_wall_thermal_inertia, + window_height, + window_width, + window_open_fraction, + beam_loc_z, + beam_position, + beam_rho, + beam_cross_section_area, + protection_k, + protection_rho, + protection_c, + protection_protected_perimeter, + protection_thickness, + fire_load_density, + fire_hrr_density, + fire_tlim, + fire_spread_speed, + fire_nft_ubound, + fire_t_alpha, + fire_gamma_fi_q, + **_ ): # Make the longest dimension between (room_depth, room_breadth) as room_depth room_depth, room_breadth = ( @@ -170,7 +161,7 @@ def b_calc_steel_temperature_in_design_fire( # Room internal surface area, total, including window openings room_total_area = (2 * room_floor_area) + ( - (room_breadth + room_depth) * 2 * room_height + (room_breadth + room_depth) * 2 * room_height ) # Fire load density related to the total surface area A_t @@ -243,9 +234,9 @@ def b_calc_steel_temperature_in_design_fire( elif fire_mode == 3: # enforced to ec parametric + travelling if ( - fire_spread_entire_room_time < burn_out_time - and 0.02 < opening_factor <= 0.2 - and 50 <= fire_load_density_total <= 1000 + fire_spread_entire_room_time < burn_out_time + and 0.02 < opening_factor <= 0.2 + and 50 <= fire_load_density_total <= 1000 ): # If fire spreads throughout compartment and ventilation is within EC limits = Parametric fire fire_temp = _fire_param(**kwargs_fire_0_paramec) @@ -259,8 +250,8 @@ def b_calc_steel_temperature_in_design_fire( elif fire_mode == 4: # enforced to german parametric + travelling if ( - fire_spread_entire_room_time < burn_out_time - and 0.125 <= (window_area / room_floor_area) <= 0.5 + fire_spread_entire_room_time < burn_out_time + and 0.125 <= (window_area / room_floor_area) <= 0.5 ): # If fire spreads throughout compartment and ventilation is within EC limits = Parametric fire fire_temp = _fire_param_ger(**kwargs_fire_2_paramdin) @@ -297,23 +288,19 @@ def b_calc_steel_temperature_in_design_fire( def c_strength_reduction_factor( - solver_steel_temperature_solved: np.ndarray, is_random_q: bool = True + solver_steel_temperature_solved: np.ndarray, is_random_q: bool = True ): - if is_random_q: epsilon_q = np.random.random_sample(len(solver_steel_temperature_solved)) else: epsilon_q = np.full(len(solver_steel_temperature_solved), 0.5) - steel_strength_reduction_factor = ky2T_probabilistic_vectorised( - T=solver_steel_temperature_solved, epsilon_q=epsilon_q - ) + steel_strength_reduction_factor = k_y_theta_prob(theta_a=solver_steel_temperature_solved, epsilon_q=epsilon_q) return steel_strength_reduction_factor def main(df_mc_params: pd.DataFrame): - dict_mc_params = df_mc_params.to_dict(orient="index") # ================================================================================================================== @@ -337,88 +324,25 @@ def main(df_mc_params: pd.DataFrame): ) for uid, kwargs in dict_mc_params.items(): - - solver_steel_temperature_solved, fire_type = b_calc_steel_temperature_in_design_fire( - protection_thickness=solver_thickness_solved, **kwargs - ) + solver_steel_temperature_solved, fire_type = b_calc_steel_temperature_in_design_fire(protection_thickness=solver_thickness_solved, **kwargs) dict_mc_params[uid]["solver_thickness_solved"] = solver_thickness_solved dict_mc_params[uid]["flag_solver_status"] = flag_solver_status dict_mc_params[uid]["solver_iteration_count"] = solver_iteration_count dict_mc_params[uid]["fire_type"] = fire_type - dict_mc_params[uid][ - "solver_steel_temperature_solved" - ] = solver_steel_temperature_solved + dict_mc_params[uid]["solver_steel_temperature_solved"] = solver_steel_temperature_solved df_out = pd.DataFrame.from_dict(dict_mc_params, orient="index") df_out.set_index("index", inplace=True) # assign 'index' column as DataFrame index - list_c_strength_reduction_factor = c_strength_reduction_factor( - df_out["solver_steel_temperature_solved"].values, True - ) + list_c_strength_reduction_factor = c_strength_reduction_factor(df_out["solver_steel_temperature_solved"].values, True) df_out["strength_reduction_factor"] = list_c_strength_reduction_factor - # TODO: HASH ITEMS BELOW TO REMOVE FROM OUTPUT CSV - df_out.pop("fire_time") - df_out.pop("fire_iso834_time") - df_out.pop("fire_iso834_temperature") - df_out.pop("window_height") - df_out.pop("window_width") - # df_mcs_in.pop("room_breadth") - # df_mcs_in.pop("room_depth") - # df_mcs_in.pop("room_height") - df_out.pop("room_wall_thermal_inertia") - df_out.pop("fire_mode") - df_out.pop("fire_time_step") - df_out.pop("fire_tlim") - df_out.pop("fire_hrr_density") - df_out.pop("fire_duration") - df_out.pop("fire_t_alpha") - df_out.pop("fire_gamma_fi_q") - df_out.pop("beam_rho") - df_out.pop("beam_cross_section_area") - df_out.pop("beam_loc_z") - df_out.pop("protection_k") - df_out.pop("protection_rho") - df_out.pop("protection_c") - df_out.pop("protection_protected_perimeter") - df_out.pop("solver_temperature_target") - df_out.pop("solver_fire_duration") - df_out.pop("solver_thickness_lbound") - df_out.pop("solver_thickness_ubound") - df_out.pop("solver_tolerance") - df_out.pop("solver_iteration_limit") - # df_mcs_in.pop("solver_thickness_solved") - # df_mcs_in.pop("flag_solver_status") - # df_mcs_in.pop("solver_iteration_count") - # df_mcs_in.pop("fire_type") + # TODO: HASH ITEMS BELOW TO REMOVE FROM OUTPUT CSV, + for k in ("fire_time", "fire_iso834_time", "fire_iso834_temperature", "window_height", "window_width", "room_wall_thermal_inertia", "fire_mode", "fire_time_step", "fire_tlim", + "fire_hrr_density", "fire_duration", "fire_t_alpha", "fire_gamma_fi_q", "beam_rho", "beam_cross_section_area", "beam_loc_z", "protection_k", "protection_rho", + "protection_c", "protection_protected_perimeter", "solver_temperature_target", "solver_fire_duration", "solver_thickness_lbound", "solver_thickness_ubound", + "solver_tolerance", "solver_iteration_limit"): + df_out.pop(k) return df_out - - -# if self.n_threads <= 1: -# results = [] -# for dict_mc_params in mc_param_list: -# results.append(calc_main(**dict_mc_params)) -# -# else: -# time_simulation_start = time.perf_counter() -# m = mp.Manager() -# q = m.Queue() -# p = mp.Pool(n_threads, maxtasksperchild=1000) -# jobs = p.map_async(calc_time_equivalence_worker, [(dict_mc_param, q) for dict_mc_param in mc_param_list]) -# n_steps = 24 # length of the progress bar -# while True: -# if jobs.ready(): -# time_simulation_consumed = time.perf_counter() - time_simulation_start -# print("{}{} {:.1f}s ".format('█' * round(n_steps), '-' * round(0), time_simulation_consumed)) -# break -# else: -# path_ = q.qsize() / n_simulations * n_steps -# print("{}{} {:03.1f}%".format('█' * int(round(path_)), '-' * int(n_steps - round(path_)), -# path_ / n_steps * 100), -# end='\r') -# time.sleep(1) -# p.close() -# p.join() -# results = jobs.get() diff --git a/sfeprapy/mcs2/__init__.py b/sfeprapy/mcs2/__init__.py new file mode 100644 index 0000000..74b2f1c --- /dev/null +++ b/sfeprapy/mcs2/__init__.py @@ -0,0 +1,83 @@ +# -*- coding: utf-8 -*- +import pandas as pd + +from sfeprapy.func.mcs_gen import dict_flatten + + +def __example_config_dict(): + return dict(n_threads=2, cwd='') + + +def __example_input_dict(): + y = { + "Standard Case 1": dict( + case_name="Standard Case 1", + n_simulations=5000, + fire_time_step=30, + fire_time_duration=18000, + fire_hrr_density=0.25, + fire_load_density=dict(dist="gumbel_r_", lbound=10, ubound=1200, mean=420, sd=126), + fire_spread_speed=dict(dist="uniform_", lbound=0.0035, ubound=0.0190), + fire_nft_limit=dict(dist="norm_", lbound=623.15, ubound=2023.15, mean=1323.15, sd=93), + fire_combustion_efficiency=1.0, + window_open_fraction=1.0, + phi_teq=dict(dist="constant_", ubound=1, lbound=1, mean=0, sd=0), + beam_cross_section_area=0.017, + beam_position_horizontal=-1, + beam_position_vertical=3.2, + beam_rho=7850, + fire_mode=0, + fire_gamma_fi_q=1, + fire_t_alpha=300, + fire_tlim=0.333, + protection_c=1700, + protection_k=0.2, + protection_protected_perimeter=2.14, + protection_rho=800, + room_floor_area=500, + room_breadth_depth_ratio=dict(dist="uniform_", lbound=0.512 - 0.2, ubound=0.512 + 0.2), + room_height=3, + room_wall_thermal_inertia=720, + solver_temperature_goal=893.15, + solver_max_iter=20, + solver_thickness_lbound=0.0001, + solver_thickness_ubound=0.0500, + solver_tol=1.0, + window_height=2.8, + window_floor_ratio=dict(dist="uniform_", lbound=0.05, ubound=0.4), + window_open_fraction_permanent=0, + timber_exposed_area=0, + timber_charring_rate=0.7, # mm/min + timber_hc=13.2, # MJ/kg + timber_density=400, # [kg/m3] + timber_solver_ilim=20, + timber_solver_tol=1, + ), + } + return y + + +def __example_input_csv(x: dict): + y = {k: dict_flatten(v) for k, v in x.items()} + y = pd.DataFrame.from_dict(y, orient="columns") + y.index.name = "PARAMETERS" + y = y.to_csv(index=True, line_terminator='\n') + return y + + +def __example_input_df(x: dict) -> pd.DataFrame: + y = {k: dict_flatten(v) for k, v in x.items()} + y = pd.DataFrame.from_dict(y, orient="columns") + y.index.name = "PARAMETERS" + return y + + +EXAMPLE_CONFIG_DICT = __example_config_dict() +EXAMPLE_INPUT_DICT = __example_input_dict() +EXAMPLE_INPUT_CSV = __example_input_csv(__example_input_dict()) +EXAMPLE_INPUT_DF = __example_input_df(__example_input_dict()) + +if __name__ == "__main__": + print(EXAMPLE_CONFIG_DICT, "\n") + print(EXAMPLE_INPUT_DICT, "\n") + print(EXAMPLE_INPUT_CSV, "\n") diff --git a/sfeprapy/mcs2/__main__.py b/sfeprapy/mcs2/__main__.py new file mode 100644 index 0000000..4cd3be8 --- /dev/null +++ b/sfeprapy/mcs2/__main__.py @@ -0,0 +1,25 @@ +import os + +from sfeprapy.mcs2.mcs2_calc import MCS2 + + +def main(fp_mcs_in: str, n_threads: int = None): + fp_mcs_in = os.path.realpath(fp_mcs_in) + + mcs = MCS2() + try: + mcs.mcs_inputs = fp_mcs_in + except Exception as e: + raise e + + try: + if n_threads: + mcs.mcs_config = ( + dict(n_threads=n_threads) + if mcs.mcs_config + else mcs.mcs_config.update(dict(n_threads=n_threads)) + ) + except KeyError: + pass + + mcs.run_mcs() diff --git a/sfeprapy/mcs2/mcs2_calc.py b/sfeprapy/mcs2/mcs2_calc.py new file mode 100644 index 0000000..14dc944 --- /dev/null +++ b/sfeprapy/mcs2/mcs2_calc.py @@ -0,0 +1,152 @@ +# -*- coding: utf-8 -*- +from sfeprapy.mcs0.mcs0_calc import MCS0 +from sfeprapy.mcs0.mcs0_calc import teq_main as teq_main_mcs0 + + +def teq_main_wrapper(args): + try: + kwargs, q = args + q.put("index: {}".format(kwargs["index"])) + return teq_main(**kwargs) + except (ValueError, AttributeError): + return teq_main(**args) + + +def teq_main( + case_name: str, + n_simulations: int, + probability_weight: float, + index: int, + beam_cross_section_area: float, + beam_position_vertical: float, + beam_position_horizontal: float, + beam_rho: float, + fire_time_duration: float, + fire_time_step: float, + fire_combustion_efficiency: float, + fire_gamma_fi_q: float, + fire_hrr_density: float, + fire_load_density: float, + fire_mode: int, + fire_nft_limit: float, + fire_spread_speed: float, + fire_t_alpha: float, + fire_tlim: float, + protection_c: float, + protection_k: float, + protection_protected_perimeter: float, + protection_rho: float, + # room_breadth: float, + # room_depth: float, + room_wall_thermal_inertia: float, + room_floor_area: float, + room_breadth_depth_ratio: float, + room_height: float, + solver_temperature_goal: float, + solver_max_iter: int, + solver_thickness_lbound: float, + solver_thickness_ubound: float, + solver_tol: float, + # window_width: float, + window_height: float, + window_floor_ratio: float, + window_open_fraction: float, + window_open_fraction_permanent: float, + phi_teq: float = 1.0, + timber_charring_rate=None, + timber_hc: float = None, + timber_density: float = None, + timber_exposed_area: float = None, + timber_solver_tol: float = None, + timber_solver_ilim: float = None, + *_, + **__, +) -> dict: + kwargs = locals() + _ = kwargs.pop('_') + __ = kwargs.pop('__') + + kwargs.pop('room_floor_area') + kwargs.pop('room_breadth_depth_ratio') + + # ----------------------------------------- + # Calculate `room_breadth` and `room_depth` + # ----------------------------------------- + # room_depth * room_breadth = room_floor_area + # room_breadth / room_depth = room_breadth_depth_ratio + + # room_breadth = room_breadth_depth_ratio * room_depth + # room_depth * room_breadth_depth_ratio * room_depth = room_floor_area + room_depth = (room_floor_area / room_breadth_depth_ratio) ** 0.5 + room_breadth = room_breadth_depth_ratio * (room_floor_area / room_breadth_depth_ratio) ** 0.5 + assert room_breadth_depth_ratio <= 1 + assert abs(room_depth * room_breadth - room_floor_area) < 1e-5 + + # ------------------------------ + # Calculate window opening width + # ------------------------------ + window_width = room_floor_area * window_floor_ratio / window_height + + # ---------------------------------- + # Calculate beam horizontal location + # ---------------------------------- + + kwargs.update(dict( + room_breadth=room_depth, + room_depth=room_breadth, + window_width=window_width, + beam_horizontal_location=0.8 * room_depth + )) + + outputs = teq_main_mcs0(**kwargs) + + return outputs + + +class MCS2(MCS0): + def __init__(self): + super().__init__() + + def mcs_deterministic_calc(self, *args, **kwargs) -> dict: + return teq_main(*args, **kwargs) + + def mcs_deterministic_calc_mp(self, *args, **kwargs) -> dict: + return teq_main_wrapper(*args, **kwargs) + + +def _test_standard_case_new(): + import copy + from sfeprapy.mcs2 import EXAMPLE_INPUT_DICT, EXAMPLE_CONFIG_DICT + from scipy.interpolate import interp1d + import numpy as np + + # increase the number of simulations so it gives sensible results + mcs_input = copy.deepcopy(EXAMPLE_INPUT_DICT) + mcs_config = copy.deepcopy(EXAMPLE_CONFIG_DICT) + for k in list(mcs_input.keys()): + mcs_input[k]["phi_teq"] = 1 + mcs_input[k]["n_simulations"] = 10000 + mcs_input[k]["probability_weight"] = 1 / 3.0 + mcs_input[k]["fire_time_duration"] = 10000 + mcs_input[k]["timber_exposed_area"] = 0 + + # increase the number of threads so it runs faster + mcs_config["n_threads"] = 1 # coverage does not support + mcs2 = MCS2() + mcs2.mcs_inputs = mcs_input + mcs2.mcs_config = mcs_config + mcs2.run_mcs() + mcs_out = mcs2.mcs_out + teq = mcs_out["solver_time_equivalence_solved"] / 60.0 + hist, edges = np.histogram(teq, bins=np.arange(0, 181, 0.5)) + x, y = (edges[:-1] + edges[1:]) / 2, np.cumsum(hist / np.sum(hist)) + teq_at_80_percentile = interp1d(y, x)(0.8) + print(teq_at_80_percentile) + # target, target_tol = 60, 2 + # assert target - target_tol < teq_at_80_percentile < target + target_tol + + +if __name__ == '__main__': + # _test_teq_phi() + # _test_standard_case() + _test_standard_case_new() diff --git a/test/mc_figures.py b/test/mc_figures.py index b31854c..60c6682 100644 --- a/test/mc_figures.py +++ b/test/mc_figures.py @@ -1,24 +1,22 @@ # -*- coding: utf-8 -*- import numpy as np +from fsetools.lib.fse_bs_en_1991_1_2_parametric_fire import temperature as _fire_param +from fsetools.lib.fse_bs_en_1993_1_2_heat_transfer_c import temperature as _steel_heat_transfer +from fsetools.lib.fse_travelling_fire import temperature_si as _fire_travelling from matplotlib import pyplot as plt -from sfeprapy.dat.steel_carbon import steel_specific_heat_carbon_steel from sfeprapy.func.fire_iso834 import fire as _fire_iso -from sfeprapy.func.fire_parametric_ec import fire as _fire_param -from sfeprapy.func.fire_travelling import fire as _fire_travelling -from sfeprapy.func.heat_transfer_protected_steel_ec import protected_steel_eurocode as _steel_heat_transfer def fire_param(ax, label=None): - - x = np.arange(0, 60*180, 1, dtype=float) + x = np.arange(0, 60 * 180, 1, dtype=float) # PARAMETRIC FIRE inputs = { "t": x, - "A_t": (10*40+40*2+2*10)*2, - "A_f": 10*40, - "A_v": 5*5, + "A_t": (10 * 40 + 40 * 2 + 2 * 10) * 2, + "A_f": 10 * 40, + "A_v": 5 * 5, "h_eq": 5, "q_fd": 420e6, "lambda_": 1, @@ -26,34 +24,33 @@ def fire_param(ax, label=None): "c": 2250000, "t_lim": 20 * 60, "temperature_initial": 273.15 - } + } y = _fire_param(**inputs) - ax.plot(x/60., y-273.15, label=label) + ax.plot(x / 60., y - 273.15, label=label) def fire_travel(ax, label=None): - x = np.arange(0, 60*180, 1, dtype=float) + x = np.arange(0, 60 * 180, 1, dtype=float) inputs = { "t": x, "T_0": 273.15, "q_fd": 420e6, - "RHRf": 0.15e6, + "hrrpua": 0.15e6, "l": 40, "w": 10, "s": 0.01, - "h_s": 4.5, - "l_s": 25, + "e_h": 4.5, + "e_l": 25, } y = _fire_travelling(**inputs) - ax.plot(x/60., y-273.15, label=label) + ax.plot(x / 60., y - 273.15, label=label) def heat_transfer_param(ax, label=None): - - x = np.arange(0, 60*180, 1, dtype=float) + x = np.arange(0, 60 * 180, 1, dtype=float) # PARAMETRIC FIRE inputs_parametric_fire = { @@ -68,77 +65,66 @@ def heat_transfer_param(ax, label=None): "c": 2250000, "t_lim": 20 * 60, "temperature_initial": 273.15 - } + } y_ = _fire_param(**inputs_parametric_fire) - # steel_prop = Thermal() - - _, y, _ = _steel_heat_transfer( - time=x, - temperature_ambient=y_, - rho_steel=7850, - # c_steel_T=steel_specific_heat_carbon_steel, - area_steel_section=0.017, - k_protection=0.2, - rho_protection=800, - c_protection=1700, - thickness_protection=0.00938, - perimeter_protected=2.14, - # is_terminate_peak=False + y = _steel_heat_transfer( + fire_time=x, + fire_temperature=y_, + beam_rho=7850, + beam_cross_section_area=0.017, + protection_k=0.2, + protection_rho=800, + protection_c=1700, + protection_thickness=0.00938, + protection_protected_perimeter=2.14, ) - ax.plot(x/60., y-273.15, label=label) + ax.plot(x / 60., y - 273.15, label=label) def heat_transfer_travel(ax, label=None): - - x = np.arange(0, 60*180, 1, dtype=float) + x = np.arange(0, 60 * 180, 1, dtype=float) inputs = { "t": x, "T_0": 273.15, "q_fd": 420e6, - "RHRf": 0.15e6, + "hrrpua": 0.15e6, "l": 40, "w": 10, "s": 0.01, - "h_s": 4.5, - "l_s": 25, + "e_h": 4.5, + "e_l": 25, } y_ = _fire_travelling(**inputs) - # steel_prop = Thermal() - - _, y, _ = _steel_heat_transfer( - time=x, - temperature_ambient=y_, - rho_steel=7850, - # c_steel_T=steel_prop.c(), - area_steel_section=0.017, - k_protection=0.2, - rho_protection=800, - c_protection=1700, - thickness_protection=0.00938, - perimeter_protected=2.14, - # is_terminate_peak=False + y = _steel_heat_transfer( + fire_time=x, + fire_temperature=y_, + beam_rho=7850, + beam_cross_section_area=0.017, + protection_k=0.2, + protection_rho=800, + protection_c=1700, + protection_thickness=0.00938, + protection_protected_perimeter=2.14, ) - ax.plot(x/60., y-273.15, label=label) + ax.plot(x / 60., y - 273.15, label=label) def fire_iso834(ax, label=None): - - x = np.arange(0, 60*180, 1) + x = np.arange(0, 60 * 180, 1) y = _fire_iso(x, 20 + 273.15) - ax.plot(x/60., y-273.15, label=label) + ax.plot(x / 60., y - 273.15, label=label) def plot_examples(list_func, list_linelabel, figname, xlabel, ylabel): - plt.style.use('seaborn-paper') - fig, ax = plt.subplots(figsize=(3.94*1.5, 2.76*1.5)) # (3.94, 2.76) for large and (2.5, 2) for small figure size + fig, ax = plt.subplots(figsize=(3.94 * 1.5, 2.76 * 1.5)) # (3.94, 2.76) for large and (2.5, 2) for small figure size for i, func in enumerate(list_func): func(ax, label=list_linelabel[i]) ax.set_xlabel(xlabel) @@ -157,11 +143,6 @@ def plot_examples(list_func, list_linelabel, figname, xlabel, ylabel): if __name__ == '__main__': plt.style.use('seaborn-paper') - # plot_examples( - # [fire_travel], [None], - # 'fire_travel.png', 'Time [minute]', 'Temperature [$^{\circ}C$]' - # ) - plot_examples( [fire_travel, heat_transfer_travel], ["Gas temperature", "Steel temperature"], 'heat_transfer_travel.png', 'Time [minute]', 'Temperature [$^{\circ}C$]' @@ -186,8 +167,3 @@ def plot_examples(list_func, list_linelabel, figname, xlabel, ylabel): [fire_iso834], [None], 'fire_iso.png', 'Time [minute]', 'Temperature [$^{\circ}C$]' ) - - # plot_examples( - # [fire_iso834, heattran], ["Gas temperature", "Steel temperature"], - # 'heat_transfer_param.png', 'Time [minute]', 'Temperature [$^{\circ}C$]' - # ) diff --git a/test/test_ec_3_1_2_kyT_probability.py b/test/test_ec_3_1_2_kyT_probability.py deleted file mode 100644 index 5aa6c01..0000000 --- a/test/test_ec_3_1_2_kyT_probability.py +++ /dev/null @@ -1,3 +0,0 @@ -from sfeprapy.dat.ec_3_1_2_kyT_probability import _test_probabilistic as test_probability - -test_probability() diff --git a/test/test_func_fire_parametric_ec.py b/test/test_func_fire_parametric_ec.py deleted file mode 100644 index 78a033f..0000000 --- a/test/test_func_fire_parametric_ec.py +++ /dev/null @@ -1,6 +0,0 @@ -# -*- coding: utf-8 -*- - - -from sfeprapy.func.fire_parametric_ec import _test_fire as test_fire_parametric_ec - -test_fire_parametric_ec() diff --git a/test/test_func_travelling.py b/test/test_func_travelling.py deleted file mode 100644 index 72bab91..0000000 --- a/test/test_func_travelling.py +++ /dev/null @@ -1,17 +0,0 @@ -# -*- coding: utf-8 -*- - - -from sfeprapy.func.fire_travelling import _test_fire as test_fire_travelling -from sfeprapy.func.fire_travelling import ( - _test_fire_backup as test_fire_travelling_backup, -) -from sfeprapy.func.fire_travelling import ( - _test_fire_multiple_beam_location as test_fire_travelling_multiple, -) -from sfeprapy.func.fire_travelling_flux import _test_fire as test_fire_travelling_flux - -test_fire_travelling() -test_fire_travelling_backup() -test_fire_travelling_multiple() -test_fire_travelling_flux() - diff --git a/test/test_mcs0.py b/test/test_mcs0.py index 815797b..40a32e6 100644 --- a/test/test_mcs0.py +++ b/test/test_mcs0.py @@ -1,9 +1,7 @@ # -*- coding: utf-8 -*- +from sfeprapy.mcs0.mcs0_calc import _test_standard_case as test_standard_case from sfeprapy.mcs0.mcs0_calc import _test_teq_phi as test_teq_phi -from sfeprapy.mcs0.mcs0_calc import _test_standard_case as test_mcs0_test_standard_case -from sfeprapy.mcs0.mcs0_calc import _test_standard_oversized_case as test_mcs0_test_standard_oversized_case test_teq_phi() -test_mcs0_test_standard_case() -test_mcs0_test_standard_oversized_case() \ No newline at end of file +test_standard_case()