diff --git a/requirements.txt b/requirements.txt index ca83735..d9e6802 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,10 +1,10 @@ -numpy>=1.17 -pandas>=0.25 -scipy>=1.3 -seaborn>=0.9 +numpy>=1.17.1 +pandas>=0.25.1 +scipy>=1.3.1 +matplotlib>=3.1.1 tqdm>=4.36 -xlrd>=1.2 -plotly>=4.2 +xlrd>=1.2.0 +plotly>=4.2.1 docopt>=0.6 pyside2>=5.12 -openpyxl>=3.0.3 \ No newline at end of file +openpyxl>=3.0.3 diff --git a/setup.py b/setup.py index bbe7f14..981e1b2 100644 --- a/setup.py +++ b/setup.py @@ -40,10 +40,12 @@ "sfeprapy.mcs1", "sfeprapy.func", "sfeprapy.dist_fit", - "sfeprapy.guilayout", - "sfeprapy.guilogic", + "sfeprapy.cli", + "sfeprapy.gui", + "sfeprapy.gui.layout", + "sfeprapy.gui.logic", ], install_requires=requirements, include_package_data=True, - entry_points={"console_scripts": ["sfeprapy=sfeprapy.cli:main"]}, + entry_points={"console_scripts": ["sfeprapy=sfeprapy.cli.__main__:main"]}, ) diff --git a/sfeprapy/__init__.py b/sfeprapy/__init__.py index 4ba4973..2199384 100644 --- a/sfeprapy/__init__.py +++ b/sfeprapy/__init__.py @@ -34,7 +34,7 @@ """ -__version__ = "0.7.2" +__version__ = "0.7.1.post6" def check_pip_upgrade(): diff --git a/sfeprapy/cli/__init__.py b/sfeprapy/cli/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/sfeprapy/cli.py b/sfeprapy/cli/__main__.py similarity index 84% rename from sfeprapy/cli.py rename to sfeprapy/cli/__main__.py index b75e2b0..7d99724 100644 --- a/sfeprapy/cli.py +++ b/sfeprapy/cli/__main__.py @@ -1,5 +1,6 @@ """SfePrapy CLI Help. Usage: + sfeprapy sfeprapy mcs0 gui sfeprapy mcs0 run [-p=] sfeprapy mcs0 figure @@ -18,7 +19,8 @@ -h --help to show this message. Commands: - mcs0 gui A GUI version of sfeprapy. + 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 . @@ -26,6 +28,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.__main__ import main as mcs0 +from sfeprapy.mcs0.__main__ import save_figure as mcs0_figure + def main(): import os @@ -38,12 +46,10 @@ def main(): if arguments["mcs0"]: if arguments["gui"]: - from sfeprapy.guilogic.main import main as main_ - main_() + main_gui() return 0 elif arguments["figure"]: - from sfeprapy.mcs0.__main__ import save_figure as mcs0_figure mcs0_figure(fp_mcs0_out=arguments[""]) @@ -57,7 +63,6 @@ def main(): f.write(EXAMPLE_INPUT_CSV) else: - from sfeprapy.mcs0.__main__ import main as mcs0 fp_mcs_in = arguments[""] n_threads = arguments["-p"] or 2 @@ -65,9 +70,6 @@ def main(): elif arguments["distfit"]: - # Prerequisites - from sfeprapy.func.stats_dist_fit import auto_fit_2 - # Default values data_type = arguments["--data_t"] or 2 distribution_list = arguments["--dist_g"] or 1 @@ -80,4 +82,4 @@ def main(): ) else: - print("instruction unclear.") + main_gui() diff --git a/sfeprapy/func/asciiplot.py b/sfeprapy/func/asciiplot.py new file mode 100644 index 0000000..9a6f3a1 --- /dev/null +++ b/sfeprapy/func/asciiplot.py @@ -0,0 +1,273 @@ +# -*- coding: utf-8 -*- + +from typing import Union, Tuple + +import numpy as np + + +class AsciiPlot: + + def __init__(self, size: Tuple[float, float] = (80, 35)): + self.__size = None + self.__xlim = None + self.__ylim = None + self.__plot_canvas = None + self.__plot_yaxis = None + self.__plot_xaxis = None + self.__plot = None + + self.size = size + + def plot(self, x, y, xlim: tuple = None, ylim: tuple = None): + + if xlim is None: + xlim = [min(x), max(x)] + if ylim is None: + ylim = [min(y), max(y)] + if ylim[0] == ylim[1] == 0: + ylim = -1., 1. + elif ylim[0] == ylim[1]: + ylim[0] -= ylim[0] * 0.1 + ylim[1] += ylim[1] * 0.11 + + yaxis_size = (None, self.size[1] - 2) + self.__plot_yaxis = self.__make_yaxis(ylim, yaxis_size) + xaxis_size = (self.size[0] - self.__plot_yaxis.shape[1], None) + self.__plot_xaxis = self.__make_xaxis(xlim, xaxis_size) + canvas_size = (self.size[0] - self.__plot_yaxis.shape[1], self.size[1] - self.__plot_xaxis.shape[0]) + self.__plot_canvas = self._make_canvas(x, y, xlim, ylim, canvas_size) + self.__plot = self.__assemble(self.__plot_canvas, self.__plot_xaxis, self.__plot_yaxis) + + return self + + def str(self): + return self.__list2str(self.__plot) + + def show(self): + try: + assert self.__plot is not None + except AssertionError: + ValueError('No plot to show') + + print(self.__list2str(self.__plot)) + + def save(self, fp: str): + try: + assert self.__plot is not None + except AssertionError: + ValueError('No plot to save') + + with open(fp, 'w+') as f: + f.write(self.__list2str(self.__plot)) + + @staticmethod + def __list2str(v: Union[list, np.ndarray]) -> str: + return '\n'.join([''.join([str(j) for j in i]) for i in v]) + + @staticmethod + def _make_canvas( + x: Union[list, np.ndarray], + y: Union[list, np.ndarray], + xlim: tuple = None, + ylim: tuple = None, + size: tuple = (80, 40), + mat_plot_canvas: np.ndarray = None + ) -> np.ndarray: + # correct data type if necessary + if not isinstance(x, np.ndarray): + x = np.array(x) + if not isinstance(y, np.ndarray): + y = np.array(y) + if mat_plot_canvas is not None: + assert mat_plot_canvas.shape[0] == size[::-1][0] + assert mat_plot_canvas.shape[1] == size[::-1][1] + else: + mat_plot_canvas = np.zeros(shape=size[::-1]) # canvas matrix containing the plot + + # workout data plot boundaries + if xlim: + x1, x2 = xlim + else: + x1, x2 = min(x), max(x) + if ylim: + y1, y2 = ylim + else: + y1, y2 = min(y), max(y) + + # workout translator arrays + # example + # mat_plot_canvas[i, j] equates location ii[i], jj[i] in data coordination system + ii = np.linspace(y1, y2, size[1]) + jj = np.linspace(x1, x2, size[0]) + dx, dy = jj[1] - jj[0], ii[1] - ii[0] + + # interpolate x and y, i.e. to increase resolution of the provided x, y + xi = np.arange(min(x), max(x) + dx / 2, dx / 2) + yi = np.interp(xi, x, y) + + # plot, assign value 1 to cells in `mat_plot_canvas` the x, y line lands + for i in range(len(xi)): + x_, y_ = xi[i], yi[i] + + # cells that the line lands + jjp = np.argwhere(np.logical_and(jj > x_ - 0.5 * dx, jj < x_ + 0.5 * dx)) + iip = np.argwhere(np.logical_and(ii > y_ - 0.5 * dy, ii < y_ + 0.5 * dy)) + + # assign value 1 to the cells + for ip in iip: + for jp in jjp: + mat_plot_canvas[ip, jp] = 1 + + return mat_plot_canvas + + @staticmethod + def __make_yaxis(ylim, size, label_fmt: str = None) -> np.ndarray: + + # workout data plot boundaries + y1, y2 = ylim + + # workout width and height per cell of the canvas matrix + ii = np.linspace(y1, y2, size[1]) + + if label_fmt is None: + max_digits = 0 + for i in ii: + if len(str(int(i))) > max_digits: + max_digits = len(str(int(i))) + length = max_digits + 4 # (potential) minus sign + max_length + decimal symbol + 2 decimal places + label_fmt = '{:' + str(int(length)) + '.2f}' + else: + length = len(label_fmt.format(0.1)) + + mat_yaxis = np.full(shape=(size[1], length + 1), fill_value=' ', dtype=' np.ndarray: + + jj = np.linspace(min(xlim), max(xlim), size[0]) + + if label_fmt is None: + max_digits = 0 + for j in jj: + if len(str(int(j))) > max_digits: + max_digits = len(str(int(j))) + length = max_digits + 4 # (potential) minus sign + max_length + decimal symbol + 2 decimal places + label_fmt = '{:' + str(int(length)) + '.2f}' + else: + length = len(label_fmt.format(0.1)) + list_x_tick_labels_exhaustive = [label_fmt.format(j) for j in jj] + + dx_label = max([len(i) for i in list_x_tick_labels_exhaustive]) + 2 + + mat_xaxis = np.full(shape=(2, size[0]), fill_value=' ', dtype=' np.ndarray: + mat_canvas = np.full_like(mat_plot_canvas, fill_value=' ', dtype=' 0] = '*' + + # construct figure matrix + mat_figure = np.full( + # shape=(mat_plot_canvas.shape[0] + mat_xaxis.shape[0], mat_yaxis.shape[1] + mat_plot_canvas.shape[1]), + shape=(mat_plot_canvas.shape[0] + mat_xaxis.shape[0], mat_plot_canvas.shape[1] + mat_yaxis.shape[1]), + fill_value=' ', dtype=' Tuple[float, float]: + return self.__size + + @size.setter + def size(self, v: Tuple[float, float]): + assert any([isinstance(v, tuple), isinstance(v, list)]) + assert len(v) == 2 + self.__size = v + + @property + def xlim(self): + return self.__xlim + + @xlim.setter + def xlim(self, v: Tuple[float, float]): + assert any([isinstance(v, tuple), isinstance(v, list)]) + assert len(v) == 2 + + self.__xlim = v + + @property + def ylim(self): + return self.__ylim + + @ylim.setter + def ylim(self, v: Tuple[float, float]): + try: + assert all([isinstance(v, tuple), isinstance(v, list), len(v) == 2]) + except: + raise ValueError('`ylim` should be a tuple with length equal to 2') + + self.__ylim = v + + +def _test_asciiplot(): + size = (80, 25) + xlim = (-2 * np.pi, 2 * np.pi) + ylim = (-1.1, 1.1) + + aplot = AsciiPlot(size=size) + + x = np.linspace(-2 * np.pi, 2 * np.pi, 50) + aplot.plot( + x=x, + y=np.sin(x), + xlim=xlim, + ylim=ylim, + ).show() + + +if __name__ == '__main__': + _test_asciiplot() diff --git a/sfeprapy/func/fire_travelling_flux.py b/sfeprapy/func/fire_travelling_flux.py new file mode 100644 index 0000000..110b822 --- /dev/null +++ b/sfeprapy/func/fire_travelling_flux.py @@ -0,0 +1,207 @@ +# -*- 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 index dda02c2..9458b2e 100644 --- a/sfeprapy/func/heat_transfer_protected_steel_ec.py +++ b/sfeprapy/func/heat_transfer_protected_steel_ec.py @@ -163,6 +163,7 @@ def protected_steel_eurocode_max_temperature( 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] diff --git a/sfeprapy/func/mcs.py b/sfeprapy/func/mcs.py new file mode 100644 index 0000000..75890f2 --- /dev/null +++ b/sfeprapy/func/mcs.py @@ -0,0 +1,244 @@ +import copy +import os +import time +from typing import Union, Callable + +import pandas as pd +from tqdm import tqdm + +from sfeprapy.func.mcs_gen import main as mcs_gen_main + + +class MCS: + """ + 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. + + To help to understand this class, a brief process of a MCS are outlined below: + + 1. stochastic problem definition (i.e. raw input parameters from the user). + 2. sample deterministic parameters from the stochastic inputs + 3. iterate the sampled parameters and run deterministic calculation + 4. summarise results + + Following above, this object consists of four primary methods, each matching one of the above step: + + `MCS.mcs_inputs` + A method to intake user inputs, restructure data so usable within MCS. + `MCS.mcs_sampler` + A method to sample deterministic parameters from the stochastic parameters, produces input to be used in + the next step. + `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` + A method to post processing results. + NOTE! This method needs to be re-defined in a child class. + """ + DEFAULT_TEMP_FOLDER_NAME = "mcs.out" + DEFAULT_MCS_OUTPUT_FILE_NAME = "mcs.out.csv" + DEFAULT_CONFIG_FILE_NAME = "config.json" + DEFAULT_CONFIG = dict(n_threads=1) + + def __init__(self): + + # ------------------------------ + # instantiate internal variables + # ------------------------------ + 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 + self._func_mcs_calc: Callable = None # monte carlo simulation deterministic calculation routine + self._func_mcs_calc_mp: Callable = None # multiprocessing version of `MCS._mcs_calc` + self.__mcs_post: Callable = None + self.__mcs_out: pd.DataFrame = None + + # assign default properties + self.func_mcs_gen = mcs_gen_main + + def mcs_deterministic_calc(self, *args, **kwargs) -> dict: + """Placeholder method. The Monte Carlo Simulation deterministic calculation routine. + :return: + """ + pass + + def mcs_deterministic_calc_mp(self, *args, **kwargs) -> dict: + """Placeholder method. The Monte Carlo Simulation deterministic calculation routine. + :return: + """ + pass + + @property + def mcs_inputs(self) -> dict: + return self.__mcs_inputs + + @mcs_inputs.setter + def mcs_inputs(self, fp_df_dict: Union[str, pd.DataFrame, dict]): + """ + :param fp_df_dict: + :return: + """ + if isinstance(fp_df_dict, str): + fp = fp_df_dict + self.path_wd = 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"): + self.__mcs_inputs = pd.read_csv(fp).set_index("PARAMETERS").to_dict() + else: + raise ValueError("Unknown input file format.") + elif isinstance(fp_df_dict, pd.DataFrame): + self.__mcs_inputs = fp_df_dict.to_dict() + elif isinstance(fp_df_dict, dict): + self.__mcs_inputs = fp_df_dict + else: + raise TypeError("Unknown input data type.") + + @property + def mcs_config(self) -> dict: + """simulation configuration""" + if self.__mcs_config is not None: + return self.__mcs_config + else: + return self.DEFAULT_CONFIG + + @mcs_config.setter + def mcs_config(self, config): + self.__mcs_config = config + if 'cwd' in config: + self.__cwd = config['cwd'] + + def run_mcs(self): + # ---------------------------- + # Prepare mcs parameter inputs + # ---------------------------- + # elevate dict dimension, for example + # { + # 'fuel_load:lbound': 100, + # 'fuel_load:ubound': 200 + # } + # became + # { + # 'fuel_load': { + # 'lbound': 100, + # 'ubound': 200 + # } + # } + x1 = self.mcs_inputs + for case_name in list(x1.keys()): + for param_name in list(x1[case_name].keys()): + if ":" in param_name: + param_name_parent, param_name_sibling = param_name.split(":") + if param_name_parent not in x1[case_name]: + x1[case_name][param_name_parent] = dict() + x1[case_name][param_name_parent][param_name_sibling] = x1[case_name].pop(param_name) + # to convert all "string numbers" to float data type + for case_name in x1.keys(): + for i, v in x1[case_name].items(): + if isinstance(v, dict): + for ii, vv in v.items(): + try: + x1[case_name][i][ii] = float(vv) + except: + pass + else: + try: + x1[case_name][i] = float(v) + except: + pass + + # ------------------------------ + # Generate mcs parameter samples + # ------------------------------ + x2 = {k: self.func_mcs_gen(v, int(v["n_simulations"])) for k, v in x1.items()} + + # ------------------ + # Run mcs simulation + # ------------------ + x3 = dict() # output container, structure: + # { + # 'case_1': df1, + # 'case_2': df2 + # } + for k, v in x2.items(): + + x3_ = self.__mcs_mp( + self.mcs_deterministic_calc, + self.mcs_deterministic_calc_mp, + x=v, + n_threads=self.mcs_config["n_threads"], + ) + if self.mcs_post: + mcs_post = self.mcs_post + x3_ = mcs_post(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", + ) + ) + + # ------------ + # 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, + ) + + def mcs_post(self, *arg, **kwargs): + pass + + @property + def mcs_out(self): + return self.__mcs_out + + @staticmethod + def __mcs_mp(func, func_mp, x: pd.DataFrame, n_threads: int) -> pd.DataFrame: + list_mcs_in = x.to_dict(orient="records") + + time.sleep(0.5) # to avoid clashes between the prints and progress bar + print("{:<24.24}: {}".format("CASE", list_mcs_in[0]["case_name"])) + print("{:<24.24}: {}".format("NO. OF THREADS", n_threads)) + print("{:<24.24}: {}".format("NO. OF SIMULATIONS", len(x.index))) + time.sleep(0.5) # to avoid clashes between the prints and progress bar + + if n_threads == 1 or func_mp is None: + mcs_out = list() + for i in tqdm(list_mcs_in, ncols=60): + mcs_out.append(func(**i)) + 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) + break + else: + if q.qsize() - pbar.n > 0: + pbar.update(q.qsize() - pbar.n) + time.sleep(1) + p.close() + p.join() + mcs_out = jobs.get() + time.sleep(0.5) + + # clean and convert results to dataframe and return + df_mcs_out = pd.DataFrame(mcs_out) + df_mcs_out.sort_values("solver_time_equivalence_solved", inplace=True) # sort base on time equivalence + return df_mcs_out diff --git a/sfeprapy/gui/__init__.py b/sfeprapy/gui/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/sfeprapy/gui/__main__.py b/sfeprapy/gui/__main__.py new file mode 100644 index 0000000..8106faf --- /dev/null +++ b/sfeprapy/gui/__main__.py @@ -0,0 +1,22 @@ + +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 new file mode 100644 index 0000000..e69de29 diff --git a/sfeprapy/guilayout/main.py b/sfeprapy/gui/layout/main.py similarity index 82% rename from sfeprapy/guilayout/main.py rename to sfeprapy/gui/layout/main.py index 674413e..7bff0e0 100644 --- a/sfeprapy/guilayout/main.py +++ b/sfeprapy/gui/layout/main.py @@ -20,7 +20,7 @@ class Ui_MainWindow(object): def setupUi(self, MainWindow): if MainWindow.objectName(): MainWindow.setObjectName(u"MainWindow") - MainWindow.resize(426, 178) + MainWindow.resize(446, 186) sizePolicy = QSizePolicy(QSizePolicy.Minimum, QSizePolicy.Minimum) sizePolicy.setHorizontalStretch(0) sizePolicy.setVerticalStretch(0) @@ -31,37 +31,37 @@ def setupUi(self, MainWindow): self.centralwidget.setObjectName(u"centralwidget") self.pushButton_run = QPushButton(self.centralwidget) self.pushButton_run.setObjectName(u"pushButton_run") - self.pushButton_run.setGeometry(QRect(305, 95, 112, 32)) + 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(5, 95, 112, 32)) + 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(5, 5, 416, 81)) + 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, 20, 112, 32)) + 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, 58, 16)) + 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, 21)) + 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, 50, 156, 21)) + 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, 50, 116, 16)) - self.line = QFrame(self.centralwidget) - self.line.setObjectName(u"line") - self.line.setGeometry(QRect(0, 85, 426, 16)) - self.line.setFrameShape(QFrame.HLine) - self.line.setFrameShadow(QFrame.Sunken) + 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, 426, 22)) + self.menubar.setGeometry(QRect(0, 0, 446, 22)) MainWindow.setMenuBar(self.menubar) self.statusbar = QStatusBar(MainWindow) self.statusbar.setObjectName(u"statusbar") @@ -82,7 +82,7 @@ def retranslateUi(self, MainWindow): 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"GroupBox", 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)) diff --git a/sfeprapy/guilayout/main.ui b/sfeprapy/gui/layout/ui/main.ui similarity index 77% rename from sfeprapy/guilayout/main.ui rename to sfeprapy/gui/layout/ui/main.ui index 4b8954a..0d618a9 100644 --- a/sfeprapy/guilayout/main.ui +++ b/sfeprapy/gui/layout/ui/main.ui @@ -6,8 +6,8 @@ 0 0 - 426 - 178 + 446 + 186 @@ -29,12 +29,17 @@ - 305 - 95 - 112 - 32 + 335 + 130 + 96 + 26 + + + 9 + + Run simulation @@ -45,12 +50,17 @@ - 5 - 95 - 112 - 32 + 230 + 130 + 96 + 26 + + + 9 + + Make and save a template input file @@ -61,22 +71,27 @@ - 5 - 5 + 15 + 15 416 - 81 + 101 + + + 9 + + - GroupBox + Inputs 300 - 20 - 112 - 32 + 25 + 96 + 26 @@ -88,8 +103,8 @@ 10 25 - 58 - 16 + 101 + 26 @@ -102,7 +117,7 @@ 132 25 156 - 21 + 26 @@ -110,9 +125,9 @@ 132 - 50 + 55 156 - 21 + 26 @@ -120,9 +135,9 @@ 10 - 50 - 116 - 16 + 55 + 101 + 26 @@ -130,26 +145,13 @@ - - - - 0 - 85 - 426 - 16 - - - - Qt::Horizontal - - 0 0 - 426 + 446 22 diff --git a/sfeprapy/guilayout/ui2py.py b/sfeprapy/gui/layout/ui2py.py similarity index 85% rename from sfeprapy/guilayout/ui2py.py rename to sfeprapy/gui/layout/ui2py.py index 466a344..6dff975 100644 --- a/sfeprapy/guilayout/ui2py.py +++ b/sfeprapy/gui/layout/ui2py.py @@ -8,7 +8,7 @@ def ui2py(): 'main.ui', ] - cwd = os.path.dirname(os.path.realpath(__file__)) + 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: diff --git a/sfeprapy/gui/logic/__init__.py b/sfeprapy/gui/logic/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/sfeprapy/guilogic/main.py b/sfeprapy/gui/logic/main.py similarity index 87% rename from sfeprapy/guilogic/main.py rename to sfeprapy/gui/logic/main.py index c2842c7..b890648 100644 --- a/sfeprapy/guilogic/main.py +++ b/sfeprapy/gui/logic/main.py @@ -1,8 +1,8 @@ import sys -from PySide2 import QtWidgets +from PySide2 import QtWidgets, QtGui, QtCore -from sfeprapy.guilayout.main import Ui_MainWindow +from sfeprapy.gui.layout.main import Ui_MainWindow from sfeprapy.mcs0 import EXAMPLE_INPUT_CSV from sfeprapy.mcs0.__main__ import main as mcs0 @@ -15,8 +15,14 @@ def __init__(self): 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) @@ -92,6 +98,6 @@ def main(app_window=None): if __name__ == "__main__": """update ui files and run the application. this is used only for testing""" - from sfeprapy.guilayout.ui2py import ui2py + from sfeprapy.gui.layout.ui2py import ui2py ui2py() main() diff --git a/sfeprapy/mcs0/__init__.py b/sfeprapy/mcs0/__init__.py index 91f25c3..0f875f9 100644 --- a/sfeprapy/mcs0/__init__.py +++ b/sfeprapy/mcs0/__init__.py @@ -5,8 +5,7 @@ def __example_config_dict(): - y = dict(n_threads=2) - return y + return dict(n_threads=2, cwd='') def __example_input_dict(): diff --git a/sfeprapy/mcs0/mcs0_calc.py b/sfeprapy/mcs0/mcs0_calc.py index 95995c2..457e226 100644 --- a/sfeprapy/mcs0/mcs0_calc.py +++ b/sfeprapy/mcs0/mcs0_calc.py @@ -7,6 +7,7 @@ import pandas as pd 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 @@ -16,6 +17,7 @@ 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): @@ -436,7 +438,7 @@ def f_(x, terminate_check_wait_time): ): # Terminate if maximum solving iteration is reached if solver_temperature_goal > y3: - solver_time_equivalence_solved = -np.inf + solver_time_equivalence_solved = 0. elif solver_temperature_goal < y3: solver_time_equivalence_solved = np.inf else: @@ -556,13 +558,24 @@ def mcs_out_post(df: pd.DataFrame) -> pd.DataFrame: x = df_res[k].values x1, x2, x3 = np.min(x), np.mean(x), np.max(x) dict_[k] = f"{x1:<9.3f} {x2:<9.3f} {x3:<9.3f}" - except (KeyError, ValueError): + except Exception: pass list_ = [f"{k:<24.24}: {v}" for k, v in dict_.items()] print("\n".join(list_), "\n") + try: + x = np.array(df_res['solver_time_equivalence_solved'].values / 60, dtype=float) + x[x == -np.inf] = 0 + x[x == np.inf] = np.amax(x[x != np.inf]) + 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()) + except Exception as e: + print(f'Failed to plot time equivalence, {e}') + return df @@ -632,6 +645,16 @@ def teq_main( + 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) @@ -643,7 +666,7 @@ def teq_main( 345.0 * np.log10((fire_time / 60.0) * 8.0 + 1.0) + 20.0 ) + 273.15 # in [K] - # Inject results, i.e. these values will be in the output + # Inject results, i.e. these input values will also be in the output file res = dict( case_name=case_name, n_simulations=n_simulations, @@ -669,8 +692,8 @@ def teq_main( # protection_k=protection_k, # protection_protected_perimeter=protection_protected_perimeter, # protection_rho=protection_rho, - # room_breadth=room_breadth, - # room_depth=room_depth, + 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, @@ -678,9 +701,9 @@ def teq_main( # solver_thickness_lbound=solver_thickness_lbound, # solver_thickness_ubound=solver_thickness_ubound, # solver_tol=solver_tol, - # window_height=window_height, + window_height=window_height, window_open_fraction=window_open_fraction, - # window_width=window_width, + window_width=window_width, phi_teq=phi_teq, # timber_charring_rate=timber_charring_rate, # timber_hc=timber_hc, @@ -776,23 +799,15 @@ def teq_main( phi_teq, ) + # additional fuel contribution from timber + if timber_exposed_area <= 0 or timber_exposed_area is None: # no timber exposed break - elif not res_solve_time_equivalence[ - "solver_convergence_status" - ]: # no time equivalence solution + elif not res_solve_time_equivalence["solver_convergence_status"]: # no time equivalence solution break - elif ( - timber_solver_iter >= timber_solver_ilim - ): # over the solver iteration limit + elif (timber_solver_iter >= timber_solver_ilim): # over the solver iteration limit 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- res_solve_time_equivalence["solver_time_equivalence_solved"])<= timber_solver_tol): # convergence sought successfully break else: timber_exposed_duration = res_solve_time_equivalence[ @@ -817,6 +832,20 @@ def teq_main( return res +class MCS0(MCS): + 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 mcs_post(self, *args, **kwargs): + return mcs_out_post(*args, **kwargs) + + def _test_teq_phi(): warnings.filterwarnings("ignore") @@ -924,3 +953,91 @@ def _test_standard_case(): print(teq_at_80_percentile) target, target_tol = 60, 2 assert target - target_tol < teq_at_80_percentile < target + target_tol + + +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 + + # 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.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 _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 + + # 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 + 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 + assert target - target_tol < teq_at_80_percentile < target + target_tol diff --git a/test/test_func_travelling.py b/test/test_func_travelling.py index fae1259..72bab91 100644 --- a/test/test_func_travelling.py +++ b/test/test_func_travelling.py @@ -8,7 +8,10 @@ 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 49096fc..815797b 100644 --- a/test/test_mcs0.py +++ b/test/test_mcs0.py @@ -2,6 +2,8 @@ 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