diff --git a/mssrc.py b/mssrc.py index 8f8f44b93e87115399ed3b9306b317e1709a29bb..b02565f8c7248c7ce0e6d4d63d03bf633dd00dca 100644 --- a/mssrc.py +++ b/mssrc.py @@ -9,8 +9,11 @@ import numpy as np from numpy import linalg import pandas as pd import time +import math import matplotlib as plt import sys +import decimal +from decimal import * try: from dataio.H5Writer import H5Writer from dataio.H5Reader import H5Reader @@ -33,7 +36,34 @@ except ImportError: """ # Define global variables -COEFF_0 = 0.3 # Eq. (43) in Shen et al. (2005): low limit constant +COEFF_0 = Decimal(0.3) # Eq. (43) in Shen et al. (2005): low limit constant + +def inverse_den(x): + """ + Calculate inverse of a number. + """ + if abs(x) < Decimal(1.e-24): + return Decimal(0.0) + else: + return Decimal(1.0) / x + +def calc_det(mat): + """ + Takes a matrix as and argument and returns the determinant + + In this function, the configuration JSON-File and the the signal data are parsed. The from the signal time series, the local velocities, the void fraction, the bubble diameters and the interfacial area density are recovered by a reconstruction algorithm. + """ + a = mat[0][0] + b = mat[0][1] + c = mat[0][2] + d = mat[1][0] + e = mat[1][1] + f = mat[1][2] + g = mat[2][0] + h = mat[2][1] + i = mat[2][2] + det = a*e*i + b*f*g + c*d*h - c*e*g - b*d*i - a*f*h + return det def main(): """ @@ -56,7 +86,7 @@ def main(): args = parser.parse_args() # Eq. (43) in Shen et al. (2005): estimated gas velocity - V_GAS = float(args.velocity) + V_GAS = Decimal(args.velocity) # Create Posix path for OS indepency path = pathlib.Path(args.path) @@ -86,13 +116,16 @@ def main(): if s_id != 0: aux_sensor_ids.append(s_id) # Set the distance vectors - S_0k[s_id] = np.asarray(sensor['relative_location']) + S_0k[s_id] = np.asarray(sensor['relative_location'],dtype='str') # Calculate magnitudes - S_0k_mag[s_id] = np.sqrt(S_0k[s_id][0]**2 - + S_0k[s_id][1]**2 - + S_0k[s_id][2]**2) + S_0k_mag[s_id] = (Decimal(S_0k[s_id][0]) * Decimal(S_0k[s_id][0]) + + Decimal(S_0k[s_id][1]) * Decimal(S_0k[s_id][1]) + + Decimal(S_0k[s_id][2]) * Decimal(S_0k[s_id][2])).sqrt() # Calculate cos of angles between vector S0k and (x,y,z)-axis - cos_eta_0k[s_id] = S_0k[s_id] / S_0k_mag[s_id] + tmp = [] + for ii in range(0,3): + tmp.append(Decimal(S_0k[s_id][ii]) / S_0k_mag[s_id]) + cos_eta_0k[s_id] = tmp # Determine valid range for lag between signal of sensors 0 # and signal of sensor k from Eq. (43) in Shen et al. (2005) max_t_k[s_id] = S_0k_mag[s_id] / (COEFF_0 * V_GAS) @@ -100,6 +133,7 @@ def main(): # Create a H5-file reader reader = H5Reader(path / 'binary_signal.h5') # Read the time vector + t_signal_str = (reader.getDataSet('time')[:]).astype(str) t_signal = reader.getDataSet('time')[:] # Read the signal time series ds_signal = reader.getDataSet('signal') @@ -115,6 +149,10 @@ def main(): # of significance t_signal_ifd = t_signal[0:(len(t_signal)-1)] \ + (t_signal[1:len(t_signal)]-t_signal[0:(len(t_signal)-1)])/2.0 + t_signal_ifd_list = [] + for ii in range(0,len(t_signal_ifd)): + t_signal_ifd_list.append(str(Decimal(t_signal_str[ii]) \ + + (Decimal(t_signal_str[ii+1])-Decimal(t_signal_str[ii]))/Decimal(2.0))) # Interface-pairing signal-processing scheme # Shen et al. (2005) @@ -131,10 +169,10 @@ def main(): for ii,idx_rise_0 in enumerate(signal_ifd_rise_0): # DataFrame with time of IFD of front (2h) and rear (2h+1) bubble # interface for each sensor - ifd_times = pd.DataFrame(index=sensor_ids, columns=['t_2h','t_2h+1'],dtype='float64') + ifd_times = pd.DataFrame(index=sensor_ids, columns=['t_2h','t_2h+1']) # Main Sensor 0 # Write t_2h of sensor 0 - ifd_times['t_2h'][0] = t_signal_ifd[idx_rise_0] + ifd_times['t_2h'][0] = t_signal_ifd_list[idx_rise_0] # Search for t_2h+1 # Indices of the falling IFD signals for main sensor 0 @@ -142,8 +180,7 @@ def main(): # Index of the falling IFD signal for main sensor 0 idx_fall_0 = min(signal_ifd_fall_0[signal_ifd_fall_0 > idx_rise_0]) # Write t_b+1 of sensor 0 - ifd_times['t_2h+1'][0] = t_signal_ifd[idx_fall_0] - + ifd_times['t_2h+1'][0] = t_signal_ifd_list[idx_fall_0] # Auxillary sensors k for k in aux_sensor_ids: idk = np.where(sensor_ids == k)[0][0] @@ -194,12 +231,12 @@ def main(): # interface for each sensor # Write t_2h of sensor k if not np.isnan(idx_rise_k): - ifd_times['t_2h'][k] = t_signal_ifd[idx_rise_k] + ifd_times['t_2h'][k] = t_signal_ifd_list[idx_rise_k] else: ifd_times['t_2h'][k] = np.nan # Write t_2h+1 of sensor k if not np.isnan(idx_fall_k): - ifd_times['t_2h+1'][k] = t_signal_ifd[idx_fall_k] + ifd_times['t_2h+1'][k] = t_signal_ifd_list[idx_fall_k] else: ifd_times['t_2h+1'][k] = np.nan @@ -208,12 +245,12 @@ def main(): # and signal of sensor k from Eq. (42) in Shen et al. (2005) for k in aux_sensor_ids: # first interface (2h) - lag_2h = abs(ifd_times['t_2h'][0] - ifd_times['t_2h'][k]) + lag_2h = abs(Decimal(ifd_times['t_2h'][0]) - Decimal(ifd_times['t_2h'][k])) if lag_2h > max_t_k[k]: # lag too large ifd_times['t_2h'][k] = np.nan # second interface (2h+1) - lag_2hp1 = abs(ifd_times['t_2h+1'][0] - ifd_times['t_2h+1'][k]) + lag_2hp1 = abs(Decimal(ifd_times['t_2h+1'][0]) - Decimal(ifd_times['t_2h+1'][k])) if lag_2hp1 > max_t_k[k]: ifd_times['t_2h+1'][k] = np.nan @@ -233,7 +270,7 @@ def main(): print(f"Detected {len(bubbles_incomplete)} incomplete bubble signals.\n") print('\nRunning reconstruction algorithm....\n') # Initialize the Interfacial Area Concentration (IAC) - iac = 0.0 + iac = Decimal(0.0) # Loop over complete bubbles for ii,bubble_props in enumerate(bubbles_complete): ifd_times = bubble_props['ifd_times'] @@ -242,8 +279,7 @@ def main(): delta_t_0k_2hp1 = {} # Time differences k=1,2,3,4 for interfaces 2h and 2h+1 delta_t_k_h = {} - delta_t_k_h[0] = ifd_times['t_2h+1'][0] - ifd_times['t_2h'][0] - + delta_t_k_h[0] = Decimal(ifd_times['t_2h+1'][0]) - Decimal(ifd_times['t_2h'][0]) # 1. Reconstruct instantaneous local velocity # Inverse of measurable displacement velocity of the h-th bubble # for sensor 0 and k @@ -264,17 +300,17 @@ def main(): for k in aux_sensor_ids: # Time differences # Eq. (2) from Shen and Nakamura (2014) - delta_t_0k_2h[k] = ifd_times['t_2h'][k] \ - - ifd_times['t_2h'][0] - delta_t_0k_2hp1[k] = ifd_times['t_2h+1'][k] \ - - ifd_times['t_2h+1'][0] + delta_t_0k_2h[k] = Decimal(ifd_times['t_2h'][k]) \ + - Decimal(ifd_times['t_2h'][0]) + delta_t_0k_2hp1[k] = Decimal(ifd_times['t_2h+1'][k]) \ + - Decimal(ifd_times['t_2h+1'][0]) # Eq. (3) from Shen and Nakamura (2014) - delta_t_k_h[k] = ifd_times['t_2h+1'][k] \ - - ifd_times['t_2h'][k] + delta_t_k_h[k] = Decimal(ifd_times['t_2h+1'][k]) \ + - Decimal(ifd_times['t_2h'][k]) # Measurable displacement velocity # Eq. (30) from Shen and Nakamura (2014) V_m0k_h_inv[k] = (delta_t_0k_2h[k] + delta_t_0k_2hp1[k]) \ - / (2.0 * S_0k_mag[k]) + / (Decimal(2.0) * Decimal(S_0k_mag[k])) A_0.append(cos_eta_0k[k]) A_01_h.append(np.array([V_m0k_h_inv[k], cos_eta_0k[k][1], cos_eta_0k[k][2]])) A_02_h.append(np.array([cos_eta_0k[k][0], V_m0k_h_inv[k], cos_eta_0k[k][2]])) @@ -282,13 +318,12 @@ def main(): # Calculate matrix determinants # Eq. (32) to (35) from Shen and Nakamura (2014) - A_0_det = linalg.det(A_0) - A_01_h_det = linalg.det(A_01_h) - A_02_h_det = linalg.det(A_02_h) - A_03_h_det = linalg.det(A_03_h) - + A_0_det = calc_det(A_0) + A_01_h_det = calc_det(A_01_h) + A_02_h_det = calc_det(A_02_h) + A_03_h_det = calc_det(A_03_h) # Calculate instantaneous velocity magnitude and components - dummy = np.sqrt(A_01_h_det**2 + A_02_h_det**2 + A_03_h_det**2) + dummy = (A_01_h_det* A_01_h_det + A_02_h_det*A_02_h_det + A_03_h_det*A_03_h_det).sqrt() # Eq. (37) from Shen and Nakamura (2014) V_bh_mag = abs(A_0_det) * inverse_den(dummy) # Eq. (36) from Shen and Nakamura (2014) @@ -310,10 +345,10 @@ def main(): for k in aux_sensor_ids: # Chord lengths of interfaces 2h and 2h+1 # Eq. (43) and (44) in Shen and Nakamura (2014) - E_0k_2h[k] = abs(V_bh_mag)**2 * delta_t_0k_2h[k] \ + E_0k_2h[k] = V_bh_mag*V_bh_mag * delta_t_0k_2h[k] \ * (delta_t_0k_2hp1[k] + delta_t_k_h[0]) / S_0k_mag[k] \ - S_0k_mag[k] - E_0k_2hp1[k] = abs(V_bh_mag)**2 * delta_t_0k_2hp1[k] \ + E_0k_2hp1[k] = V_bh_mag*V_bh_mag * delta_t_0k_2hp1[k] \ * (delta_t_0k_2h[k] - delta_t_k_h[0]) / S_0k_mag[k] \ - S_0k_mag[k] # Eq. (46) and (48) in Shen and Nakamura (2014) @@ -325,17 +360,20 @@ def main(): B_03_2hp1.append(np.array([cos_eta_0k[k][0], cos_eta_0k[k][1], E_0k_2hp1[k]])) # Calculate matrix determinants # Eq. (46) and (48) in Shen and Nakamura (2014) - B_01_2h_det = linalg.det(B_01_2h) - B_02_2h_det = linalg.det(B_02_2h) - B_03_2h_det = linalg.det(B_03_2h) - B_01_2hp1_det = linalg.det(B_01_2hp1) - B_02_2hp1_det = linalg.det(B_02_2hp1) - B_03_2hp1_det = linalg.det(B_03_2hp1) + B_01_2h_det = calc_det(B_01_2h) + B_02_2h_det = calc_det(B_02_2h) + B_03_2h_det = calc_det(B_03_2h) + B_01_2hp1_det = calc_det(B_01_2hp1) + B_02_2hp1_det = calc_det(B_02_2hp1) + B_03_2hp1_det = calc_det(B_03_2hp1) # Calculate the bubble diameter - dummy_2h = np.sqrt(B_01_2h_det**2 + B_02_2h_det**2 + B_03_2h_det**2) - dummy_2hp1 = np.sqrt(B_01_2hp1_det**2 + B_02_2hp1_det**2 \ - + B_03_2hp1_det**2) + dummy_2h = (B_01_2h_det*B_01_2h_det \ + + B_02_2h_det*B_02_2h_det \ + + B_03_2h_det*B_03_2h_det).sqrt() + dummy_2hp1 = (B_01_2hp1_det*B_01_2hp1_det \ + + B_02_2hp1_det*B_02_2hp1_det \ + + B_03_2hp1_det*B_03_2hp1_det).sqrt() # Eq. (50) from Shen and Nakamura (2014) D_h_2h = dummy_2h / A_0_det D_h_2hp1 = dummy_2hp1 / A_0_det @@ -343,12 +381,12 @@ def main(): # Calculate interfacial normal unit vectors # Eq. (51) in Shen and Nakamura (2014) - cos_eta_i_2h = np.ones(3) * A_0_det \ + cos_eta_i_2h = np.array([Decimal(1.0), Decimal(1.0), Decimal(1.0)]) * A_0_det \ * inverse_den(abs(A_0_det)) \ * inverse_den(dummy_2h) cos_eta_i_2h = cos_eta_i_2h \ * np.array([B_01_2h_det, B_02_2h_det, B_03_2h_det]) - cos_eta_i_2hp1 = np.ones(3) * A_0_det \ + cos_eta_i_2hp1 = np.array([Decimal(1.0), Decimal(1.0), Decimal(1.0)]) * A_0_det \ * inverse_den(abs(A_0_det)) \ * inverse_den(dummy_2hp1) cos_eta_i_2hp1 = cos_eta_i_2hp1 \ @@ -359,14 +397,16 @@ def main(): # Calculate the local Interfacial Area Concentration (IAC) # Eq. (54) in Shen and Nakamura (2014) # Sampling duration Omega - Omega = max(t_signal) - min(t_signal) - dummy_A = A_01_h_det**2 + A_02_h_det**2 + A_03_h_det**2 - dummy_B = np.sqrt(B_01_2h_det**2 + B_02_2h_det**2 + B_03_2h_det**2) + Omega = Decimal(t_signal[-1]) - Decimal(t_signal[0]) + dummy_A = A_01_h_det*A_01_h_det + A_02_h_det*A_02_h_det + A_03_h_det*A_03_h_det + dummy_B = (B_01_2h_det*B_01_2h_det \ + + B_02_2h_det*B_02_2h_det \ + + B_03_2h_det*B_03_2h_det).sqrt() dummy_AB = abs(A_0_det * (A_01_h_det*B_01_2h_det \ + A_02_h_det*B_02_2h_det \ + A_03_h_det*B_03_2h_det)) iac_h = dummy_A * dummy_B * inverse_den(dummy_AB) - iac += 2.0 / Omega * iac_h + iac += Decimal(2.0) / Omega * iac_h # Display progress if args.progress: printProgressBar(ii, len(bubbles_complete), prefix = 'Progress:', suffix = 'Complete', length = 50) @@ -436,7 +476,7 @@ def main(): # Add the attributes ds_d.attrs['labels'] = ['D_h_2h', 'D_h_2hp1'] # Create the IAC and void_fraction datasets - writer.writeDataSet('IAC', np.array([iac]), 'float64') + writer.writeDataSet('IAC', np.array([iac],dtype='float64'), 'float64') writer.writeDataSet('voidFraction', np.array([np.mean(signal)]), 'float64') writer.close() time2 = time.time()