diff --git a/modules/createEVENT/IsolatedBuildingCFD/CMakeLists.txt b/modules/createEVENT/IsolatedBuildingCFD/CMakeLists.txt index 2b2b53d9a..5047cc44e 100644 --- a/modules/createEVENT/IsolatedBuildingCFD/CMakeLists.txt +++ b/modules/createEVENT/IsolatedBuildingCFD/CMakeLists.txt @@ -1,7 +1,7 @@ #Python scripts for pre and post processing simcenter_add_python_script(SCRIPT foam_dict_reader.py) simcenter_add_python_script(SCRIPT setup_case.py) -simcenter_add_python_script(SCRIPT process_output_data.py) +simcenter_add_python_script(SCRIPT post_process_output.py) simcenter_add_python_script(SCRIPT IsolatedBuildingCFD.py) add_subdirectory(templateOF10Dicts) diff --git a/modules/createEVENT/IsolatedBuildingCFD/post_process_output.py b/modules/createEVENT/IsolatedBuildingCFD/post_process_output.py new file mode 100644 index 000000000..461187f04 --- /dev/null +++ b/modules/createEVENT/IsolatedBuildingCFD/post_process_output.py @@ -0,0 +1,452 @@ +# -*- coding: utf-8 -*- +# Copyright (c) 2016-2017, The Regents of the University of California (Regents). +# All rights reserved. +# +# Redistribution and use in source and binary forms, with or without +# modification, are permitted provided that the following conditions are met: +# +# 1. Redistributions of source code must retain the above copyright notice, this +# list of conditions and the following disclaimer. +# 2. Redistributions in binary form must reproduce the above copyright notice, +# this list of conditions and the following disclaimer in the documentation +# and/or other materials provided with the distribution. +# +# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR +# ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +# +# The views and conclusions contained in the software and documentation are those +# of the authors and should not be interpreted as representing official policies, +# either expressed or implied, of the FreeBSD Project. +# +# REGENTS SPECIFICALLY DISCLAIMS ANY WARRANTIES, INCLUDING, BUT NOT LIMITED TO, +# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. +# THE SOFTWARE AND ACCOMPANYING DOCUMENTATION, IF ANY, PROVIDED HEREUNDER IS +# PROVIDED "AS IS". REGENTS HAS NO OBLIGATION TO PROVIDE MAINTENANCE, SUPPORT, +# UPDATES, ENHANCEMENTS, OR MODIFICATIONS. + +# +# Contributors: +# Abiy Melaku + + +# +# This script reads OpenFOAM output and plot the characteristics of the +# approaching wind. For now, it read and plots only velocity field data and +# pressure on predicted set of probes. +# + +import sys +import os +import subprocess +import json +import stat +import shutil +from pathlib import Path +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +from scipy import signal +from scipy.interpolate import interp1d +from scipy.interpolate import UnivariateSpline +from scipy import stats +import pandas as pd +import plotly.graph_objects as go +from plotly.subplots import make_subplots +import argparse +import re + +def read_bin_forces(fileName): + """ + Reads binData measured at the center of each bin. + + Reads force data integrated over a surface from Open file and returns + bin heights, time, and the forces and moments vector on the bins for each + time step. + + """ + forces = [] + moments = [] + time = [] + nbins = 0 + + with open(fileName, "r") as f: + for line in f: + if line.startswith('#'): + #Read the origin where the force are integrated + if line.startswith('# bins'): + line = line.replace('# bins', '') + line = line.replace(':', '') + line = line.split() + nbins = int(line[0]) + coords = np.zeros((nbins, 3)) + elif line.startswith('# x co-ords'): + line = line.replace('# x co-ords', '') + line = line.replace(':', '') + line = line.split() + for i in range(nbins): + coords[i,0] = line[i] + elif line.startswith('# y co-ords'): + line = line.replace('# y co-ords', '') + line = line.replace(':', '') + line = line.split() + for i in range(nbins): + coords[i,1] = line[i] + elif line.startswith('# z co-ords'): + line = line.replace('# z co-ords', '') + line = line.replace(':', '') + line = line.split() + for i in range(nbins): + coords[i,2] = line[i] + else: + continue + + # Read only the pressure part + else: + line = line.replace('(','') + line = line.replace(')','') + line = line.split() + time.append(float(line[0])) + story_force = np.zeros((nbins, 3)) + story_moments = np.zeros((nbins, 3)) + for i in range(nbins): + start = 12*i + 1 + + #Take only pressure part + story_force[i,0] = line[start] + story_force[i,1] = line[start + 1] + story_force[i,2] = line[start + 2] + + story_moments[i,0] = line[start + 6] + story_moments[i,1] = line[start + 6 + 1] + story_moments[i,2] = line[start + 6 + 2] + + forces.append(story_force) + moments.append(story_moments) + + time = np.asarray(time, dtype=np.float32) + forces = np.asarray(forces, dtype=np.float32) + moments = np.asarray(moments, dtype=np.float32) + + return coords, time, forces, moments + +def read_forces(fileName): + """ + Reads force data integrated over a surface from OpenFOAM file and returns + origin(center of rotation), time, and the forces and moments vector for + each time step. + + """ + origin = np.zeros(3) + forces = [] + moments = [] + time = [] + + with open(fileName, "r") as f: + for line in f: + if line.startswith('#'): + #Read the origin where the force are integrated + if line.startswith('# CofR'): + line = line.replace('(','') + line = line.replace(')','') + line = line.split() + origin[0] = line[3] # x-coordinate + origin[1] = line[4] # y-coordinate + origin[2] = line[5] # z-coordinate + else: + continue + # Read only the pressure part of force and moments. + # Viscous and porous are ignored + else: + line = line.replace('(','') + line = line.replace(')','') + line = line.split() + time.append(float(line[0])) + forces.append([float(line[1]), float(line[2]), float(line[3])]) + moments.append([float(line[1 + 9]), float(line[2 + 9]), float(line[3 + 9])]) + + time = np.asarray(time, dtype=np.float32) + forces = np.asarray(forces, dtype=np.float32) + moments = np.asarray(moments, dtype=np.float32) + + return origin, time, forces, moments + + +def readPressureProbes(fileName): + """ + Created on Wed May 16 14:31:42 2018 + + Reads pressure probe data from OpenFOAM and return the probe location, time, and the pressure + for each time step. + + @author: Abiy + """ + probes = [] + p = [] + time = [] + + with open(fileName, "r") as f: + for line in f: + if line.startswith('#'): + if line.startswith('# Probe'): + line = line.replace('(','') + line = line.replace(')','') + line = line.split() + probes.append([float(line[3]),float(line[4]),float(line[5])]) + else: + continue + else: + line = line.split() + time.append(float(line[0])) + p_probe_i = np.zeros([len(probes)]) + for i in range(len(probes)): + p_probe_i[i] = float(line[i + 1]) + p.append(p_probe_i) + + probes = np.asarray(probes, dtype=np.float32) + time = np.asarray(time, dtype=np.float32) + p = np.asarray(p, dtype=np.float32) + + return probes, time, p + +def read_pressure_data(file_names): + """ + This functions takes names of different OpenFOAM pressure measurements and connect + them into one file removing overlaps if any. All the probes must be in the same + location, otherwise an error might show up. + + Parameters + ---------- + *args + List of file pashes of pressure data to be connected together. + + Returns + ------- + time, pressure + Returns the pressure time and pressure data of the connected file. + """ + no_files = len(file_names) + connected_time = [] # Connected array of time + connected_p = [] # connected array of pressure. + + time1 = [] + p1 = [] + time2 = [] + p2 = [] + probes= [] + + for i in range(no_files): + probes, time2, p2 = readPressureProbes(file_names[i]) + + if i==0: + connected_time = time2 + connected_p = p2 + else: + try: + index = np.where(time2 > time1[-1])[0][0] + # index += 1 + + except: + # sys.exit('Fatal Error!: the pressure filese have time gap') + index = 0 # Joint them even if they have a time gap + + connected_time = np.concatenate((connected_time, time2[index:])) + connected_p = np.concatenate((connected_p, p2[index:])) + + time1 = time2 + p1 = p2 + return probes, connected_time, connected_p + + +class PressureData: + """ + A class that holds a pressure data and performs the following operations: + - mean and rms pressure coefficients + - peak pressure coefficients + """ + def __init__(self, path, u_ref=0.0, rho=1.25, p_ref=0.0, + start_time=None, end_time=None): + self.path = path + self.u_ref = u_ref + self.p_ref = p_ref + self.rho = rho + self.start_time = start_time + self.end_time = end_time + self.read_cfd_data() + self.set_time() + self.Nt = len(self.time) + self.T = self.time[-1] + self.z = self.probes[:,2] + self.y = self.probes[:,1] + self.x = self.probes[:,0] + self.dt = np.mean(np.diff(self.time)) + self.probe_count = np.shape(self.probes)[0] + + def read_cfd_data (self): + if os.path.isdir(self.path): + print("Reading from path : %s" % (self.path)) + time_names = os.listdir(self.path) + sorted_index = np.argsort(np.float_(time_names)).tolist() + # print(sorted_index) + # print("\tTime directories: %s" %(time_names)) + file_names = [] + + for i in range(len(sorted_index)): + file_name = os.path.join(self.path, time_names[sorted_index[i]],'p') + file_names.append(file_name) + + # print(file_names) + self.probes, self.time, self.p = read_pressure_data(file_names) + self.p = self.rho*np.transpose(self.p) # OpenFOAM gives p/rho + p_dyn = 0.5*self.rho*self.u_ref**2.0 + + if self.u_ref != 0.0: + self.cp = self.p/p_dyn + + + # self.p = np.transpose(self.p) # OpenFOAM gives p/rho + else: + print("Cannot find the file path: %s" % (self.path)) + + + def set_time (self): + if(self.start_time != None): + start_index = int(np.argmax(self.time > self.start_time)) + self.time = self.time[start_index:] + # self.cp = self.cp[:,start_index:] + try: + self.p = self.p[:,start_index:] + self.cp = self.cp[:,start_index:] + except: + pass + + + if(self.end_time != None): + end_index = int(np.argmax(self.time > self.end_time)) + self.time = self.time[:end_index] + # self.cp = self.cp[:,:end_index] + try: + self.p = self.p[:,:end_index] + self.cp = self.cp[:,:end_index] + except: + pass + + +if __name__ == '__main__': + """" + Entry point to read the simulation results from OpenFOAM case and post-process it. + """ + + #CLI parser + parser = argparse.ArgumentParser(description="Get EVENT file from OpenFOAM output") + parser.add_argument('-c', '--case', help="OpenFOAM case directory", required=True) + + arguments, unknowns = parser.parse_known_args() + + + case_path = arguments.case + # case_path = "C:\\Users\\fanta\\Documents\\WE-UQ\\LocalWorkDir\\IsolatedBuildingCFD" + + print("Case full path: ", case_path) + + + #Read JSON data + json_path = os.path.join(case_path, "constant", "simCenter", "input", "IsolatedBuildingCFD.json") + with open(json_path) as json_file: + json_data = json.load(json_file) + + # Returns JSON object as a dictionary + rm_data = json_data["resultMonitoring"] + wc_data = json_data["windCharacteristics"] + duration = json_data["numericalSetup"]["duration"] + + load_output_path = os.path.join(case_path, "constant", "simCenter", "output", "windLoads") + + #Check if it exists and remove files + if os.path.exists(load_output_path): + shutil.rmtree(load_output_path) + + #Create new path + Path(load_output_path).mkdir(parents=True, exist_ok=True) + + + #Read and write the story forces + force_file_name = os.path.join(case_path, "postProcessing", "storyForces", "0", "forces_bins.dat") + origin, time, forces, moments = read_bin_forces(force_file_name) + start_time = 1000 + time = time[start_time:] + forces = forces[start_time:,:,:] + moments = moments[start_time:,:,:] + + + print(force_file_name) + + num_times = len(time) + num_stories = rm_data["numStories"] + + + storyLoads = np.zeros((num_times, num_stories*3 + 1)) + + + storyLoads[:, 0] = time + + + for i in range(num_stories): + storyLoads[:, 3*i + 1] = forces[:,i,0] + storyLoads[:, 3*i + 2] = forces[:,i,1] + storyLoads[:, 3*i + 3] = moments[:,i,2] + + np.savetxt(os.path.join(load_output_path, "storyLoad.txt"), storyLoads, delimiter='\t') + + + #Write base loads + if rm_data["monitorBaseLoad"]: + force_file_name = os.path.join(case_path, "postProcessing", "baseForces", "0", "forces.dat") + origin, time, forces, moments = read_forces(force_file_name) + + time = time[start_time:] + forces = forces[start_time:,:] + moments = moments[start_time:,:] + + print(force_file_name) + + num_times = len(time) + + baseLoads = np.zeros((num_times, 3*2 + 1)) + + baseLoads[:, 0] = time + baseLoads[:, 1:4] = forces + baseLoads[:, 4:7] = moments + + np.savetxt(os.path.join(load_output_path, "baseLoad.txt"), baseLoads, delimiter='\t') + + + #Write base loads + if rm_data["monitorSurfacePressure"]: + p_file_name = "" + if rm_data["importPressureSamplingPoints"]: + p_file_name = os.path.join(case_path, "postProcessing", "importedPressureSamplingPoints") + else: + p_file_name = os.path.join(case_path, "postProcessing", "generatedPressureSamplingPoints") + + + wind_speed = wc_data["referenceWindSpeed"] + print(p_file_name) + + cfd_p = PressureData(p_file_name, u_ref=wind_speed, rho=1.25, p_ref=0.0, start_time=duration*0.1, end_time=None) + + num_times = len(cfd_p.time) + + pressureData = np.zeros((num_times, cfd_p.probe_count + 1)) + + pressureData[:, 0] = cfd_p.time + pressureData[:, 1:] = np.transpose(cfd_p.cp) + + np.savetxt(os.path.join(load_output_path, "pressureData.txt"), pressureData, delimiter='\t') diff --git a/modules/createEVENT/IsolatedBuildingCFD/setup_case.py b/modules/createEVENT/IsolatedBuildingCFD/setup_case.py index b45bb39aa..9f468a3a9 100644 --- a/modules/createEVENT/IsolatedBuildingCFD/setup_case.py +++ b/modules/createEVENT/IsolatedBuildingCFD/setup_case.py @@ -1246,10 +1246,20 @@ def write_controlDict_file(input_json_path, template_dict_path, case_path): # n monitor_base_load = rm_data['monitorBaseLoad'] monitor_surface_pressure = rm_data['monitorSurfacePressure'] +<<<<<<< HEAD + monitor_vtk_planes = rm_data['monitorVTKPlane'] + vtk_planes = rm_data['vtkPlanes'] + + # Need to change this for + max_delta_t = 10*time_step + + #Write 10 times +======= # Need to change this for max_delta_t = 10 * time_step # Write 10 times +>>>>>>> upstream/master write_frequency = 10.0 write_interval_time = duration / write_frequency write_interval_count = int(write_interval_time / time_step) @@ -1319,6 +1329,24 @@ def write_controlDict_file(input_json_path, template_dict_path, case_path): # n added_part = ' #includeFunc baseForces\n' dict_lines.insert(start_index, added_part) +<<<<<<< HEAD + #Write VTK sampling sampling points + if monitor_vtk_planes: + added_part = "" + for pln in vtk_planes: + added_part += " #includeFunc {}\n".format(pln["name"]) + dict_lines.insert(start_index, added_part) + + + + #Write edited dict to file + write_file_name = case_path + "/system/controlDict" + + if os.path.exists(write_file_name): + os.remove(write_file_name) + + output_file = open(write_file_name, "w+") +======= # Write pressure sampling points if monitor_surface_pressure: added_part = ' #includeFunc generatedPressureSamplingPoints\n' @@ -1332,6 +1360,7 @@ def write_controlDict_file(input_json_path, template_dict_path, case_path): # n os.remove(write_file_name) # noqa: PTH107 output_file = open(write_file_name, 'w+') # noqa: SIM115, PTH123 +>>>>>>> upstream/master for line in dict_lines: output_file.write(line) @@ -1980,7 +2009,109 @@ def write_boundary_data_files(input_json_path, case_path): foam.write_foam_field(wind_profiles[:, 8:17], bd_path + 'L') +<<<<<<< HEAD +def write_vtk_plane_file(input_json_path, template_dict_path, case_path): + + #Read JSON data + with open(input_json_path + "/IsolatedBuildingCFD.json") as json_file: + json_data = json.load(json_file) + + # Returns JSON object as a dictionary + rm_data = json_data["resultMonitoring"] + ns_data = json_data["numericalSetup"] + solver_type = ns_data['solverType'] + time_step = ns_data['timeStep'] + + + vtk_planes = rm_data['vtkPlanes'] + write_interval = rm_data['vtkWriteInterval'] + + if rm_data['monitorVTKPlane'] == False: + return + + if len(vtk_planes)==0: + return + + #Write dict files for wind profiles + for pln in vtk_planes: + #Open the template file (OpenFOAM file) for manipulation + dict_file = open(template_dict_path + "/vtkPlaneTemplate", "r") + + dict_lines = dict_file.readlines() + dict_file.close() + + #Write writeControl + start_index = foam.find_keyword_line(dict_lines, "writeControl") + if solver_type=="pimpleFoam": + dict_lines[start_index] = " writeControl \t{};\n".format("adjustableRunTime") + else: + dict_lines[start_index] = " writeControl \t{};\n".format("timeStep") + + #Write writeInterval + start_index = foam.find_keyword_line(dict_lines, "writeInterval") + if solver_type=="pimpleFoam": + dict_lines[start_index] = " writeInterval \t{:.6f};\n".format(write_interval*time_step) + else: + dict_lines[start_index] = " writeInterval \t{};\n".format(write_interval) + + #Write start and end time for the section + start_time = pln['startTime'] + end_time = pln['endTime'] + start_index = foam.find_keyword_line(dict_lines, "timeStart") + dict_lines[start_index] = " timeStart \t\t{:.6f};\n".format(start_time) + + start_index = foam.find_keyword_line(dict_lines, "timeEnd") + dict_lines[start_index] = " timeEnd \t\t{:.6f};\n".format(end_time) + + #Write name of the profile + name = pln["name"] + start_index = foam.find_keyword_line(dict_lines, "planeName") + dict_lines[start_index] = "{}\n".format(name) + + #Write field type + field_type = pln["field"] + start_index = foam.find_keyword_line(dict_lines, "fields") + + if field_type=="Velocity": + dict_lines[start_index] = " fields \t\t({});\n".format("U") + if field_type=="Pressure": + dict_lines[start_index] = " fields \t\t({});\n".format("p") + + #Write normal and point coordinates + point_x = pln["pointX"] + point_y = pln["pointY"] + point_z = pln["pointZ"] + + normal_axis = pln["normalAxis"] + + start_index = foam.find_keyword_line(dict_lines, "point") + dict_lines[start_index] = "\t point\t\t({:.6f} {:.6f} {:.6f});\n".format(point_x, point_y, point_z) + + start_index = foam.find_keyword_line(dict_lines, "normal") + if normal_axis=="X": + dict_lines[start_index] = "\t normal\t\t({} {} {});\n".format(1, 0, 0) + if normal_axis=="Y": + dict_lines[start_index] = "\t normal\t\t({} {} {});\n".format(0, 1, 0) + if normal_axis=="Z": + dict_lines[start_index] = "\t normal\t\t({} {} {});\n".format(0, 0, 1) + + #Write edited dict to file + write_file_name = case_path + "/system/" + name + + if os.path.exists(write_file_name): + os.remove(write_file_name) + + output_file = open(write_file_name, "w+") + for line in dict_lines: + output_file.write(line) + output_file.close() + + +if __name__ == '__main__': + +======= if __name__ == '__main__': +>>>>>>> upstream/master input_args = sys.argv # Set filenames @@ -2027,12 +2158,18 @@ def write_boundary_data_files(input_json_path, case_path): # Write results to be monitored write_base_forces_file(input_json_path, template_dict_path, case_path) write_story_forces_file(input_json_path, template_dict_path, case_path) +<<<<<<< HEAD + write_generated_pressure_probes_file(input_json_path, template_dict_path, case_path) + write_imported_pressure_probes_file(input_json_path, template_dict_path, case_path) + write_vtk_plane_file(input_json_path, template_dict_path, case_path) +======= write_generated_pressure_probes_file( input_json_path, template_dict_path, case_path ) write_imported_pressure_probes_file( input_json_path, template_dict_path, case_path ) +>>>>>>> upstream/master # Write fvSolution dict write_fvSolution_file(input_json_path, template_dict_path, case_path) @@ -2057,3 +2194,9 @@ def write_boundary_data_files(input_json_path, case_path): # Write TInf files write_boundary_data_files(input_json_path, case_path) +<<<<<<< HEAD + + + +======= +>>>>>>> upstream/master diff --git a/modules/createEVENT/IsolatedBuildingCFD/templateOF10Dicts/CMakeLists.txt b/modules/createEVENT/IsolatedBuildingCFD/templateOF10Dicts/CMakeLists.txt index 4ec43b536..ded7069b3 100644 --- a/modules/createEVENT/IsolatedBuildingCFD/templateOF10Dicts/CMakeLists.txt +++ b/modules/createEVENT/IsolatedBuildingCFD/templateOF10Dicts/CMakeLists.txt @@ -24,3 +24,4 @@ simcenter_add_file(NAME fvSolutionTemplate) simcenter_add_file(NAME probeTemplate) simcenter_add_file(NAME storyForcesTemplate) simcenter_add_file(NAME baseForcesTemplate) +simcenter_add_file(NAME vtkPlaneTemplate) diff --git a/modules/createEVENT/IsolatedBuildingCFD/templateOF10Dicts/vtkPlaneTemplate b/modules/createEVENT/IsolatedBuildingCFD/templateOF10Dicts/vtkPlaneTemplate new file mode 100644 index 000000000..367f06859 --- /dev/null +++ b/modules/createEVENT/IsolatedBuildingCFD/templateOF10Dicts/vtkPlaneTemplate @@ -0,0 +1,37 @@ +/*--------------------------------*- C++ -*----------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Version: 10 + \\/ M anipulation | +\*---------------------------------------------------------------------------*/ + +planeName +{ + type surfaces; + libs ("libsampling.so"); + + writeControl timeStep; + writeInterval 10; + timeStart 1; + timeEnd 3; + + surfaceFormat vtk; + fields (U); + + interpolationScheme cellPoint; + + surfaces + ( + sectionPlane + { + type cuttingPlane; + planeType pointAndNormal; + point (0 0 0); + normal (0 1 0); + interpolate true; + } + ); +} + +// ************************************************************************* //