Machine Learning in Aerodynamics, Propulsion and Fracture Mechanics!
SolidWorks Videos
Wednesday, March 27, 2024
Teaching the Navier-Stokes Equations to Sand: Data-Less Physics Informed Neural Network (Includes PyTorch Code, Validated)
Here is a simple code I wrote for Data-less đ Physics Informed Neural đ§ Network đ¸️ to solve the steady-state 2D incompressible Navier-Stokes equations in Python using PyTorch. Already validated results will be uploaded here in good time. The code has all the comments đ¨đĢ you could possibly need. If you have questions, well... đš
Please note that a properly cooled đ§ GPU is recommended for running this code. If your CPU melts đĨ, don't blame me đ¤. Furthermore, if you classify yourself as a poor peasant đ¨đž, use: device = torch.device("cpu") # cuda for gpu, cpu for cpu đ.
The New Code: (Fall 2025)
After successfully đ completing the validation for the SIMB, yours truly has been developing in the abundant spare time đ; it is now the perfect time to revisit PINNs. On a fine morning, yours truly decided to separate đ¯đđ§♡ the tensors for boundary conditions and the equations of motion. It turns out đą that making separate functions for the boundary conditions and the equations of motion with separate tensors īŽŠŲ¨ŲīŽŠīŽŠŲ¨Ų♡īŽŠŲ¨ŲīŽŠīŽŠŲ¨Ų for each boundary makes the loss function converge to the global minimum đŗ️ much faster ⚡.
Copyright <2025> <Fahad Butt>
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
The helper function to define Navier-Stokes equations is defined first, navier_stokes_equation. This function calculates the the model output at interior points and then calculates the derivatives using autograd to assemble the equations of motion.
The boundary condition helper function is defined next, boundary_conditions. This function evaluates the model output at each boundary separately and then enforces the boundary conditions. As a result of these modifications đ ️, the loss đ is also rewritten.
def navier_stokes_equation(ns):
ns.requires_grad_(True) # watch ns for gradient computation
out = model(ns) # output from model at interior
u = out[:, 0] # u interior
v = out[:, 1] # v interior
p = out[:, 2] # p interior
du = torch.autograd.grad(u, ns, torch.ones_like(u), create_graph=True)[0] # du/dns
dv = torch.autograd.grad(v, ns, torch.ones_like(v), create_graph=True)[0] # dv/dns
nn.MSELoss()(v_bottom_bc, torch.zeros_like(v_bottom_bc)) # v momentum boundary condition loss
bc_loss = u_bc_loss + v_bc_loss # total boundary condition loss
mse_loss = x_mom_loss + y_mom_loss + cont_loss # total equation loss
total_loss = mse_loss + bc_loss # total loss
return total_loss
Instead of creating a single tensor đ§ for the entire input, five 𧎠tensors are created for this new, at least to yours truly, method. Of course, the training loop also has to be modified đŠ.
B_left = torch.tensor(cp.vstack((x_l, y_l)).T, dtype=torch.float32).to(device) # stack x, y together
optimizer.zero_grad() # clear gradients of all optimized variables
loss_value = loss_fn() # compute prediction by passing inputs to model
loss_value.backward() # compute gradient of loss with respect to model parameters
optimizer.step() # perform parameter up date
The Code:
#%% license
#Copyright <2024> <Fahad Butt>
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#%% import necessary libraries
import os
import numpy as np
import random as rn
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.nn.init as init
import torch.optim as optim
# import torch_directml # enable for any gpu except CUDA gpu
import glob
device = torch.device("cpu") # cuda for CUDA gpu, cpu for cpu, torch_directml.device() for any other gpu
#%% 2D navier-stokes equations and boundary conditions
def navier_stokes_equation(ns):
ns.requires_grad_(True) # watch ns for gradient computation
out = model(ns) # output from model
u = out[:, 0] # x velocity
v = out[:, 1] # y velocity
p = out[:, 2] # pressure
du = torch.autograd.grad(u, ns, torch.ones_like(u), create_graph=True)[0] # du/dns
dv = torch.autograd.grad(v, ns, torch.ones_like(v), create_graph=True)[0] # dv/dns
num_trainable_params = count_parameters(model) # show trainable parameters
print("Number of trainable parameters:", num_trainable_params)
#%% prepare input data
train_points = np.vstack((x, y)).T # stack x, y together
train_points = torch.tensor(train_points, dtype=torch.float32).to(device) # convert to tensor and move training points to GPU / CPU, use float64 if mroe precision is required
#%% train using Adam
resume_training = False # True if resuming training
start_epoch = 1
optimizer = optim.Adam(model.parameters(), lr = 0.00001) # select optimizer and parameters
checkpoint_dir = './pytorch_checkpoints' # directory to save checkpoints
os.makedirs(checkpoint_dir, exist_ok=True) # create a directory to save checkpoints
x_m_te = np.random.uniform(0, L, 50000) # mid field
y_m_te = np.random.uniform(0, D, 50000)
x_te = np.concatenate([x_l_te, x_r_te, x_b_te, x_t_te, x_m_te]) # x co-ordinates
y_te = np.concatenate([y_l_te, y_r_te, y_b_te, y_t_te, y_m_te]) # y co-ordinates
test_points = np.vstack((x_te, y_te)).T # stack x, y together
test_points = torch.tensor(test_points, dtype=torch.float32).to(device) # convert to tensor and move training points to GPU / CPU, use float64 if mroe precision is required
predict = model(test_points).detach().cpu().numpy() # predict and move back to CPU
u = predict[:, 0] * Uinf # predict and bring back to dimensional form
sc = plt.scatter(x_te, y_te, c = u, cmap = 'jet', s = 1, edgecolor = 'none', alpha = 1)
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.axis('off')
plt.show()
plt.figure(dpi = 300)
sc = plt.scatter(x_te, y_te, c = v, cmap = 'jet', s = 1, edgecolor = 'none', alpha = 1)
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.axis('off')
plt.show()
plt.figure(dpi = 300)
sc = plt.scatter(x_te, y_te, c = p, cmap = 'jet', s = 1, edgecolor = 'none', alpha = 1)
plt.colorbar(orientation = 'vertical')
plt.gca().set_aspect('equal', adjustable = 'box')
plt.xticks([0, L])
plt.yticks([0, D])
plt.xlabel('x [m]')
plt.ylabel('y [m]')
plt.axis('off')
plt.show()
Flow between Flat - Plates
The boundary conditions are taken from [2]. The results presented in Fig. 1 are from flow between two parallel flat plates. For these boundary conditions, analytical solutions to the Navier-Stokes equations are available. As this is "just another blog" đ; just by eye-balling đ, the results look very promising. I mean from Fig. 1, the velocities at left side are captured well-ish!
Fig. 1, post processing
Lid-Driven Cavity
The case of flow between two flat plates shown in Fig. 1 is self-validating as area times velocity should be same for inlet and outlet, which it is. Lid-Driven Cavity however, is a benchmark problem for many numerical methods in fluid dynamics because of complex vortex dynamics. Here I present, the results of Lid-Driven Cavity problem. The results are very good đ. The results are presented in Fig. 2. The comparison of PINN results with [1] is shown in Fig. 3. The validation of CFD is available to be read here. While the CFD code is available here.
Fig. 2, Post processing
Fig. 3, A comparison
External Fluid Dynamics
Here is a video of the PINN learning flow around obstacle i.e. a square cylinder.
Fig. 4, The animation
The PINN is now capable to predict flow around curved objects. The flow around circular cylinder is mentioned. I use polar coordinates for applying Neuman boundary conditions for pressure to satisfy the no-slip wall. The results are shown in Fig. 5. The code is also made available for fellow researchers to improve upon and use in the research, please cite properly. The results show u and v velocities at the left portion while pressure and streamlines are shown on the right side. As compared to the commercial software that shall not be named đ, the results are very accurate. The Reynolds number is kept at 10to avoid making the loss landscape complex. Please refer to this post for more details!
Fig. 5, The post processing
Code
# Copyright <2024> <FAHAD BUTT>
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
#%% import necessary libraries
import os
import numpy as np
import random as rn
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
import torch.nn.init as init
import torch.optim as optim
# import torch_directml # enable for AMD and Intel gpu
import glob
from scipy.interpolate import griddata
device = torch.device("cpu") # cuda for CUDA gpu, cpu for cpu, mps for Apple GPU, torch_directml.device() for any other gpu
#%% 2D navier-stokes equations and boundary conditions
def navier_stokes_equation(ns):
ns.requires_grad_(True) # watch ns for gradient computation
out = model(ns) # output from model
u = out[:, 0] # x velocity
v = out[:, 1] # y velocity
p = out[:, 2] # pressure
du = torch.autograd.grad(u, ns, torch.ones_like(u), create_graph = True)[0] # du/dns
dv = torch.autograd.grad(v, ns, torch.ones_like(v), create_graph = True)[0] # dv/dns
num_trainable_params = count_parameters(model) # show trainable parameters
print("Number of trainable parameters:", num_trainable_params)
#%% prepare input data
train_points = np.vstack((x_train, y_train)).T # stack x, y together
train_points = torch.tensor(train_points, dtype=torch.float32).to(device) # convert to tensor and move training points to GPU / CPU, use float64 if more precision is required
#%% train using Adam
resume_training = False # True if resuming training
start_epoch = 1
optimizer = optim.Adam(model.parameters(), lr = 0.0001) # select optimizer and parameters
checkpoint_dir = './pytorch_checkpoints' # directory to save checkpoints
os.makedirs(checkpoint_dir, exist_ok = True) # create a directory to save checkpoints
x_te = np.concatenate([x_l_te, x_r_te, x_b_te, x_t_te, x_m_te]) # x co-ordinates
y_te = np.concatenate([y_l_te, y_r_te, y_b_te, y_t_te, y_m_te]) # y co-ordinates
distance_te = np.sqrt((x_te - center_x) ** 2 + (y_te - center_y) ** 2) # distance of each point from center
outside_circle_te = distance_te >= radius # filter out points inside circle
x_te_filtered = x_te[outside_circle_te]
y_te_filtered = y_te[outside_circle_te]
test_points = np.vstack((x_te_filtered, y_te_filtered)).T # stack x, y together
test_points = torch.tensor(test_points, dtype=torch.float32).to(device) # convert to tensor and move training points to GPU / CPU, use float64 if more precision is required
predict = model(test_points).detach().cpu().numpy() # predict and move back to CPU
u = predict[:, 0] * Uinf # predict and bring back to dimensional form
Thank you for reading! If you want to hire me as your next PostDoc researcher, reach out for collaboration in research, please do so!
References
[1] U. Ghia, K. N. Ghia, and C. T. Shin, “High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method”, Journal of Computational Physics 48, 387-411 (1982)
[2] Cahya Amalinadhi, Pramudita S. Palar, Rafael Stevenson and Lavi Zuhal. "On Physics-Informed Deep Learning for Solving Navier-Stokes Equations," AIAA 2022-1436. AIAA SCITECH 2022 Forum. January 2022. https://doi.org/10.2514/6.2022-1436
No comments:
Post a Comment