Same (or very similar) time domain postprocessing settings in all too…

…ls -> downsampling option added, uq_results_analysis -> predefined colors option added + pie charts sobol indices

use_cases/stability_analysis/config/alaska/config_alaska_case_study_1.txt

@@ -78,7 +78,7 @@
     [[postprocessor]]
         # signal selection
         var_name = torsion_b1, torsion_b2, torsion_b3
-        interp_radpos = 70,
+        interp_radpos = 50, 55, 60, 65, 70, 75
 
         # windowing method:
         # 0 -> windowing based on user defined crop_window
@@ -97,12 +97,13 @@
         # options are: DMD, log_decr, exp_fit, linear_fit
         damping_method = DMD
         # number of svd modes to be retained in DMD, 0 -> automatic determination by pyDMD
-        svd_rank = 0
+        svd_rank = 20
 
         # enable or disable MBC transformation
-        MBC_transf = False
+        MBC_transf = True
 
         plotting = True
 
-
+        downsampling_flag = True
+        downsampling_freq = 8  # Hz
 

use_cases/stability_analysis/config/alaska/config_alaska_case_study_2.txt

@@ -145,12 +145,12 @@
         # options are: DMD, log_decr, exp_fit, linear_fit
         damping_method = DMD
         # number of svd modes to be retained in DMD, 0 -> automatic determination by pyDMD
-        svd_rank = 0
+        svd_rank = 20
 
         # enable or disable MBC transformation
-        MBC_transf = false
-        MBC_only_collective = false
-        MBC_only_cyclic = false
+        MBC_transf = True
 
         plotting = True
 
+        downsampling_flag = True
+        downsampling_freq = 8  # Hz

use_cases/stability_analysis/config/bladed/config_bladed_case_study_1.txt

@@ -77,8 +77,12 @@
 
         # damping determination method
         damping_method = DMD
-        svd_rank = 0
+        svd_rank = 20
 
         plotting = True
 
-        MBC_transf = False
+        # enable or disable MBC transformation
+        MBC_transf = True
+
+        downsampling_flag = False
+        downsampling_freq = 8  # Hz

use_cases/stability_analysis/config/bladed/config_bladed_case_study_2.txt

@@ -126,8 +126,11 @@
 
         # damping determination method
         damping_method = DMD
-        svd_rank = 0
+        svd_rank = 20
 
         plotting = True
 
-        MBC_transf = False
+        MBC_transf = True
+
+        downsampling_flag = False
+        downsampling_freq = 8  # Hz

use_cases/stability_analysis/config/hawc2/config_hawc2_case_study_1.txt

@@ -3,7 +3,7 @@
     [[uncertainty]]
         base_directory = ./reference_data/hawc2
         run_directory = case_study_1
-        n_CPU = 1
+        n_CPU = 8
 
     [[parameters]]
         [[[flap_stiffness]]]
@@ -94,11 +94,13 @@
 
         # damping determination method
         damping_method = DMD
-        svd_rank = 0
+        svd_rank = 20
 
-        plotting = False
-
-        MBC_transf = False
+        plotting = True
 
+        MBC_transf = True
+
+        downsampling_flag = True
+        downsampling_freq = 8  # Hz
 
 

use_cases/stability_analysis/config/hawc2/config_hawc2_case_study_2.txt

@@ -146,8 +146,11 @@
 
         # damping determination method
         damping_method = DMD
-        svd_rank = 0
+        svd_rank = 20
 
-        plotting = False
+        plotting = True
 
-        MBC_transf = False
+        MBC_transf = True
+
+        downsampling_flag = True
+        downsampling_freq = 8  # Hz

use_cases/stability_analysis/config/simpack/config_simpack_case_study_1.txt

@@ -87,4 +87,7 @@
 
         plotting = True
 
-        MBC_transf = False
+        MBC_transf = True
+
+        downsampling_flag = True
+        downsampling_freq = 8  # Hz

use_cases/stability_analysis/config/simpack/config_simpack_case_study_2.txt

@@ -138,4 +138,7 @@
 
         plotting = True
 
-        MBC_transf = False
+        MBC_transf = True
+
+        downsampling_flag = True
+        downsampling_freq = 8  # Hz

use_cases/stability_analysis/postprocessor.py

@@ -8,7 +8,7 @@
 import matplotlib.pyplot as plt
 import numpy as np
 import os 
-from scipy.signal import find_peaks  # new feature (scipy 1.1.0)
+from scipy.signal import find_peaks, resample  # find peaks = new feature (scipy 1.1.0)
 from scipy.optimize import curve_fit
 from scipy.interpolate import interp1d
 import logging
@@ -69,7 +69,8 @@ def __init__(self, resultdict, time=np.array([]), sampling_time=float(), signal=
         self.signal = signal
         self.result_dir = result_dir
 
-    def read_results(self, var_names, interp_radpos=None, nr_radpos=None):
+    def read_results(self, var_names, interp_radpos=None, nr_radpos=None,
+                     downsampling_flag=False, downsampling_frequency=None):
         """
         Read result dictionary
 
@@ -82,6 +83,10 @@ def read_results(self, var_names, interp_radpos=None, nr_radpos=None):
             Radial positions for interpolation of a 2D data set
         nr_radpos : int, optional
             Number of equidistantly spaced radial positions to be used for each of the variables
+        downsampling_flag : bool, optional
+            Flag to indicate if the time series have to be downsampled
+        downsampling_frequency : float, optional
+            Sample rate to downsample to (Hz!)
 
         Raises
         ------
@@ -96,20 +101,14 @@ def read_results(self, var_names, interp_radpos=None, nr_radpos=None):
             - (radpos: only if var_name is a 2D signal)
         """
         logger = logging.getLogger('quexus.uq.postprocessor.read_results')
-        # read time and determine sampling time
-        self.time = np.squeeze(self.resultdict['time'])
-        # determine sampling time from time signal
-        # sometimes the time stamps are not accurate enough (in Bladed f.ex.: 0.0, 0.009990785, 0.02, 0.03, ...), this
-        # might lead to problems later on
-        dt1 = self.time[1] - self.time[0]
-        dtfull = (self.time[-1] - self.time[0]) / (np.size(self.time) - 1)
-        if dt1 != dtfull:
-            logger.debug('sampling time determined as time[1] - time[0] is not the same as '
-                         'time[end] - time[0] / total time')
-            logger.debug(str(dtfull) + 's used as sampling time')
-        self.sampling_time = dtfull
-        logger.debug('Sampling time = {}'.format(self.sampling_time))
-
+
+        # DOWNSAMPLING
+        if downsampling_flag is True:
+            logger.debug('Time series will be downsampled to sampling frequency {} Hz'.format(downsampling_frequency))
+            self.downsample_resultdict(downsampling_frequency)
+
+        self.time, self.sampling_time = self.get_time_details_resultdict()
+
         radpos_list = []
         signal_list = []
         for var_name in var_names:
@@ -152,7 +151,71 @@ def read_results(self, var_names, interp_radpos=None, nr_radpos=None):
             self.radpos = np.array(radpos_list)
         else:
             logger.debug('only 1D arrays used in input, no radpos available')
-            self.radpos = np.array([])      
+            self.radpos = np.array([])
+
+    def get_time_details_resultdict(self):
+        """
+        determine characteristics of the time series in the resultdict
+
+        Returns
+        -------
+        time : 1D array
+            interpolated signal at the interp_radpos
+        sampling_time : float
+            sampling time of the time signal in the resultdict
+        """
+        logger = logging.getLogger('quexus.uq.postprocessor.get_time_details_resultdict')
+
+        time = np.squeeze(self.resultdict['time'])
+
+        # determine sampling time from time signal
+        # sometimes the time stamps are not accurate enough (in Bladed f.ex.: 0.0, 0.009990785, 0.02, 0.03, ...), this
+        # might lead to problems later on
+        dt1 = time[1] - time[0]
+        dtfull = (time[-1] - time[0]) / (np.size(time) - 1)
+        if dt1 != dtfull:
+            logger.debug('sampling time determined as time[1] - time[0] is not the same as '
+                         'time[end] - time[0] / total time')
+            logger.debug(str(dtfull) + 's used as sampling time')
+        sampling_time = dtfull
+        logger.debug('Sampling time = {}'.format(sampling_time))
+
+        return time, sampling_time
+
+    def downsample_resultdict(self, downsampling_frequency):
+        """
+        Downsample the time series which are provided in the resultdict
+
+        Parameters
+        ----------
+        downsampling_frequency : float
+            sampling frequency (in Hz) to downsample to
+
+        Returns
+        -------
+        """
+        orig_time, orig_sampling_time = self.get_time_details_resultdict()
+
+        nr_new_samples = round(np.size(orig_time) * downsampling_frequency / (1/orig_sampling_time))
+
+        for key in self.resultdict:
+            if key == 'time':
+                continue
+
+            if np.size(orig_time) not in self.resultdict[key].shape:
+                continue
+            else:
+                dim = find_corresponding_axis(self.resultdict[key], np.size(orig_time), 3)
+                self.resultdict[key] = resample(self.resultdict[key], nr_new_samples, axis=dim)
+
+        self.resultdict['time'] = np.linspace(np.squeeze(self.resultdict['time'])[0],
+                                              np.squeeze(self.resultdict['time'])[-1],
+                                              nr_new_samples, endpoint=False)
+
+        # self.resultdict['time'] = self.resultdict['time'][::10]
+        # self.resultdict['torsion_b1'] = self.resultdict['torsion_b1'][:, ::10]
+        # self.resultdict['torsion_b2'] = self.resultdict['torsion_b2'][:, ::10]
+        # self.resultdict['torsion_b3'] = self.resultdict['torsion_b3'][:, ::10]
 
     @staticmethod
     def interpolate_signals(signal, radpos, interp_radpos=None, nr_radpos=None):

use_cases/stability_analysis/preprocessor.py

@@ -252,6 +252,8 @@ def modify_structural_model(self, uncertain_parameter, spline_object, method):
 
             elif uncertain_parameter == 'shear_center_x':
                 if method == 'multiplicative':
+                    # This spline application is chosen very unlucky... Now the input changes sign, so when the
+                    # spline is >1 the delta actually moves in negative direction...
                     delta_sc_x = (1-self.apply_spline(spline_object, css.arclength)[0]) * \
                                  self.get_local_chord(css.arclength)
                 elif method == 'additive':
@@ -261,6 +263,8 @@ def modify_structural_model(self, uncertain_parameter, spline_object, method):
 
             elif uncertain_parameter == 'shear_center_y':
                 if method == 'multiplicative':
+                    # This spline application is chosen very unlucky... Now the input changes sign, so when the
+                    # spline is >1 the delta actually moves in negative direction...
                     delta_sc_y = (1-self.apply_spline(spline_object, css.arclength)[0]) * \
                                  self.get_local_chord(css.arclength)
                 elif method == 'additive':
@@ -270,6 +274,8 @@ def modify_structural_model(self, uncertain_parameter, spline_object, method):
 
             elif uncertain_parameter == 'cog_x':
                 if method == 'multiplicative':
+                    # This spline application is chosen very unlucky... Now the input changes sign, so when the
+                    # spline is >1 the delta actually moves in negative direction...
                     delta_cog_x = (1-self.apply_spline(spline_object, css.arclength)[0]) * \
                                   self.get_local_chord(css.arclength)
                 elif method == 'additive':
@@ -279,6 +285,8 @@ def modify_structural_model(self, uncertain_parameter, spline_object, method):
 
             elif uncertain_parameter == 'cog_y':
                 if method == 'multiplicative':
+                    # This spline application is chosen very unlucky... Now the input changes sign, so when the
+                    # spline is >1 the delta actually moves in negative direction...
                     delta_cog_y = (1-self.apply_spline(spline_object, css.arclength)[0]) * \
                                   self.get_local_chord(css.arclength)
                 elif method == 'additive':

use_cases/stability_analysis/result_comparison/config_case_study_1.txt

@@ -1,10 +1,12 @@
 
 base_directory = .\reference_data
-uncertainpy_results = bladed-lin\case_study_1\uq_results\StabilityAnalysisModel.h5, bladed\case_study_1\uq_results\StabilityAnalysisModel.h5, alaska\case_study_1\uq_results\StabilityAnalysisModel.h5, hawc2\case_study_1\uq_results\StabilityAnalysisModel.h5, hawcstab2\case_study_1\uq_results\StabilityAnalysisModel.h5, simpack\case_study_1\uq_results\StabilityAnalysisModel.h5
+uncertainpy_results = hawc2\case_study_1\uq_results\StabilityAnalysisModel.h5, hawcstab2\case_study_1\uq_results\StabilityAnalysisModel.h5, bladed\case_study_1\uq_results\StabilityAnalysisModel.h5, bladed-lin\case_study_1\uq_results\StabilityAnalysisModel.h5, alaska\case_study_1\uq_results\StabilityAnalysisModel.h5, simpack\case_study_1\uq_results\StabilityAnalysisModel.h5
 
 model_name = StabilityAnalysisModel
 
-label_names = Bladed (lin.), Bladed, alaska/Wind, HAWC2, HAWCStab2 (lin.), Simpack
+label_names = HAWC2, HAWCStab2 (lin.), Bladed, Bladed (lin.), alaska/Wind, Simpack
+
+predefined_colors = C2, C5, C1, C6, C0, C4
 
 output_dir = result_comparison\comparison_case_study_1
 

use_cases/stability_analysis/result_comparison/config_case_study_2.txt

@@ -1,12 +1,14 @@
 
 base_directory = .\reference_data
-uncertainpy_results = bladed-lin\case_study_2\uq_results\StabilityAnalysisModel.h5, bladed\case_study_2\uq_results\StabilityAnalysisModel.h5, alaska\case_study_2\uq_results\StabilityAnalysisModel.h5, hawc2\case_study_2\uq_results\StabilityAnalysisModel.h5, hawcstab2\case_study_2\uq_results\StabilityAnalysisModel.h5, simpack\case_study_2\uq_results\StabilityAnalysisModel.h5
+uncertainpy_results = hawc2\case_study_2\uq_results\StabilityAnalysisModel.h5, hawcstab2\case_study_2\uq_results\StabilityAnalysisModel.h5, bladed\case_study_2\uq_results\StabilityAnalysisModel.h5, bladed-lin\case_study_2\uq_results\StabilityAnalysisModel.h5, alaska\case_study_2\uq_results\StabilityAnalysisModel.h5, simpack\case_study_2\uq_results\StabilityAnalysisModel.h5
 
 model_name = StabilityAnalysisModel
 
-label_names = Bladed (lin.), Bladed, alaska/Wind, HAWC2, HAWCStab2 (lin.), Simpack
+label_names = HAWC2, HAWCStab2 (lin.), Bladed, Bladed (lin.), alaska/Wind, Simpack
 
-output_dir = result_comparison\comparison_case_study_1
+predefined_colors = C2, C5, C1, C6, C0, C4
+
+output_dir = result_comparison\comparison_case_study_2
 
 bokeh_plot = True
 
@@ -16,4 +18,4 @@ bokeh_plot = True
 edge_stiffness_fixed_distr = Edge Stiffness
 cog_x_fixed_distr = Chordwise COG
 shear_center_x_fixed_distr = Chordwise SC
-principle_axis_orientation_cp01_y = PA Orientation
+principle_axis_orientation_cp01_y = PA Orientation

use_cases/stability_analysis/run_analysis.py

@@ -5,7 +5,7 @@
 import os
 import logging
 
-#from preprocessor import PreProcessor
+from preprocessor import PreProcessor
 from postprocessor import PostProcessor
 from libDMD import DMDanalysis
 from libMBC import MBCTransformation
@@ -111,8 +111,12 @@ def _postprocessor(self, resultdict, iteration_run_directory):
         # args:
         # - variable name (key of resultdict),
         # - interp_radpos = optional radpos positions for interpolation of 2D array
+        # - downsampling = downsample time series
+        # - downsampling_frequency = frequency to downsample to
         postprocess.read_results(self.postprocessor_config['var_name'],
-                                 interp_radpos=interp_radpos)
+                                 interp_radpos=interp_radpos,
+                                 downsampling_flag=self.postprocessor_config['downsampling_flag'],
+                                 downsampling_frequency=self.postprocessor_config['downsampling_freq'])
 
         # do some pre-processing on the signals: de-mean, crop_window
         postprocess.prepare_signals(window_method=window_method,

wtuq_framework/helperfunctions.py

@@ -227,7 +227,7 @@ def find_corresponding_axis(arr, size, max_dim):
     """
 
     if len(arr.shape)>max_dim:
-        raise ValueError('Only' + str(max_dim) + '-dimensional arrays are allowed')
+        raise ValueError('Only ' + str(max_dim) + '-dimensional arrays are allowed')
 
     if arr.shape.count(size) == 0:
         print('WARNING: Dimensions of arrays to be correlated do not match')