adds new rom preprocessing for additional split runner cfd data (6a9a18f9) · Commits · dt-hydro / dt-hydro-cfd

dthydro_cfd/rom_cfd/dthydro_dash_app.py

0 → 100644

+193 −0

Original line number	Diff line number	Diff line
		"""
		Visualization for the propor Orthogonal decomposition (POD)
		reduced order model (ROM) of the 3D cfd data through the
		turbine chamber.
		This Visualization script depends on the plotly and dash
		python packages for fast live-updating 3D plots.

		NOTE: You must FIRST run the pod_cfd_model.py script
		before executing this script.
		The data produced by
		the pod_cfd_model.py script is visualized by this script
		seperately. This is to decouple the model from the plotting
		procedure so that they may run on seperate processes.

		Author: William Gurecky. wll@ornl.gov. June 2023
		"""
		import numpy as np
		import h5py
		import time
		import argparse
		# plotly imports
		import plotly.graph_objects as go
		import plotly.express as px
		# dash imports
		import dash
		from dash import dcc, html
		from dash.dependencies import Input, Output
		try:
		from dthydro_cfd.rom_cfd.pod_rom import PODInterpModel
		except ImportError:
		from pod_rom import PODInterpModel


		# Connect to h5 file
		parser = argparse.ArgumentParser()
		parser.add_argument("-i", help="Input h5 file", type=str, default="db_dynamic_inelastic.h5")
		args = parser.parse_args()
		db_pod_file = args.i
		f_db = h5py.File(db_pod_file, 'r', libver='latest', swmr=True)

		app = dash.Dash(__name__)
		# html.Div(id='live-vane-angle-text'),
		app.layout = html.Div(
		html.Div([
		html.H4('Reduced Order Model (ROM) Live Feed'),
		dcc.Graph(id='live-update-graph-p'),
		dcc.Graph(id='live-update-graph-v'),
		dcc.Graph(id='live-update-graph-vane'),
		dcc.Graph(id='live-update-graph-flow'),
		dcc.Interval(
		id='interval-component',
		interval=500, # in milliseconds
		n_intervals=0
		)
		])
		)

		X_POINTS = f_db["X"][:].flatten()
		Y_POINTS = f_db["Y"][:].flatten()
		Z_POINTS = f_db["Z"][:].flatten()
		# data thinning
		PTX = np.random.randint(0, len(X_POINTS), size=6000)
		P_ = None
		from collections import deque
		vane_history = deque(maxlen=100)
		flow_history = deque(maxlen=100)
		time_history = deque(maxlen=100)


		# Multiple components can update everytime interval gets fired.
		@app.callback(Output('live-update-graph-p', 'figure'),
		Input('interval-component', 'n_intervals'))
		def update_graph_live(n):
		global P_
		# Read data from database
		P_ = f_db["P"][:]
		# scale pressure to be consisent with model turbine
		# inlet pressure
		p_inlet = 500.e3 # Pa
		p_scale = p_inlet - np.max(P_)

		fig = px.scatter_3d(
		x=X_POINTS[PTX],
		y=Y_POINTS[PTX],
		z=Z_POINTS[PTX],
		color=P_[PTX] + p_scale,
		range_color=[101.e3, 400.e3],
		opacity=0.95,
		size_max=4,
		title="Pressure (Pa)",
		)
		fig.update_traces(marker=dict(size=3))
		fig.update_layout(
		coloraxis_colorbar=dict(title="Pressure (Pa)",),
		)

		# Draw new figure
		# fig = go.Figure(data=go.Scatter3d(
		# x=X_POINTS[PTX],
		# y=Y_POINTS[PTX],
		# z=Z_POINTS[PTX],
		# mode='markers',
		# marker=dict(size=3, color=P_[PTX], opacity=0.95,
		# colorbar=dict(thickness=5, title="Pressure (Pa)", ),
		# ),
		# customdata=P_[PTX],
		# hovertemplate=('%{customdata}%'),
		# ))
		# preserve zoom, pan settings
		# from: https://stackoverflow.com/questions/68798315/how-to-update-plotly-plot-and-keep-ui-settings
		fig['layout']['height'] = 700 # px
		fig['layout']['uirevision'] = 'const'
		return fig


		# Multiple components can update everytime interval gets fired.
		@app.callback(Output('live-update-graph-v', 'figure'),
		Input('interval-component', 'n_intervals'))
		def update_graph_live_v(n):
		# Read data from database
		pred_vx = f_db["VX"][:][PTX]
		pred_vy = f_db["VY"][:][PTX]
		pred_vz = f_db["VZ"][:][PTX]
		v_mag = np.sqrt(pred_vx 2. + pred_vy 2. + pred_vz ** 2.)

		fig_v = go.Figure(data=[go.Cone(
		x=X_POINTS[PTX],
		y=Y_POINTS[PTX],
		z=Z_POINTS[PTX],
		u=pred_vx,
		v=pred_vy,
		w=pred_vz,
		# sizemode="absolute",
		sizeref=3,
		# colorbar=dict(thickness=8, title="Velocity (m/s)"),
		),
		go.Scatter3d(
		x=X_POINTS[PTX],
		y=Y_POINTS[PTX],
		z=Z_POINTS[PTX],
		mode='markers',
		marker=dict(size=2, color=v_mag, opacity=0.65,
		colorbar=dict(thickness=5, title="Velocity (m/s)", ),
		),
		)
		]
		)
		fig_v['layout']['height'] = 700 # px
		fig_v['layout']['uirevision'] = 'const'
		return fig_v


		# Multiple components can update everytime interval gets fired.
		@app.callback(Output('live-update-graph-vane', 'figure'),
		Input('interval-component', 'n_intervals'))
		def update_graph_live_vane(n):
		global vane_history
		global time_history
		t = f_db["t"][:]
		va = f_db["vane_angle"][:]
		vane_history.append(va[0])
		time_history.append(t[0])

		fig_vane = go.Figure(data=go.Scatter(
		x=np.asarray(list(time_history)),
		y=np.asarray(list(vane_history)),
		))
		fig_vane['layout']['uirevision'] = 'const'
		return fig_vane

		# Multiple components can update everytime interval gets fired.
		@app.callback(Output('live-update-graph-flow', 'figure'),
		Input('interval-component', 'n_intervals'))
		def update_graph_live_vane(n):
		global flow_history
		global time_history
		t = f_db["t"][:]
		flow = f_db["flow"][:]
		flow_history.append(flow[0])
		print(flow[0])
		time_history.append(t[0])

		fig = go.Figure(data=go.Scatter(
		x=np.asarray(list(time_history)),
		y=np.asarray(list(flow_history)),
		))
		fig['layout']['uirevision'] = 'const'
		return fig


		if __name__ == "__main__":
		# init the dash app
		app.run_server(debug=True)

dthydro_cfd/rom_cfd/pod_cfd_model.py

+21 −12

Original line number	Diff line number	Diff line
		@@ -23,8 +23,17 @@ try:
		except ImportError:
		from pod_rom import PODInterpModel

		# globals
		v_i_name = "Velocity[i] (m/s)"
		v_j_name = "Velocity[j] (m/s)"
		v_k_name = "Velocity[k] (m/s)"
		p_name = "Static Pressure (Pa)"
		v_i_name = "Velocity in Rotating[i] (m/s)"
		v_j_name = "Velocity in Rotating[j] (m/s)"
		v_k_name = "Velocity in Rotating[k] (m/s)"

		def construct_pod_model(snapshot_h5file, data_name="Pressure (Pa)", n_pod_modes=5):

		def construct_pod_model(snapshot_h5file, data_name="Static Pressure (Pa)", n_pod_modes=5):
		# read the snapshot data from the hdf5 file
		h5f = h5py.File(snapshot_h5file, mode='r')
		xyz_points = h5f["/points"][:]
		@@ -58,7 +67,7 @@ class PodCfdModel(object):
		Creates a Reduced Order Model (ROM) of the 3D CFD data
		of the turbine chamber.
		"""
		def __init__(self, snapshot_h5file: str, data_name: str="Pressure (Pa)", n_pod_modes: int=7):
		def __init__(self, snapshot_h5file: str, data_name: str="Static Pressure (Pa)", n_pod_modes: int=7):
		pod_model, xyz_points = construct_pod_model(snapshot_h5file, data_name, n_pod_modes)
		self._pod_model = pod_model
		self._xyz_points = xyz_points
		@@ -87,16 +96,16 @@ def test_pod_cfd_model(use_plotly=True,
		n_pod_modes=7):
		# create the reduced order model(s)
		pod_model = PodCfdModel(
		snapshot_h5file, data_name="Pressure (Pa)",
		snapshot_h5file, data_name=p_name,
		n_pod_modes=n_pod_modes)
		xyz_points = pod_model.xyz_points

		pod_model_vx = PodCfdModel(
		snapshot_h5file, data_name="Velocity[i] (m/s)")
		snapshot_h5file, data_name=v_i_name)
		pod_model_vy = PodCfdModel(
		snapshot_h5file, data_name="Velocity[j] (m/s)")
		snapshot_h5file, data_name=v_j_name)
		pod_model_vz = PodCfdModel(
		snapshot_h5file, data_name="Velocity[k] (m/s)")
		snapshot_h5file, data_name=v_k_name)

		# evaluate the pod model at an vane angle not in training set
		ts = time.time()
		@@ -124,7 +133,7 @@ def test_pod_cfd_model(use_plotly=True,
		z=Z.flatten(),
		mode='markers',
		marker=dict(size=4, color=values, opacity=0.95,
		colorbar=dict(thickness=5, title="Pressure (Pa)")),
		colorbar=dict(thickness=5, title="Static Pressure (Pa)")),
		customdata=values,
		hovertemplate=('%{customdata}%'),
		))
		@@ -138,7 +147,7 @@ def test_pod_cfd_model(use_plotly=True,
		v=pred_vy,
		w=pred_vz,
		sizemode="absolute",
		sizeref=150,
		# sizeref=150,
		colorbar=dict(thickness=8, title="Velocity (m/s)"),
		))
		fig_v.show()
		@@ -163,16 +172,16 @@ def test_pod_cfd_model(use_plotly=True,
		def test_dynamic_pod_model(snapshot_h5file="cfd_pod_snapshots.h5", n_pod_modes=7):
		# create the reduced order model(s)
		pod_model = PodCfdModel(
		snapshot_h5file, data_name="Pressure (Pa)",
		snapshot_h5file, data_name=p_name,
		n_pod_modes=n_pod_modes)
		xyz_points = pod_model.xyz_points

		pod_model_vx = PodCfdModel(
		snapshot_h5file, data_name="Velocity[i] (m/s)")
		snapshot_h5file, data_name=v_i_name)
		pod_model_vy = PodCfdModel(
		snapshot_h5file, data_name="Velocity[j] (m/s)")
		snapshot_h5file, data_name=v_j_name)
		pod_model_vz = PodCfdModel(
		snapshot_h5file, data_name="Velocity[k] (m/s)")
		snapshot_h5file, data_name=v_k_name)

		# update vane angle and store data in tmp hdf5 file
		# acting as a database

dthydro_cfd/rom_cfd/preprocess_cfd.py

+159 −25

Original line number	Diff line number	Diff line
		import numpy as np
		import pandas as pd
		import argparse
		import h5py
		import os
		from scipy.interpolate import griddata
		from six import iteritems


		#run_csv_files = [
		# "fixed_3pt95.csv",
		# "fixed_4pt48.csv",
		# "fixed_5pt01.csv",
		# "fixed_7pt50.csv",
		# "fixed_8pt51.csv",
		# "fixed_9pt52.csv",
		# ]
		# run_parameters = [
		# [3.95,],
		# [4.48,],
		# [5.01,],
		# [7.50,],
		# [8.51,],
		# [9.52,],
		# ]
		# NTNU CFD FILES
		run_dir = "./ntnu/"
		run_csv_files = [
		"alpha-3pt95_fixed.csv",
		"alpha-4pt48_fixed.csv",
		@@ -49,7 +36,100 @@ run_parameters = [
		[9.01],
		]

		def conv2hdf(csv_fname, outfile_name="fixed_cfd.h5", case_params=[], dset_names=None, **kwargs):
		# NEW CFD DATA FILES
		run_dir = "./alder_u12/"
		run_csv_files = [
		"alpha-01_runner.csv",
		"alpha-02_runner.csv",
		"alpha-03_runner.csv",
		"alpha-04_runner.csv",
		"alpha-05_runner.csv",
		"alpha-06_runner.csv",
		"alpha-07_runner.csv",
		"alpha-08_runner.csv",
		"alpha-09_runner.csv",
		"alpha-10_runner.csv",
		]
		run_csv_files = [
		"alpha-01_spiral-casing-and-draft-tube.csv",
		"alpha-02_spiral-casing-and-draft-tube.csv",
		"alpha-03_spiral-casing-and-draft-tube.csv",
		"alpha-04_spiral-casing-and-draft-tube.csv",
		"alpha-05_spiral-casing-and-draft-tube.csv",
		"alpha-06_spiral-casing-and-draft-tube.csv",
		"alpha-07_spiral-casing-and-draft-tube.csv",
		"alpha-08_spiral-casing-and-draft-tube.csv",
		"alpha-09_spiral-casing-and-draft-tube.csv",
		"alpha-10_spiral-casing-and-draft-tube.csv",
		]
		run_parameters = [
		[1.3875],
		[2.775],
		[4.1625],
		[5.55],
		[6.9375],
		[8.325],
		[9.7125],
		[11.1],
		[12.4875],
		[13.875],
		]


		def plot_cfd_vel_field(mesh_x, mesh_y, mesh_z, v_i, v_j, v_k):
		import plotly.graph_objects as go
		import plotly

		n_base_points = len(mesh_x)
		n_subsample = 100000
		idxs = np.arange(n_base_points, dtype=int)
		subsample_idxs = np.random.choice(idxs, size=min(n_subsample, n_base_points), replace=False)

		v_mag = np.sqrt(v_i 2 + v_j 2 + v_k ** 2)
		print("min vmag:", np.min(v_mag))
		print("max vmag:", np.max(v_mag))
		fig_v = go.Figure(data=go.Cone(
		x=mesh_x[subsample_idxs].flatten(),
		y=mesh_y[subsample_idxs].flatten(),
		z=mesh_z[subsample_idxs].flatten(),
		u=v_i[subsample_idxs],
		v=v_j[subsample_idxs],
		w=v_k[subsample_idxs],
		sizemode="absolute",
		sizeref=150,
		colorbar=dict(thickness=8, title="CFD Velocity (m/s)"),
		))
		fig_v.show()
		fig_p = go.Figure(data=go.Scatter3d(
		x=mesh_x[subsample_idxs].flatten(),
		y=mesh_y[subsample_idxs].flatten(),
		z=mesh_z[subsample_idxs].flatten(),
		mode='markers',
		marker=dict(size=4, color=v_mag[subsample_idxs], opacity=0.95,
		colorbar=dict(thickness=5, title="CFD Vel Magnitude (m/s)")),
		customdata=v_mag[subsample_idxs],
		hovertemplate=('%{customdata}%'),
		))
		fig_p.show()


		def plot_cfd_pres_field(mesh_x, mesh_y, mesh_z, p_s):
		import plotly.graph_objects as go
		import plotly
		fig_p = go.Figure(data=go.Scatter3d(
		x=mesh_x.flatten(),
		y=mesh_y.flatten(),
		z=mesh_z.flatten(),
		mode='markers',
		marker=dict(size=4, color=p_s, opacity=0.95,
		colorbar=dict(thickness=5, title="CFD Static Pressure (Pa)")),
		customdata=p_s,
		hovertemplate=('%{customdata}%'),
		))
		fig_p.show()


		def conv2hdf(csv_fname, outfile_name="fixed_cfd.h5", case_params=[], dset_names=None, make_plot=False, **kwargs):
		"""
		Read data from csv and convert to compressed hdf5 format
		for faster processing.
		@@ -80,6 +160,18 @@ def conv2hdf(csv_fname, outfile_name="fixed_cfd.h5", case_params=[], dset_names=
		if case_params:
		hdf.create_dataset("/params", data=np.array([case_params]).flatten())

		# Plot raw cfd velocities
		if make_plot:
		xx = hdf["/X (m)"][:]
		yy = hdf["/Y (m)"][:]
		zz = hdf["/Z (m)"][:]
		v_i = hdf["/Velocity[i] (m/s)"][:]
		v_j = hdf["/Velocity[j] (m/s)"][:]
		v_k = hdf["/Velocity[k] (m/s)"][:]
		p_s = hdf["/Static Pressure (Pa)"][:]
		plot_cfd_vel_field(xx, yy, zz, v_i, v_j, v_k)
		plot_cfd_pres_field(xx, yy, zz, p_s)

		# clean up
		hdf.close()

		@@ -110,7 +202,7 @@ def reinterpolate_cfd(
		n_base_points = len(base_x)
		print("Number of base points: %d" % n_base_points)
		idxs = np.arange(n_base_points, dtype=int)
		subsample_idxs = np.random.choice(idxs, size=n_subsample, replace=False)
		subsample_idxs = np.random.choice(idxs, size=min(n_subsample, n_base_points), replace=False)
		sub_base_x = base_x[subsample_idxs]
		sub_base_y = base_y[subsample_idxs]
		sub_base_z = base_z[subsample_idxs]
		@@ -129,6 +221,7 @@ def reinterpolate_cfd(

		hdf_out.create_dataset("/points", data=sub_base_xyz, compression=True)

		data_summaries = []
		for h5f_name in h5_fnames:
		print("Interpolating data from: %s" % h5f_name)
		hdf = h5py.File(h5f_name, mode='r')
		@@ -152,8 +245,22 @@ def reinterpolate_cfd(

		# write selected, interpolated field data to h5
		h5f_name.split(sep='.')
		print("Max field val: %0.5e" % np.max(case_sub_field_data))
		print("Min field val: %0.5e" % np.min(case_sub_field_data))

		# Summary data
		if "Velocity" in field_name:
		summary_field_data = np.abs(case_sub_field_data)
		else:
		summary_field_data = case_sub_field_data
		data_summaries.append([
		h5f_name + ", " + str(field_name),
		np.max(summary_field_data),
		np.min(summary_field_data),
		np.mean(summary_field_data),
		np.var(summary_field_data),
		np.quantile(summary_field_data, [0.05, 0.3, 0.5, 0.7, 0.95]),
		])

		# Write the re-interpolated data to hdf5
		hdf_out.create_dataset(
		"/snapshots/" + h5f_name + "/" + field_name, data=case_sub_field_data,
		compression=True)
		@@ -168,19 +275,46 @@ def reinterpolate_cfd(
		# clean up
		hdf.close()
		hdf_out.close()
		print("Fname, field, Max, Min, Mean, Var")
		for summary in data_summaries:
		print("%s, %0.3e, %0.3e, %0.3e, %0.3e" % (summary[0], summary[1], summary[2], summary[3], summary[4]))
		print("Fname, field, Quantiles")
		for summary in data_summaries:
		print("%s, %0.3e, %0.3e, %0.3e, %0.3e, %0.3e" % \
		(summary[0], summary[-1][0], summary[-1][1], summary[-1][2], summary[-1][3], summary[-1][4]))


		if __name__ == "__main__":
		outfile_name = ""
		parser = argparse.ArgumentParser()
		parser.add_argument("-o", help="output h5 file name", type=str, default='cfd_pod_snapshots.h5')
		args = parser.parse_args()
		h5_fnames = []
		field_names = ["Pressure (Pa)","Velocity[i] (m/s)",
		"Velocity[j] (m/s)", "Velocity[k] (m/s)",]
		# field_names = [
		# "Static Pressure (Pa)",
		# "Velocity[i] (m/s)",
		# "Velocity[j] (m/s)",
		# "Velocity[k] (m/s)",
		# ]
		field_names = [
		"Static Pressure (Pa)",
		"Velocity in Rotating[i] (m/s)",
		"Velocity in Rotating[j] (m/s)",
		"Velocity in Rotating[k] (m/s)",
		"Velocity[i] (m/s)",
		"Velocity[j] (m/s)",
		"Velocity[k] (m/s)",
		]
		loc_names = ["X (m)","Y (m)","Z (m)"]
		csv_col_names = field_names + loc_names
		i = 0
		for csv_file, params in zip(run_csv_files, run_parameters):
		print("Preprocessing file: %s" % csv_file)
		out_file = csv_file + ".h5"
		conv2hdf(csv_file, outfile_name=out_file, case_params=params,
		make_plot = False
		if i == 7:
		make_plot = True
		conv2hdf(run_dir + csv_file, outfile_name=out_file, case_params=params, make_plot=make_plot,
		csv_col_names=csv_col_names)
		h5_fnames.append(out_file)
		reinterpolate_cfd(h5_fnames[0], h5_fnames=h5_fnames, field_names=field_names)
		i += 1
		reinterpolate_cfd(h5_fnames[0], h5_fnames=h5_fnames, field_names=field_names, out_fname=args.o)

examples/alder_inelastic_penstock.py

+77 −5

Original line number	Diff line number	Diff line
		@@ -400,6 +400,65 @@ def alder_inelastic(t_end=20., h_in=60., h_out=0.0, pid_control=True, args, *k
		plt.close()


		def write_to_h5_db(f_db, res_dict):
		t = res_dict["time"]
		current_vane_angle = np.rad2deg(res_dict["turbine_vane_angle"])
		pred_p = res_dict['3d_pred_p']
		pred_vx = res_dict['3d_pred_vx']
		pred_vy = res_dict['3d_pred_vy']
		pred_vz = res_dict['3d_pred_vz']
		pred_flow = res_dict['flow_rate']
		p_s = res_dict['penstock_pressures']
		s = res_dict['penstock_pressure_points']
		turb_power = res_dict['turbine_shaft_power']
		turb_w = res_dict['turbine_rotation_rate']

		# 3D points for turbine chamber plots
		xyz_points = res_dict['3d_xyz_points']
		X = xyz_points[:, 0]
		Y = xyz_points[:, 1]
		Z = xyz_points[:, 2]

		try:
		# overwrite datasets
		f_db["P"][...] = pred_p
		f_db["VX"][...] = pred_vx
		f_db["VY"][...] = pred_vy
		f_db["VZ"][...] = pred_vz
		# 1D system vars
		f_db["vane_angle"][...] = [current_vane_angle]
		f_db["t"][...] = [float(t)]
		f_db["P_s"][...] = p_s
		f_db["s"][...] = s
		f_db["flow"][...] = [pred_flow]
		f_db["turb_w"][...] = [turb_w]
		f_db["turb_power"][...] = [turb_power]
		except:
		# create new datasets
		f_db.create_dataset("X", data=X)
		f_db.create_dataset("Y", data=Y)
		f_db.create_dataset("Z", data=Z)
		# turbine chamber 3d
		f_db.create_dataset("P", data=pred_p)
		f_db.create_dataset("VX", data=pred_vx)
		f_db.create_dataset("VY", data=pred_vy)
		f_db.create_dataset("VZ", data=pred_vz)
		# 1D system vars
		f_db.create_dataset("vane_angle", data=[current_vane_angle])
		f_db.create_dataset("t", data=[float(t)])
		# 1D pressure along penstock
		f_db.create_dataset("P_s", data=p_s)
		# 1D locations along penstock (x-axis var)
		f_db.create_dataset("s", data=s)
		# flow rate
		f_db.create_dataset("flow", data=[pred_flow])
		# turbine rotation rate
		f_db.create_dataset("turb_w", data=[turb_w])
		# turbine power
		f_db.create_dataset("turb_power", data=[turb_power])
		f_db.flush()


		if __name__ == "__main__":
		import argparse
		parser = argparse.ArgumentParser()
		@@ -408,11 +467,24 @@ if __name__ == "__main__":
		parser.add_argument("-h_out", type=float, help="tailwater depth", default=0.0) # 6.87
		parser.add_argument("--pod_rom_snapshots", type=str, help="file name of CFD POD snapshots", default="")
		args = parser.parse_args()
		for res_dict in alder_inelastic(

		# Open live hdf5 file for logging
		db_pod_file = "db_dynamic_inelastic.h5"
		if os.path.exists(db_pod_file):
		os.remove(db_pod_file)
		f_db = h5py.File(db_pod_file, 'w', libver='latest')
		f_db.swmr_mode = True
		f_db.flush()

		for i, res_dict in enumerate(alder_inelastic(
		t_end=args.t_f,
		h_in=args.h_in,
		h_out=args.h_out,
		pod_rom_snapshots=args.pod_rom_snapshots):
		# print(res_dict)
		# Could do realtime additional plotting here
		pass
		pod_rom_snapshots=args.pod_rom_snapshots)):
		# write results to database every 100 steps
		if i % 100 == 0:
		write_to_h5_db(f_db, res_dict)
		time.sleep(0.25)

		# clean up
		f_db.close()

Admin message