Merge pull request #10 from GEMScienceTools/slf-generator

Add SLF generation functionality

Author: mouayed-nafeh

.github/workflows/linux_test.yml

@@ -40,6 +40,10 @@ jobs:
       - name: Install vmtk Package
         run: |
           source ~/openquake/bin/activate
+          PY_VER=`echo py${{ matrix.python-version }} | tr -d .`
+          echo $PY_VER
+          export PIP_DEFAULT_TIMEOUT=100
+          pip install -r requirements-$PY_VER-linux.txt
           pip install -e .
       - name: Run tests
         run: |

.github/workflows/windows_test.yml

@@ -22,7 +22,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [windows-2022]
-        python-version: ["3.10", "3.11", "3.12"]
+        python-version: ["3.11", "3.12"]
     steps:
       - uses: actions/checkout@v4
       - name: Set up Python  ${{ matrix.python-version }}
@@ -37,6 +37,10 @@ jobs:
         run: |
           Set-ExecutionPolicy -ExecutionPolicy RemoteSigned -Scope CurrentUser
           C:\Users\runneradmin\openquake\Scripts\activate.ps1
+          $PY_VER="py${{ matrix.python-version }}"
+          $py = $PY_VER.replace(".","")
+          set PIP_DEFAULT_TIMEOUT=100
+          python -m pip install -r requirements-$py-win64.txt
           pip install -e .
       - name: Run tests
         run: |

.virtual_documents/demos/example_2.ipynb

@@ -0,0 +1,305 @@
+
+
+
+
+
+
+import os
+import sys
+import shutil
+import numpy as np
+import pandas as pd
+import matplotlib.pyplot as plt
+
+# Import the classes necessary for structural analysis
+from openquake.vmtk.units         import units              # oq-vtmk units class
+from openquake.vmtk.calibration   import calibrate_model    # oq-vmtk sdof-to-mdof calibration class
+from openquake.vmtk.modeller      import modeller           # oq-vmtk numerical modelling class
+from openquake.vmtk.postprocessor import postprocessor      # oq-vtmk postprocessing class
+from openquake.vmtk.plotter       import plotter            # oq-vmtk plotting class
+from openquake.vmtk.utilities     import sorted_alphanumeric, import_from_pkl, export_to_pkl # oq-vmtk utility class
+
+
+
+
+
+# Define the directory of the ground-motion records
+gm_directory  = './in/records'            
+
+# Define the main output directory
+nrha_directory = './out/nltha'  
+os.makedirs(nrha_directory, exist_ok=True)
+
+# Define directory for temporary analysis outputs: it is used to store temporary .txt files used as accelerations recorders
+temp_nrha_directory = os.path.join(nrha_directory,'temp')
+os.makedirs(temp_nrha_directory, exist_ok=True)
+
+
+
+
+
+# Import the intensity measure dictionary (output from example 1)
+ims           = import_from_pkl(os.path.join(gm_directory, 'imls_esrm20.pkl'))  
+
+
+
+
+
+
+
+
+# Number of storeys
+number_storeys = 2
+
+# Relative floor heights list
+floor_heights = [2.80, 2.80]
+
+# First-mode based participation factor
+gamma = 1.33
+
+# SDOF capacity (First row are Spectral Displacement [m] values - Second row are Spectral Acceleration [g] values)
+sdof_capacity = np.array([[0.00060789, 0.00486316, 0.02420000, 0.04353684], 
+                          [0.10315200, 0.20630401, 0.12378241, 0.12502023]]).T
+# Frame flag
+isFrame = False
+
+# Soft-storey mechanism flag
+isSOS = False
+
+# Degradation flag 
+mdof_degradation = True
+
+# Inherent damping 
+mdof_damping = 0.05
+
+
+
+
+
+# Intensity measures to use for postprocessing cloud analyses
+IMTs      = ['PGA', 'SA(0.3s)', 'SA(0.6s)', 'SA(1.0s)','AvgSA']
+
+# Damage thresholds (maximum peak storey drift values in rad)
+damage_thresholds    =  [0.00150, 0.00545, 0.00952, 0.0135] 
+
+# The lower limit to be applied for censoring edp values (below 0.1 the minimum threshold for slight damage is considered a negligible case)
+lower_limit = 0.1*damage_thresholds[0]
+
+# The upper limit to be applied for consoring edp values (above 1.5 the maximum threshold is considered a collapse case) 
+censored_limit = 1.5*damage_thresholds[-1]   
+
+# Define consequence model to relate structural damage to a decision variable (i.e., expected loss ratio) 
+consequence_model = [0.05, 0.20, 0.60, 1.00] # damage-to-loss ratios
+
+
+
+
+
+
+
+
+
+
+
+# Calibrate the model using the Lu et al. (2020) method
+floor_masses, storey_disps, storey_forces, mdof_phi = calibrate_model(number_storeys, gamma, sdof_capacity, isFrame, isSOS)
+
+print('The mass of each floor (in tonnes):', floor_masses)
+print('The first-mode shape used for calibration:', mdof_phi)
+
+# Plot the capacities to visualise the outcome of the calibration
+for i in range(storey_disps.shape[0]):
+   plt.plot(np.concatenate(([0.0], storey_disps[i,:])), np.concatenate(([0.0], storey_forces[i,:]*9.81)), label = f'Storey #{i+1}')
+plt.plot(np.concatenate(([0.0], sdof_capacity[:,0])), np.concatenate(([0.0], sdof_capacity[:,1]*9.81)), label = 'SDOF Capacity')
+plt.xlabel('Storey Deformation [m]', fontsize= 16)
+plt.ylabel('Storey Shear [kN]', fontsize = 16)
+plt.legend(loc = 'lower right')
+plt.grid(visible=True, which='major')
+plt.grid(visible=True, which='minor')
+plt.xlim([0.00, 0.03])
+plt.show()
+
+
+
+
+
+# Initialise MDOF storage lists
+conv_index_list = []               # List for convergence indices
+peak_disp_list  = []               # List for peak floor displacement (returns all peak values along the building height)
+peak_drift_list = []               # List for peak storey drift (returns all peak values along the building height)
+peak_accel_list = []               # List for peak floor acceleration (returns all peak values along the building height)
+max_peak_drift_list = []           # List for maximum peak storey drift (returns the maximum value) 
+max_peak_drift_dir_list = []       # List for maximum peak storey drift directions
+max_peak_drift_loc_list = []       # List for maximum peak storey drift locations
+max_peak_accel_list = []           # List for maximum peak floor acceleration (returns the maximum value)
+max_peak_accel_dir_list = []       # List for maximum peak floor acceleration directions 
+max_peak_accel_loc_list = []       # List for maximum peak floor acceleration locations
+
+# Loop over ground-motion records, compile MDOF model and run NLTHA
+gmrs = sorted_alphanumeric(os.listdir(os.path.join(gm_directory,'acc')))                         # Sort the ground-motion records alphanumerically
+dts  = sorted_alphanumeric(os.listdir(os.path.join(gm_directory,'dts')))                         # Sort the ground-motion time-step files alphanumerically
+
+# Run the analysis
+for i in range(len(gmrs)):
+    ### Print post-processing iteration
+    print('================================================================')
+    print('============== Analysing: {:d} out of {:d} =================='.format(i+1, len(gmrs)))
+    print('================================================================')
+
+    ### Compile the MDOF model    
+    model = modeller(number_storeys,
+                     floor_heights,
+                     floor_masses,
+                     storey_disps,
+                     storey_forces*units.g,
+                     mdof_degradation)                                                                # Initialise the class (Build the model)
+    
+    model.compile_model()                                                                             # Compile the MDOF model
+    
+    if i==0:
+        model.plot_model()                                                                            # Visualise the model (only on first iteration)        
+    model.do_gravity_analysis()                                                                       # Do gravity analysis
+
+    if number_storeys == 1:
+        num_modes = 1
+    else:
+        num_modes = 3
+    T, phi = model.do_modal_analysis(num_modes = num_modes)                                           # Do modal analysis and get period of vibration (Essential step for running NLTHA)
+
+    ### Define ground motion objects
+    fnames = [os.path.join(gm_directory,'acc',f'{gmrs[i]}')]                                          # Ground-motion record names
+    fdts   =  os.path.join(gm_directory,'dts',f'{dts[i]}')                                            # Ground-motion time-step names 
+    dt_gm = pd.read_csv(fdts, header=None)[pd.read_csv(fdts,header=None).columns[0]].loc[1]-\
+            pd.read_csv(fdts, header=None)[pd.read_csv(fdts,header=None).columns[0]].loc[0]           # Ground-motion time-step
+    t_max = pd.read_csv(fdts)[pd.read_csv(fdts).columns[0]].iloc[-1]                                  # Ground-motion duration
+    
+    ### Define analysis params and do NLTHA
+    dt_ansys = dt_gm                                                            # Set the analysis time-step
+    sf = units.g                                                                # Set the scaling factor (if records are in g, a scaling factor of 9.81 m/s2 must be used to be consistent with opensees) 
+    control_nodes, conv_index, peak_drift, peak_accel, max_peak_drift, max_peak_drift_dir, max_peak_drift_loc, max_peak_accel, max_peak_accel_dir, max_peak_accel_loc, peak_disp = model.do_nrha_analysis(fnames, 
+                                                                                                                                                                                                          dt_gm, 
+                                                                                                                                                                                                          sf, 
+                                                                                                                                                                                                          t_max, 
+                                                                                                                                                                                                          dt_ansys,
+                                                                                                                                                                                                          temp_nrha_directory,
+                                                                                                                                                                                                          pflag=False,
+                                                                                                                                                                                                          xi = mdof_damping)
+
+    ### Store the analysis
+    conv_index_list.append(conv_index)
+    peak_drift_list.append(peak_drift)
+    peak_accel_list.append(peak_accel)
+    peak_disp_list.append(peak_disp)
+    max_peak_drift_list.append(max_peak_drift)
+    max_peak_drift_dir_list.append(max_peak_drift_dir)
+    max_peak_drift_loc_list.append(max_peak_drift_loc)
+    max_peak_accel_list.append(max_peak_accel)
+    max_peak_accel_dir_list.append(max_peak_accel_dir)
+    max_peak_accel_loc_list.append(max_peak_accel_loc)
+
+# Remove the temporary directory
+shutil.rmtree(f'{temp_nrha_directory}')
+
+# Store the analysis results in a dictionary
+ansys_dict = {}
+labels = ['T','control_nodes', 'conv_index_list',
+          'peak_drift_list','peak_accel_list',
+          'max_peak_drift_list', 'max_peak_drift_dir_list', 
+          'max_peak_drift_loc_list','max_peak_accel_list',
+          'max_peak_accel_dir_list','max_peak_accel_loc_list',
+          'peak_disp_list']
+
+for i, label in enumerate(labels):
+    ansys_dict[label] = vars()[f'{label}']
+# Export the analysis output variable to a pickle file using the "export_to_pkl" function from "utilities"
+export_to_pkl(os.path.join(nrha_directory,'ansys_out.pkl'), ansys_dict) 
+
+print('ANALYSIS COMPLETED!')
+
+
+
+
+
+# Initialise the postprocessor class
+pp = postprocessor()
+
+# Initialise the plotter class
+pl = plotter()
+
+# Loop over the intensity measure types and perform cloud regression to fit the probabilistic seismic demand-capacity model
+for _, current_imt in enumerate(IMTs):
+    
+    # Import the current intensity measure type
+    imls = ims[f'{current_imt}']                   
+
+    # Import the engineering demand parameters (i.e., mpsd) from the analysis dictionary (processed from example 2)
+    edps = ansys_dict['max_peak_drift_list']  
+    
+    # Process cloud analysis results using the "do_cloud_analysis" function called from "postprocessor" 
+    # The output will be automatically stored in a dictionary
+    cloud_dict = pp.do_cloud_analysis(imls,
+                                      edps,
+                                      damage_thresholds,
+                                      lower_limit,
+                                      censored_limit) 
+        
+    ## Visualise the cloud analysis results
+    pl.plot_cloud_analysis(cloud_dict, 
+                          output_directory = None, 
+                          plot_label = f'cloud_analysis_{current_imt}',
+                          xlabel = f'{current_imt} [g]', 
+                          ylabel = r'Maximum Peak Storey Drift, $\theta_{max}$ [%]') # The y-axis values of drift are converted to % automatically by the plotter
+
+    ## Visualise the fragility functions
+    pl.plot_fragility_analysis(cloud_dict,
+                               output_directory = None,
+                               plot_label = f'fragility_{current_imt}',
+                               xlabel = f'{current_imt}')
+
+    ## Visualise the seismic demands
+    pl.plot_demand_profiles(ansys_dict['peak_drift_list'], 
+                            ansys_dict['peak_accel_list'], 
+                            ansys_dict['control_nodes'], 
+                            output_directory = None,
+                            plot_label="seismic_demand_profiles") # The y-axis values of drift and acceleration are converted to % and g automatically by the plotter
+        
+    ## Visualise the entire set of results using subplots
+    pl.plot_ansys_results(cloud_dict,
+                          ansys_dict['peak_drift_list'],
+                          ansys_dict['peak_accel_list'],
+                          ansys_dict['control_nodes'],
+                          output_directory = None,
+                          plot_label = f'analysis_output_{current_imt}',
+                          cloud_xlabel = f'{current_imt}',
+                          cloud_ylabel = r'Maximum Peak Storey Drift, $\theta_{max}$ [%]')
+
+
+
+
+
+# In this example, since the latest iteration of the previous cell uses 'AvgSA' as the intensity measure,
+# then all variables stored inside the "cloud_dict" dictionary correspond to that same IM. Hence, 
+# the vulnerability function derived here will represent the continuous relationship of the expected 
+# structural loss ratio conditioned on increasing levels of ground-shaking expressed in terms of the 
+# average spectral acceleration (in g)
+
+structural_vulnerability = pp.get_vulnerability_function(cloud_dict['poes'],
+                                                         consequence_model,
+                                                         uncertainty=True)
+
+
+# Plot the structural vulnerability function
+pl.plot_vulnerability_analysis(structural_vulnerability['IMLs'],
+                               structural_vulnerability['Loss'],
+                               structural_vulnerability['COV'],
+                               'SA(1.0s)',
+                               'Structural Loss Ratio',
+                               output_directory = None,
+                               plot_label = 'Structural Vulnerability')
+
+
+# The output is a DataFrame with three keys: IMLs (i.e., intensity measure levels), Loss and COV
+print(structural_vulnerability)
+
+
+