import numpy as np import torch import torch.nn as nn import matplotlib.pyplot as plt from collections import OrderedDict import scipy.io as sio import scipy.stats.qmc as qmc import pandas as pd import os from matplotlib.colors import Normalize from matplotlib.gridspec import GridSpec from matplotlib.patches import Circle ''' #### Get data, use consistent variables ## Keep these variables consistent between files!!!!!!!! u0 = 1 rho0 = 1 nu0 = 0.01 r=0.5 test_ver = "remake_files" N_domain = 4000 ### 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) ### Code to remake PINN model ##### Define Sub-model ## CREATING THE MINI MODEL FOR u,v AND p AND THEN COMBINE THEM INTO ONE MODEL class submodel(nn.Module): def __init__( self, N_input, N_hidden_arr, N_output, activation = nn.Tanh ): super(submodel, self).__init__() # Create network # Create input layer w/ activation function layers = [('Input', nn.Linear(N_input, N_hidden_arr[0]))] layers.append(('Input activation', activation())) # Create hidden layers for i in range(len(N_hidden_arr)-1): layers.append( ("Hidden %d" % (i+1), nn.Linear(N_hidden_arr[i], N_hidden_arr[i+1])) ) layers.append(('Hidden activation %d' % (i+1), activation())) layers.append(('Output', nn.Linear(N_hidden_arr[-1], N_output))) layerdict = OrderedDict(layers) self.layers = nn.Sequential(layerdict) def forward(self, x): y = self.layers(x) return y ### Define class for one model, containing several sub-models class Net(nn.Module): def __init__(self, N_input, N_hidden_arr, N_output, activation = nn.Tanh): super(Net, self).__init__() # Create network # creates three models using submodel as the blueprint self.model_u = submodel(N_input, N_hidden_arr ,N_output, activation) self.model_v = submodel(N_input, N_hidden_arr, N_output, activation) self.model_p = submodel(N_input, N_hidden_arr, N_output, activation) # combine the outputs of all the models into a single output def forward(self, xy): out_u = self.model_u(xy) out_v = self.model_v(xy) out_p = self.model_p(xy) combined = torch.cat((out_u, out_v, out_p), dim=1) return combined ### This model is stripped back, only needed for a forward pass to be made class PINN: def __init__(self): ### Need to change this to check if gpu available as well device = torch.device("cpu") print("Using CPU") self.model = Net( N_input=2, N_hidden_arr=[32,16,16,32], N_output = 1 ).to(device) Tensors = torch.load(f"PINN_results/{test_ver}.pt") # DATA CREATION self.X, self.BC_X, self.BC_Y = Tensors["X_col"], Tensors["BC_X"], Tensors["BC_Y"] self.X_D, self.U_D, self.PXY_D, self.P_D = Tensors["X_D"], Tensors["U_D"], Tensors["PXY_D"], Tensors["P_D"] # copy and seperate tensors into format for loss calculations self.X = self.X.clone().detach().requires_grad_(True) self.BC_X = self.BC_X.clone().detach().requires_grad_(True) self.X_D = self.X_D.clone().detach().requires_grad_(True) self.PXY_D = self.PXY_D.clone().detach().requires_grad_(True) # OPTIMISERS self.optimiser = torch.optim.LBFGS( params=self.model.parameters(), lr=1.0, max_iter = 25*10**3, max_eval = 25*10**3, history_size=50, tolerance_change=1e-7, tolerance_grad=1e-7, line_search_fn="strong_wolfe" ) self.adam = torch.optim.Adam(self.model.parameters()) self.loss_fn = nn.MSELoss() # Counter for printing loss self.iter =1 def eval(self): self.model.eval() ##### def Domain_Generation(N_col=N_domain, r=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 = torch.tensor(xy_outside, dtype=torch.float32) return X X = Domain_Generation() X = X.clone().detach().requires_grad_(True) ## Try saving model state?? def pressure_distribution(r, x0, y0, model): d = 0.01 theta_u = np.arange(0, np.pi+d, d) theta_l = np.arange(np.pi, 2*np.pi+d, d) x_u = (r * np.cos(theta_u) + x0) y_u = (r * np.sin(theta_u) + y0) x_l = (r * np.cos(theta_l) + x0) y_l = (r * np.sin(theta_l) + y0) xy_u = np.column_stack([x_u, y_u]) xy_l = np.column_stack([x_l, y_l]) X_u = torch.Tensor(xy_u) X_l = torch.Tensor(xy_l) preds_u = model(X_u) preds_l = model(X_l) p_u = preds_u[:, 2:3] #pressure p_np_u = p_u.detach().numpy() p_l = preds_l[:, 2:3] p_np_l = p_l.detach().numpy() cp_u = p_np_u/(0.5*rho0*u0**2) cp_l = p_np_l/(0.5*rho0*u0**2) return cp_u, cp_l, xy_u, xy_l, theta_u, theta_l # Test: Create new PINN instance and load the model new = PINN() # new.model.load_state_dict(torch.load(f"PINN_results/Model ver - {test_ver}.pt")) # This is the test model referenced in PINN_results/ directory new.model.load_state_dict(torch.load("model_data_lim_new_BC.txt")) ## Here is the model trained before using Cylinder_flow__3.ipynb, losses are not tracked for this one # forward pass of the model !!! fwd pass outside eval/inference is okay as long as you dont .backward() # we have to call .model() because our model is nested inside class PINN # y_preds = test.model(X) y_preds = new.model(X) ''' # Trying to save model state to a text file --- it works! def Plotting(model,Tensors, results_path, version_name): y_preds = model.model(Tensors["X_col"]) u = y_preds[:, 0:1] # u velocity v = y_preds[:, 1:2] p = y_preds[:, 2:3] #pressure # convert to array for postprocessing p_np = p.detach().numpy() u_np = u.detach().numpy() v_np = v.detach().numpy() x_np = Tensors["X_col"][:, 0:1].detach().numpy() y_np = Tensors["X_col"][:, 1:2].detach().numpy() # create data points to plot a circle on the plots xc,yc, r = 0,0,0.5 theta = np.linspace(0, 2*np.pi, 400) cylinder_x = (r * np.cos(theta) + xc) cylinder_y = (r * np.sin(theta) + yc) def tricontour(gs, x, y, z, title): """ this is an explanation for tripcontour() and tripplot() which are two different methods for plotting contour type plots Args: grid: plot position in subplot x: x-array (x coords) y: y-array (y coords) z: z-array (engineering value at [x,y]) title: title """ plt.subplot(gs) tcf = plt.tricontourf(x, y, z, levels=50, cmap='rainbow') # Add the circle plt.fill(cylinder_x,cylinder_y, color='blue', alpha=0.75) plt.colorbar(tcf) plt.gca().set_aspect('equal') plt.title(title) plt.xlabel('x') plt.ylabel('y') # note we have to squeeze the arrays as they are in column vector format(ML format) gs = GridSpec(2,2) tricontour(gs[0,0], np.squeeze(x_np), np.squeeze(y_np), np.squeeze(p_np), 'Pressure contour') tricontour(gs[0,1], np.squeeze(x_np), np.squeeze(y_np), np.squeeze(u_np), 'u velocity contour') tricontour(gs[1,0], np.squeeze(x_np), np.squeeze(y_np), np.squeeze(v_np), 'v velocity contour') plt.tight_layout() plt.savefig(f"{results_path}/{version_name}/contour.jpg") print(f"a contour has been saved to: {results_path}/{version_name}/contour.jpg") plt.close() #plt.show() def tripplot(gs, x, y, z, title): plt.subplot(gs) plt.tripcolor(x, y, z, cmap='rainbow') # Add the circle plt.fill(cylinder_x,cylinder_y, color='blue', alpha=0.75) plt.ylim(-2.5,2.5) plt.xlim(-2.5,5) plt.colorbar() plt.title(title) plt.xlabel('x') plt.ylabel('y') gs2 = GridSpec(2,2) tripplot(gs2[0,0], np.squeeze(x_np), np.squeeze(y_np), np.squeeze(p_np), 'Pressure contour') tripplot(gs2[0,1], np.squeeze(x_np), np.squeeze(y_np), np.squeeze(u_np), 'u velocity contour') tripplot(gs2[1,0], np.squeeze(x_np), np.squeeze(y_np), np.squeeze(v_np), 'v velocity contour') plt.tight_layout() plt.savefig(f"{results_path}/{version_name}/zoomed_contour.jpg") print(f"a zoomed in contour has been saved to: {results_path}/{version_name}/zoomed_contour.jpg") plt.close() #plt.show() ''' # Get pressure distribution for cylinder cp_u, cp_l, xy_u, xy_l, th_u, th_l = pressure_distribution(0.5, 0, 0, new.model) # Plot stuff x_u = xy_u[:,0] y_u = xy_u[:,1] x_l, y_l = xy_l[:,0], xy_l[:,1] plt.figure() plt.plot(x_u, cp_u, '-b', label="Upper Surface") plt.plot(x_l, cp_l, '-r', label="Lower Surface") plt.xlabel("x [m]") plt.ylabel("Cp") plt.legend() plt.grid() plt.show() #### get surface cps cp_u_sin = np.array([]) cp_l_sin = np.array([]) cp_u_cos = np.array([]) cp_l_cos = np.array([]) # had to do this weirdly idk why, tried doing cp*np.sin(th) but just got an obscenely large matrix for i in range(0, len(cp_u)): val = cp_u[i]*np.sin(th_u[i]) cp_u_sin = np.append(cp_u_sin, val) print() for i in range(0, len(cp_l)): val = cp_l[i]*np.sin(th_l[i]) cp_l_sin = np.append(cp_l_sin, val) for i in range(0, len(cp_u)): val = cp_u[i]*np.cos(th_u[i]) cp_u_cos = np.append(cp_u_cos, val) for i in range(0, len(cp_l)): val = cp_l[i]*np.cos(th_l[i]) cp_l_cos = np.append(cp_l_cos, val) ### Find CD by taking horizontal component of Cp along upper and lower surface plt.figure() plt.plot(th_u, cp_u_sin, '-b', label="Upper surface vert. component") plt.plot(th_l, cp_l_sin, '-r', label="Lower surface vert. component") plt.xlabel("theta [rad]") plt.ylabel("Cp*sin(theta)") plt.legend() plt.grid() plt.show() plt.figure() plt.plot(th_u, cp_u_cos, '-b', label="Upper surface vert. component") plt.plot(th_l, cp_l_cos, '-r', label="Lower surface vert. component") plt.xlabel("theta [rad]") plt.ylabel("Cp*cos(theta)") plt.legend() plt.grid() plt.show() print(cp_u_cos) CD = -1*(np.trapz(cp_u_cos, th_u)+np.trapz(cp_l_cos, th_l)) CL = np.trapz(cp_u_sin, th_u)+np.trapz(cp_l_sin, th_l) print(f"Lift Coefficient CL: {CL:.5f}") print(f"Drag Coefficient CD: {CD:.5f}") ''' # Need to find a way to get drag coefficient???? # Or at least calculate drag ### Trying formulae given in source: https://community.ptc.com/sejnu66972/attachments/sejnu66972/PTCMathcad/41405/1/AeroLab1.pdf