diff --git a/postprocessing_h5py/postprocessing_common_h5py.py b/postprocessing_h5py/postprocessing_common_h5py.py index 2c237c59..abe26955 100755 --- a/postprocessing_h5py/postprocessing_common_h5py.py +++ b/postprocessing_h5py/postprocessing_common_h5py.py @@ -19,6 +19,9 @@ """ This script contains a number of helper functions to create visualizations outsside of fenics. + +All functions authored by David Bruneau unless specified + """ def read_command_line(): @@ -137,6 +140,7 @@ def get_interface_ids(meshFile): def get_sampling_constants(df,start_t,end_t): ''' + Author: Daniel Macdonald T = period, in seconds, nsamples = samples per cycle fs = sample rate @@ -162,6 +166,10 @@ def read_csv_file(filepath): return df def filter_SPI(U, W_low_cut, tag): + """ + Author: Mehdi Najafi + + """ if tag=="withmean": U_fft = fft(U) else: @@ -185,6 +193,10 @@ def filter_SPI(U, W_low_cut, tag): def filter_SPI_print(U, W_low_cut, tag): + """ + Author: Mehdi Najafi + + """ if tag=="withmean": U_fft = fft(U) else: @@ -208,6 +220,10 @@ def filter_SPI_print(U, W_low_cut, tag): return Power_25Hz/Power_0Hz def calculate_spi(case_name, df, output_folder, meshFile,start_t,end_t,low_cut,high_cut, dvp): + """ + Author: Mehdi Najafi and David Bruneau + + """ # cut, thresh # Get wall and fluid ids fluidIDs, wallIDs, allIDs = get_domain_ids(meshFile) diff --git a/postprocessing_h5py/spectrograms.py b/postprocessing_h5py/spectrograms.py index ad233084..55ad3c65 100755 --- a/postprocessing_h5py/spectrograms.py +++ b/postprocessing_h5py/spectrograms.py @@ -27,7 +27,7 @@ except: print("Could not import pyvista and/or vtk, install these packages or use 'RandomPoint' sampling for spectrograms.") """ -This library contains helper functions for creating spectrograms. +This library contains helper functions for creating spectrograms. All functions authored by David Bruneau unless otherwise specified. """ def read_spec_config(config_file,dvp): @@ -210,19 +210,8 @@ def read_spectrogram_data(case_path, mesh_name, save_deg, stride, start_t, end_t elif sampling_method == "Spatial": # This method only works with a specified sphere mesh, surf = postprocessing_common_pv.assemble_mesh(mesh_path) - #bounds=surf.bounds bounds = [x_sphere-r_sphere, x_sphere+r_sphere, y_sphere-r_sphere, y_sphere+r_sphere, z_sphere-r_sphere, z_sphere+r_sphere] - #def generate_points(bounds, subdivisions=50): - # x_points=np.linspace(bounds[0], bounds[1],num=subdivisions) - # y_points=np.linspace(bounds[2], bounds[3],num=subdivisions) - # z_points=np.linspace(bounds[4], bounds[5],num=subdivisions) - # points = np.array([[x_points[0], y_points[0], z_points[0]],[x_points[0], y_points[0], z_points[1]]]) - # for i in range(subdivisions): - # for j in range(subdivisions): - # for k in range(subdivisions): - # points=np.append(points,[[x_points[i], y_points[j], z_points[k]]], axis=0) - # return points[2:,:] - + def generate_points(bounds, subdivisions=50): x_points=np.linspace(bounds[0], bounds[1],num=subdivisions) @@ -254,10 +243,6 @@ def generate_points(bounds, subdivisions=50): t = time.time() print("shaped point cloud") print(t - start) - #coords_sphere = coords[region_ids] - #print(coords_sphere) - #print("Length of coords array") - #print(coords.shape) tree = KDTree(coords) _, idx = tree.query(sphere_sel.points) #find closest node to the points in the equispaced points @@ -266,21 +251,10 @@ def generate_points(bounds, subdivisions=50): print("queried kdtree") print(t - start) - #print("Unsampled equispaced data") - #print(len(idx)) - #print(idx) - # + # Make sure this is the correct order equispaced_fluid_ids = [x for x in idx if x in region_ids] - #equispaced_fluid_ids=list(set(region_ids).intersection(idx)) - #print("Unsampled equispaced fluid data") - #print(len(equispaced_fluid_ids)) - #print(equispaced_fluid_ids) - idx_sampled = np.random.choice(equispaced_fluid_ids, n_samples) - #print("Sampled equispaced fluid data") - #print(len(idx_sampled)) - #print(idx_sampled) idx_sampled_set = set(idx_sampled) # compare the length and print if the list contains duplicates @@ -290,17 +264,13 @@ def generate_points(bounds, subdivisions=50): print("No duplicates found in the list") #sampled_point_cloud=pv.PolyData(coords[idx_sampled]) -## #plotter= pv.Plotter(off_screen=True) - # + #plotter.add_mesh(surf, # color='red', # show_scalar_bar=False, # opacity=0.05) -### #plotter.add_points(sampled_point_cloud, render_points_as_spheres=True, point_size=0.03) - # - # #plotter.show(auto_close=False) #plotter.show(screenshot=imageFolder + "/points"+str(idx_sampled[0])+".png") @@ -325,32 +295,6 @@ def get_location_from_pointID(probeID,polydata): return nearestPointLoc -#def find_nearest_points_in_radius(probe_pos,pl,num_points,polyData,radius): -# """ -# probeID is actual ID of nearest point -# nearestPointLoc is xyz location of probeID -# nearestNeighbours is all the points within the radius -# nearest_points is sample of nearest_Neighbours -# NOTE: nearest_points is likely unneccesary after implementing masking -# function. -# """ -# probeID = pl.FindClosestPoint(probe_pos) -# nearestPointLoc = [10000,10000,10000] -# polyData.GetPoint(probeID,nearestPointLoc) -# nearestNeighbours = vtk.vtkIdList() -# pl.FindPointsWithinRadius(radius,probe_pos,nearestNeighbours) -# num_Ids = nearestNeighbours.GetNumberOfIds() -# print('There are',num_Ids,'in radius',radius) -# idList = [] -# if num_Ids >= 1: -# if num_points > num_Ids: -# num_points = num_Ids -# random_index = random.sample(range(num_Ids), num_points-1) -# idList = [nearestNeighbours.GetId(idx) for idx in random_index] -# nearest_points = [probeID] + idList -# return nearest_points - - def find_points_in_sphere(cent,rad,coords): # Calculate vector from center to each node in the mesh @@ -371,42 +315,9 @@ def find_points_in_sphere(cent,rad,coords): return points_in_sphere - -#def vtk_read_unstructured_grid(FILE): -# ''' -# Reads vtu file. Create a more general reader, particularly for h5 mesh -# going forward. -# ''' -# reader = vtk.vtkXMLUnstructuredGridReader() -# reader.SetFileName(FILE) -# reader.Update() -# return reader.GetOutput() - -#def vtk_mask_points(polyData, max_points=1000): -# ''' -# Randomly mask points. -# SetRandomModeType(2) uses stratified spatial sampling. -# ''' -# ptMask = vtk.vtkMaskPoints() -# ptMask.SetInputData(polyData) -# ptMask.RandomModeOn() -# ptMask.SetRandomModeType(2) -# ptMask.SetMaximumNumberOfPoints(max_points) -# ptMask.Update() -# return ptMask.GetOutput() - -#def get_original_ids(masked_polydata, polyData): -# n = masked_polydata.GetPoints().GetNumberOfPoints() -# original_ids = np.zeros(n) -# pll = vtk.vtkPointLocator() -# pll.SetDataSet(polyData) -# for pt in range(n): -# location = masked_polydata.GetPoints().GetPoint(pt) -# original_ids[pt] = pll.FindClosestPoint(location) -# return original_ids - def shift_bit_length(x): ''' + Author: Dan MacDonald round up to nearest pwr of 2 https://stackoverflow.com/questions/14267555/find-the-smallest-power-of-2-greater-than-n-in-python ''' @@ -433,6 +344,7 @@ def get_psd(dfNearest,fsamp,scaling="density"): def get_spectrogram(dfNearest,fsamp,nWindow,overlapFrac,window,start_t,end_t, scaling='spectrum', interpolate = False): ''' + Author: Daniel Macdonald Calculates spectrogram input dfNearest is a pandas df of shape (num_points, num_timesteps) fsamp is sampling frequency @@ -475,20 +387,8 @@ def get_spectrogram(dfNearest,fsamp,nWindow,overlapFrac,window,start_t,end_t, sc Pxx_mean[Pxx_mean<0] = 1e-16 return Pxx_mean, freqs, bins -def spectrogram_scaling(Pxx_mean,thresh_percent): - Pxx_scaled = np.log(Pxx_mean) - max_val = np.max(Pxx_scaled) - min_val = np.min(Pxx_scaled) - thresh = max_val-(max_val-min_val)*thresh_percent - print('Pxx_scaled max', max_val) - print('Pxx_scaled max', min_val) - print('Pxx threshold', thresh) - - Pxx_threshold_indices = Pxx_scaled < thresh - Pxx_scaled[Pxx_threshold_indices] = thresh - return Pxx_scaled, max_val, min_val, thresh - -def spectrogram_scaling_old(Pxx_mean,lower_thresh): +def spectrogram_scaling(Pxx_mean,lower_thresh): + # Author: Daniel Macdonald Pxx_scaled = np.log(Pxx_mean) max_val = np.max(Pxx_scaled) min_val = np.min(Pxx_scaled) @@ -523,25 +423,26 @@ def butter_bandpass(lowcut, highcut, fs, order=5, btype='band'): return b, a def butter_bandpass_filter(data, lowcut=25.0, highcut=15000.0, fs=2500.0, order=5, btype='band'): + # Author: Daniel Macdonald b, a = butter_bandpass(lowcut, highcut, fs, order=order,btype=btype) y = filtfilt(b, a, data) return y def filter_time_data(df,fs,lowcut=25.0,highcut=15000.0,order=6,btype='highpass'): + # Author: Daniel Macdonald df_filtered = df.copy() for row in range(df.shape[0]): df_filtered.iloc[row] = butter_bandpass_filter(df.iloc[row],lowcut=lowcut,highcut=highcut,fs=fs,order=order,btype=btype) return df_filtered def compute_average_spectrogram(df, fs, nWindow,overlapFrac,window,start_t,end_t,thresh, scaling="spectrum", filter_data=False,thresh_method="new"): + # Author: Daniel Macdonald if filter_data == True: df = filter_time_data(df,fs) Pxx_mean, freqs, bins = get_spectrogram(df,fs,nWindow,overlapFrac,window,start_t,end_t, scaling) # mode of the spectrogram - if thresh_method == "new": + if thresh_method == "old": Pxx_scaled, max_val, min_val, lower_thresh = spectrogram_scaling(Pxx_mean,thresh) - elif thresh_method == "old": - Pxx_scaled, max_val, min_val, lower_thresh = spectrogram_scaling_old(Pxx_mean,thresh) elif thresh_method == "log_only": Pxx_scaled = np.log(Pxx_mean) max_val = np.max(Pxx_scaled) @@ -557,8 +458,6 @@ def compute_average_spectrogram(df, fs, nWindow,overlapFrac,window,start_t,end_t return bins, freqs, Pxx_scaled, max_val, min_val, lower_thresh def plot_spectrogram(fig1,ax1,bins,freqs,Pxx,ylim_,title=None,convert_a=0.0,convert_b=0.0,x_label=None,color_range=None): - #plt.figure(figsize=(14,7)) #fig size same as before - if color_range == None: im = ax1.pcolormesh(bins, freqs, Pxx, shading = 'gouraud')#, vmin = -30, vmax = -4) # look up in matplotlib to add colorbar @@ -586,37 +485,8 @@ def time_convert(x): ax2.set_xlabel(x_label) -def plot_spectrogram_1ax(bins,freqs,Pxx,ylim_,title=None,path=None,convert_a=1.0,convert_b=0.0,x_label=None,color_range=None): - #plt.figure(figsize=(14,7)) #fig size same as before - - bins = bins*convert_a+convert_b - fig1, ax1 = plt.subplots() - if color_range == None: - im = ax1.pcolormesh(bins, freqs, Pxx, shading = 'gouraud')#, vmin = -30, vmax = -4) # look up in matplotlib to add colorbar - else: - im = ax1.pcolormesh(bins, freqs, Pxx, shading = 'gouraud', vmin = color_range[0], vmax = color_range[1]) # look up in matplotlib to add colorbar - - - fig1.colorbar(im, ax=ax1) - #fig1.set_size_inches(8, 6) #fig1.set_size_inches(10, 7) - - ax1.set_xlabel(x_label) - ax1.set_ylabel('Frequency [Hz]') - - ax1.set_ylim([0,ylim_]) # new - ax1.set_title('{}'.format(title),y=1.08) - - if path != None: - fig1.savefig(path) - path_csv = path.replace(".png",".csv") - #freqs_txt = np.array2string(freqs, precision=2, separator=',',) - data_csv = np.append(freqs[np.newaxis].T,Pxx, axis=1) - bins_txt = "freq (Hz)/bins (s)," + np.array2string(bins, precision=2, separator=',',).replace("[","").replace("]","") - np.savetxt(path_csv, data_csv,header=bins_txt, delimiter=",") - - def chromagram_from_spectrogram(Pxx,fs,n_fft,n_chroma=24,norm=True): - + # Author: Daniel Macdonald # Calculate chroma filterbank chromafb = chroma_filterbank( sr=fs, @@ -639,6 +509,7 @@ def chromagram_from_spectrogram(Pxx,fs,n_fft,n_chroma=24,norm=True): return chroma def calc_chroma_entropy(chroma,n_chroma): + # Author: Daniel Macdonald chroma_entropy = -np.sum(chroma * np.log(chroma), axis=0) / np.log(n_chroma) #print(np.sum(chroma,axis=0)) #chroma_entropy = entropy(normed_power) / np.log(n_chroma) @@ -672,6 +543,7 @@ def plot_chromagram(fig1,ax1,bins,chroma,title=None,path=None,convert_a=1.0,conv def get_sampling_constants(df,start_t,end_t): + # Author: Daniel Macdonald ''' T = period, in seconds, nsamples = samples per cycle @@ -682,23 +554,8 @@ def get_sampling_constants(df,start_t,end_t): fs = nsamples/T return T, nsamples, fs -#def sample_polydata(df, polydata_geo, num_points=2000): -# polydata_masked = vtk_mask_points(polydata_geo, num_points) -# #print(polydata_masked) -# ''' -# Find the indices of the df that we'll use, filter df -# ''' -# sampled_points_og_ids = get_original_ids(polydata_masked, polydata_geo) -# df_sampled = df.filter(sampled_points_og_ids.astype(int), axis = 0) -# ''' -# Reset index here allows us to easily search the masked polydata -# ''' -# df_sampled.reset_index(level=0, inplace=True) -# df_sampled.index.names = ['NewIds'] -# df_sampled.drop('Ids', 1, inplace = True) -# return df_sampled, polydata_masked - -# The following functions are taken from librosa to generate chroma filterbank ############################################################ + +# The remaining functions are taken from librosa to generate chroma filterbank def octs_to_hz(octs, tuning=0.0, bins_per_octave=12): """Convert octaves numbers to frequencies. diff --git a/src/fsipy/simulations/cylinder.py b/src/fsipy/simulations/cylinder.py index 89f84f4b..bffd27bb 100644 --- a/src/fsipy/simulations/cylinder.py +++ b/src/fsipy/simulations/cylinder.py @@ -1,3 +1,194 @@ """ Propblem file for tiny cylinder FSI simulation """ +import numpy as np +from turtleFSI.problems import * +from dolfin import HDF5File, Mesh, MeshFunction, assemble, UserExpression, FacetNormal, ds, \ + DirichletBC, Measure, inner, parameters + +# set compiler arguments +parameters["form_compiler"]["quadrature_degree"] = 6 +parameters["form_compiler"]["optimize"] = True +# The "ghost_mode" has to do with the assembly of form containing the facet +# normals n('+') within interior boundaries (dS). for 3D mesh the value should +# be "shared_vertex", for 2D mesh "shared_facet", the default value is "none" +parameters["ghost_mode"] = "shared_vertex" +_compiler_parameters = dict(parameters["form_compiler"]) + + +def set_problem_parameters(default_variables, **namespace): + # Overwrite default values + E_s_val = 1e6 # Young modulus (Pa) + nu_s_val = 0.45 + mu_s_val = E_s_val / (2 * (1 + nu_s_val)) # 0.345E6 + lambda_s_val = nu_s_val * 2.0 * mu_s_val / (1.0 - 2.0 * nu_s_val) + + default_variables.update( + dict( + # Temporal parameters + T=0.1, # Simulation end time + dt=0.001, # Time step size + theta=0.501, # Theta scheme (implicit/explicit time stepping) + save_step=1, # Save frequency of files for visualisation + checkpoint_step=50, # Save frequency of checkpoint files + # Linear solver parameters + linear_solver="mumps", + atol=1e-6, # Absolute tolerance in the Newton solver + rtol=1e-6, # Relative tolerance in the Newton solver + recompute=20, # Recompute the Jacobian matix within time steps + recompute_tstep=20, # Recompute the Jacobian matix over time steps + # boundary condition parameters + mesh_path="mesh/cylinder.h5", + inlet_id=2, # inlet id for the fluid + inlet_outlet_s_id=11, # inlet and outlet id for solid + fsi_id=22, # id for fsi interface + rigid_id=11, # "rigid wall" id for the fluid and mesh problem + outer_wall_id=33, # id for the outer surface of the solid + # Fluid parameters + rho_f=1.025e3, # Fluid density [kg/m3] + mu_f=3.5e-3, # Fluid dynamic viscosity [Pa.s] + dx_f_id=1, # ID of marker in the fluid domain. When reading the mesh, the fuid domain is assigned with a 1. + v_max_final=0.75, # Steady state, maximum centerline velocity of parabolic profile + P_final=10000, # Steady State pressure applied to wall + # Solid parameters + rho_s=1.0e3, # Solid density [kg/m3] + mu_s=mu_s_val, # Solid shear modulus or 2nd Lame Coef. [Pa] + nu_s=nu_s_val, # Solid Poisson ratio [-] + lambda_s=lambda_s_val, # Solid 1st Lame Coef. [Pa] + dx_s_id=2, # ID of marker in the solid domain + # mesh lifting parameters (see turtleFSI for options) + extrapolation="laplace", # laplace, elastic, biharmonic, no-extrapolation + extrapolation_sub_type="constant", # ["constant","small_constant","volume","volume_change","bc1","bc2"] + folder="cylinder_results", # output folder generated for simulation + save_deg=1 # save_deg=1 saves corner nodes only, save_deg=2 saves corner + mid-point nodes for viz + ) + ) + + return default_variables + + +def get_mesh_domain_and_boundaries(mesh_path, folder, **namespace): + print("Obtaining mesh, domains and boundaries...") + mesh = Mesh() + hdf = HDF5File(mesh.mpi_comm(), mesh_path, "r") + hdf.read(mesh, "/mesh", False) + boundaries = MeshFunction("size_t", mesh, 2) + hdf.read(boundaries, "/boundaries") + domains = MeshFunction("size_t", mesh, 3) + hdf.read(domains, "/domains") + + return mesh, domains, boundaries + + +class VelInPara(UserExpression): + def __init__(self, t, t_ramp, v_max_final, n, dsi, mesh, **kwargs): + self.t = t + self.t_ramp = t_ramp + self.v_max_final = v_max_final + self.v = 0.0 + self.n = n + self.dsi = dsi + self.d = mesh.geometry().dim() + self.x = SpatialCoordinate(mesh) + # Compute area of boundary tesselation by integrating 1.0 over all facets + self.A = assemble(Constant(1.0, name="one") * self.dsi) + # Compute barycenter by integrating x components over all facets + self.c = [assemble(self.x[i] * self.dsi) / self.A for i in range(self.d)] + # Compute radius by taking max radius of boundary points + self.r = np.sqrt(self.A / np.pi) + super().__init__(**kwargs) + + def update(self, t): + self.t = t + # apply a sigmoid ramp to the pressure + if self.t < self.t_ramp: + ramp_factor = -0.5 * np.cos(np.pi * self.t / self.t_ramp) + 0.5 + else: + ramp_factor = 1.0 + self.v = ramp_factor * self.v_max_final + + if MPI.rank(MPI.comm_world) == 0: + print("v (centerline, at inlet) = {} m/s".format(self.v)) + + def eval(self, value, x): + r2 = ( + (x[0] - self.c[0]) ** 2 + (x[1] - self.c[1]) ** 2 + (x[2] - self.c[2]) ** 2 + ) # radius**2 + fact_r = 1 - (r2 / self.r**2) + + value[0] = -self.n[0] * (self.v) * fact_r # *self.t + value[1] = -self.n[1] * (self.v) * fact_r # *self.t + value[2] = -self.n[2] * (self.v) * fact_r # *self.t + + def value_shape(self): + return (3,) + + +class InnerP(UserExpression): + def __init__(self, t, t_ramp, P_final, **kwargs): + self.t = t + self.t_ramp = t_ramp + self.P_final = P_final + self.P = 0.0 + super().__init__(**kwargs) + + def update(self, t): + self.t = t + # apply a sigmoid ramp to the pressure + if self.t < self.t_ramp: + ramp_factor = -0.5 * np.cos(np.pi * self.t / self.t_ramp) + 0.5 + else: + ramp_factor = 1.0 + self.P = ramp_factor * self.P_final + + if MPI.rank(MPI.comm_world) == 0: + print("P = {} Pa".format(self.P)) + + def eval(self, value, x): + value[0] = self.P + + def value_shape(self): + return () + + +def create_bcs(dvp_, DVP, mesh, boundaries, P_final, v_max_final, fsi_id, inlet_id, + inlet_outlet_s_id, rigid_id, psi, F_solid_linear, **namespace): + + # Apply pressure at the fsi interface by modifying the variational form + p_out_bc_val = InnerP(t=0.0, t_ramp=0.1, P_final=P_final, degree=2) + dSS = Measure("dS", domain=mesh, subdomain_data=boundaries) + n = FacetNormal(mesh) + # defined on the reference domain + F_solid_linear += (p_out_bc_val * inner(n("+"), psi("+")) * dSS(fsi_id)) + + # Fluid velocity BCs + dsi = ds(inlet_id, domain=mesh, subdomain_data=boundaries) + n = FacetNormal(mesh) + ndim = mesh.geometry().dim() + ni = np.array([assemble(n[i] * dsi) for i in range(ndim)]) + n_len = np.sqrt(sum([ni[i] ** 2 for i in range(ndim)])) # Should always be 1!? + normal = ni / n_len + + # Parabolic Inlet Velocity Profile + u_inflow_exp = VelInPara(t=0.0, t_ramp=0.1, v_max_final=v_max_final, + n=normal, dsi=dsi, mesh=mesh, degree=3) + u_inlet = DirichletBC(DVP.sub(1), u_inflow_exp, boundaries, inlet_id) + u_inlet_s = DirichletBC(DVP.sub(1), ((0.0, 0.0, 0.0)), boundaries, inlet_outlet_s_id) + + # Solid Displacement BCs + d_inlet = DirichletBC(DVP.sub(0), ((0.0, 0.0, 0.0)), boundaries, inlet_id) + d_inlet_s = DirichletBC(DVP.sub(0), ((0.0, 0.0, 0.0)), boundaries, inlet_outlet_s_id) + d_rigid = DirichletBC(DVP.sub(0), ((0.0, 0.0, 0.0)), boundaries, rigid_id) + + # Assemble boundary conditions + bcs = [u_inlet, d_inlet, u_inlet_s, d_inlet_s, d_rigid] + + return dict(bcs=bcs, u_inflow_exp=u_inflow_exp, p_out_bc_val=p_out_bc_val, + F_solid_linear=F_solid_linear) + + +def pre_solve(t, u_inflow_exp, p_out_bc_val, **namespace): + # Update the time variable used for the inlet boundary condition + u_inflow_exp.update(t) + p_out_bc_val.update(t) + return dict(u_inflow_exp=u_inflow_exp, p_out_bc_val=p_out_bc_val)