## Domain Setup File import numpy as np import torch import torch.nn as nn from collections import OrderedDict import scipy.io as sio import scipy.stats.qmc as qmc import matplotlib.pyplot as plt import os ''' #### VERSION/ITERATION NAME test_ver = "remake_files" # This will be used when saving tensors for references in other files. Keep this consistent pls :) #### ### Setup working directory - makes it so it's git_main/PINN for ease curr_cwd_path = os.getcwd() working_path = curr_cwd_path+"\PINN" os.chdir(working_path) # print(os.getcwd()) # Data path for where Cyl100 is stored. This will change depending on the user, pls change this if needed pls #data_path = "C://Users//Harry Harris//Documents//GDP Masters Project//Data//Cyl100" #### INPUT PARAMETERS # defining velocity, density, and viscocity for fluid u0 = 1 # lid velocity rho0 = 1 nu0 = 0.01 N_col = 4400 # number of collocation data points [x,y] N_BC = 500 N_circ = 400 r = 0.5 N_data = 1000 ''' def Uniform_with_Cylinder(u0, N_col, N_BC, N_circ, r): """ Create collocation points inside a 1by1 box with a cylinder in the middle Args: u0 : Inlet velocity N_col : number of collocation points # note the acc number is slightly less as we are removing points in the cylinder N_circ : Number of points along cylinder boundary r : Cylinder radius Returns: Tensors of: X_col : Collocation points [x, y] BC_X : Boundary points [x, y] BC_Y : Boundary condition values [u, v] """ # setting up domain size xc, yc = 0, 0 # cylinder center # CREATING collocation points # set random points across the domain using the random normal distribution function hence *3 or *10 to change standard dev and get distribution centred round cylinder y = 2*np.random.randn(N_col,1) x = 4*np.random.randn(N_col,1) xy = np.hstack((x, y)) # Distance from cylinder center dist = np.sqrt((xy[:, 0] - xc)**2 + (xy[:, 1] - yc)**2) # Keep only points outside the cylinder. Boolean mask for cylinder and BCs valid_col = (dist > r) & (xy[:,0]>-10) & (xy[:,1]>-7) & (xy[:,1]< 7) xy_outside = xy[valid_col] X_col = torch.tensor(xy_outside, dtype=torch.float32) # CREATING Boundary condition locations # domain currently set y = -7 to 7 : x = -10 to 25 # evenly distribute BC points along the sides # N_BC = int(np.sqrt(N_col) // 2) y = np.linspace(-7, 7, N_BC) x = np.linspace(-10, 35, N_BC) # turn those points into x,y co-ords b_left = np.column_stack([(np.ones_like(y)*-10), y]) b_right = np.column_stack([np.ones_like(y)*35, y]) b_bottom = np.column_stack([x, np.ones_like(x)*-7]) b_top = np.column_stack([x, np.ones_like(x)*7]) # cylinder BCs theta = np.linspace(0, 2*np.pi, N_circ) cylinder_x = (r * np.cos(theta) + xc) cylinder_y = (r * np.sin(theta) + yc) cylinder_xy = np.column_stack([cylinder_x, cylinder_y]) # Add all BC points BC_X_np = np.vstack([b_left, cylinder_xy]) BC_X = torch.tensor(BC_X_np, dtype=torch.float32) # CREATING the BC_Y essentially 0s everywhere except inlet # tensor of zeros for all walls except inlet # walls_Y = torch.zeros_like(torch.tensor(np.vstack([b_top, b_bottom, b_right]), dtype=torch.float32)) # inlet velocity inlet = torch.zeros_like(torch.tensor(b_left, dtype=torch.float32)) inlet[:, 0] = u0 # Cylinder no slip so 0 cylinder_uv = torch.zeros_like(torch.tensor(cylinder_xy, dtype=torch.float32)) # Combine boundary velocities BC_Y = torch.cat([inlet, cylinder_uv], dim=0) return X_col, BC_X, BC_Y def Data(N, path): engine = qmc.LatinHypercube(d=1) """ loads data from benchmark file Args: Returns: As Tensors (dim = Nx2) X_D : coordinates of data points [x,y] for u,v, dim = [N*2] U_D : training output [u*, v*] at data points X_D, dim = [N*2] PXY_D : coordinates of data points [x,y] for p, dim = [N*2] P_D : training output [p*] at data points PXY_D Notes about data: - flow around a cylinder - Re = 100 , all data has been non-dimensionalised?? , inlet velocity u = 1 """ ## Read file holding training data xy_data = sio.loadmat(path + "/xstar.mat") uv_data = sio.loadmat(path + "/ustar_Aug24.mat") # pressure is calculated at different data points p_data = sio.loadmat(path + "/pstar.mat") pxy_data = sio.loadmat(path + "/xpstar.mat") uv_data_np = np.array(uv_data["ustar"][:,:,0]) u_np = uv_data_np[:,0] v_np = uv_data_np[:,1] xy_data_np = np.array(xy_data["xstar"]) x_np = xy_data_np[:,0] y_np = xy_data_np[:,1] p_data_np = np.array(p_data["pstar"]) pxy_data_np = np.array(pxy_data["xpstar"]) # print(pxy_data_np[:,0]) # Generate indicies of random points in dataset idx = engine.integers(l_bounds=0, u_bounds=len(x_np), n=N) p_idx = engine.integers(l_bounds=0, u_bounds=len(pxy_data_np[:,0]), n=N) x = torch.tensor(x_np[idx], dtype=torch.float32) y = torch.tensor(y_np[idx], dtype=torch.float32) u = torch.tensor(u_np[idx], dtype=torch.float32) v = torch.tensor(v_np[idx], dtype=torch.float32) X_D = torch.column_stack((x, y)) U_D = torch.column_stack((u,v)) PXY_D = torch.tensor(pxy_data_np[p_idx], dtype=torch.float32) P_D = torch.tensor(p_data_np[:,0:1][p_idx], dtype=torch.float32) #print(X_D.size(), U_D.size(),PXY_D.shape(), P_D.shape()) #print(PXY_D.size(), P_D.size()) return X_D, U_D, PXY_D, P_D ''' ## Create Tensors, display domain with different data points, BCs, etc X_col, BC_X, BC_Y= Uniform_with_Cylinder(u0, N_col, N_BC, N_circ, r) X_D, U_D, PXY_D, P_D = Data(N_data, data_path) # print(PXY_D.size(), P_D.size()) ''' ## Plot figure, save it to folder PINN_results def Domain_Plot(X_col, BC_X, X_D, main_path, version_name): plt.figure() plt.title(f"Domain plot insert version type") plt.scatter(X_col[:,0].numpy(),X_col[:,1].numpy(),marker='.', label="Collocation Points") plt.scatter(BC_X[:,0].numpy(),BC_X[:,1].numpy(),marker='x', label="Boundary Condition Points") plt.scatter(X_D[:,0].numpy(), X_D[:,1].numpy(),marker='.', color='r', label="Data Points") plt.xlabel("x [m]") plt.ylabel("y [m]") plt.legend() plt.savefig(f"{main_path}/{version_name}/domain.jpg") print(f"your domain image has been saved to: {main_path}/{version_name}/domain.jpg") plt.close() #plt.show() ''' ## Saving these tensors is just done to see if this works, disregard them when it comes to reproducing PINN, the models are already saved dict = {"X_col" : X_col, "BC_X": BC_X, "BC_Y": BC_Y, "X_D": X_D, "U_D": U_D, "PXY_D": PXY_D, "P_D": P_D} torch.save(dict, f"PINN_results/{test_ver}.pt") ## load tensors just to test tens_dict = torch.load(f"PINN_results/{test_ver}.pt") test_bc_y = tens_dict["BC_Y"] print(test_bc_y.size(), BC_Y.size()) '''