Motivation¶
The three main methods used to solve scientific equations and problems are analytical, numerical, and experimental. However, we cannot always apply the experimental technique due to financial and time restrictions. Analytical techniques are the classic approaches to problem-solving. However, we cannot solve equations analytically due to limitations imposed by complex geometry, boundary conditions, and other factors.
As a result, we have been encouraging numerical approaches for several years since they may generate answers that are virtually as trustworthy as those of analytical methods in a lot shorter and simpler amount of time. For instance, the well-known Navier-Stoke equation may be swiftly solved using Numerical Schemes even though it has never been analytically solved. On the other hand, the drawbacks of conventional numerical schemes are their difficulty in expressing irregular borders, the fact that they are not designed for unstructured meshes, and the fact that momentum, energy, and mass are not conserved. Lastly, they are oriented toward edges and one-dimensional physics and unsuitable for substantial turbulent flow issues.
A revolutionary deep learning paradigm called physics-informed neural networks (PINNs) is well suited for tackling forward and inverse issues of nonlinear partial differential equations (PDEs). Feedforward neural networks are trained as surrogate models for approximation solutions to the PDEs without the requirement for label data by including the physical information defined by PDEs in them. Numerous PINN-based approaches have been created to solve various problems, including integer-order PDEs, fractional PDEs, stochastic PDEs, and integrodifferential equations. This is due to neural networks' incredible capacity to describe complicated relationships.
In this article, I investigate the usage of PINNs as a solution approximation for 1D Burgers' Equation. I first solve the equation using Finite Difference Method and the Solve the exact same equation using a PINN.
A Simple 1D Burgers' Equation:¶
You can read about Burgers' Equation on its wikipedia page.
Burgers' equation in one spatial dimension looks like this:
$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial ^2u}{\partial x^2}$$As you can see, it is a combination of non-linear convection and diffusion. It is surprising how much you learn from this neat little equation!
$$x\in[0,2]$$$$t\in[0,0.48]$$$$\nu → Viscosity = 0.01/\pi $$Initial Condition:¶
$$y(x,0)=sin(\pi x)$$Boundary Conditions:¶
$$y(-1,t)=0$$$$y(1,t)=0$$Finite Difference Method¶
The second-order derivative can be represented geometrically as the line tangent to the curve given by the first derivative. We will discretize the second-order derivative with a Central Difference scheme: a combination of Forward Difference and Backward Difference of the first derivative. Consider the Taylor expansion of $u_{i+1}$ and $u_{i-1}$ around $u_i$:
$u_{i+1} = u_i + \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i + \frac{\Delta x^3}{3!} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)$
$u_{i-1} = u_i - \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i - \frac{\Delta x^3}{3!} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)$
If we add these two expansions, you can see that the odd-numbered derivative terms will cancel each other out. If we neglect any terms of $O(\Delta x^4)$ or higher (and really, those are very small), then we can rearrange the sum of these two expansions to solve for our second-derivative.
$u_{i+1} + u_{i-1} = 2u_i+\Delta x^2 \frac{\partial ^2 u}{\partial x^2}\bigg|_i + O(\Delta x^4)$
Then rearrange to solve for $\frac{\partial ^2 u}{\partial x^2}\bigg|_i$ and the result is:
$$\frac{\partial ^2 u}{\partial x^2}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2} + O(\Delta x^2)$$We can now write the discretized version of the Burgers equation in 1D:
$$ \frac{u^{n+1}_i - u^n_i}{\Delta t} + u_i^n \frac{u^{n}_i - u^n_{i-1}}{\Delta x} = \nu \frac{u^{n}_{i+1} -2u^n_i + u^n_{i-1}}{\Delta x^2} $$As before, we notice that once we have an initial condition, the only unknown is $u_{i}^{n+1}$, so we re-arrange the equation solving for our unknown:
$$ u^{n+1}_i = u^n_i - u^n_i \frac{\Delta t}{\Delta x} (u^n_i - u^n_{i-1}) + \frac{\nu \Delta t}{\Delta x^2}(u^{n}_{i+1} - 2u^n_i + u^n_{i-1}) $$! pip install pyDOE --quiet
from pyDOE import lhs
#We will use Latin Hypercube Sampling from this library
import torch # Pytorch
import torch.autograd as autograd # computation graph
from torch import Tensor
import torch.nn as nn # neural networks
import torch.optim as optim # optimizers
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
import time
import pandas as pd
import matplotlib.pyplot as plt
import time, sys
from tqdm.notebook import tqdm_notebook
#Set default dtype to float32
torch.set_default_dtype(torch.float)
#PyTorch random number generator
torch.manual_seed(1234)
# Random number generators in other libraries
np.random.seed(1234)
# Device configuration
device = torch.device('cuda' if torch.cuda.is_available() else 'cgpu')
print("The neural network will be trainied on",device)
if device == 'cuda':
print(torch.cuda.get_device_name())
%matplotlib inline
%config InlineBackend.print_figure_kwargs = {'bbox_inches':None}
def plot3D(x,t,y):
x_plot =x.squeeze(1)
t_plot =t.squeeze(1)
X,T= torch.meshgrid(x_plot,t_plot,indexing='ij')
u_xt = y
fig = plt.figure()
ax=fig.subplots(1,1)
cp = ax.contourf(T,X,u_xt,20,cmap=cm.rainbow) #)levels = np.linspace(-1.0,1.0,12))
fig.colorbar(cp) # Add a colorbar to a plot
ax.set_title('u')
ax.set_xlabel('t')
ax.set_ylabel('x')
plt.show()
ax = plt.axes(projection='3d')
surf = ax.plot_surface(T.numpy(), X.numpy(), u_xt.numpy(),cmap=cm.rainbow, antialiased=False)
ax.set_xlabel('t')
ax.set_ylabel('x')
ax.set_zlabel('u')
#ax.set_zlim3d(-1, 1)
plt.show()
#x∈[0,2]
x_min=0
x_max=2
#t∈[0,0.48]
t_min=0
t_max=0.48
viscosity = 0.01/np.pi # Given
#Discretization points for x and t
total_points_x=1001
total_points_t=1000
dx = (x_max-x_min)/(total_points_x-1)
dt = (t_max-t_min)/(total_points_t)
dx,dt, dt/dt
#Implementing Finite Difference Method to solve the 1D Diffusion Equation
def u_fem(x,t):
un = torch.ones(total_points_x)
rec = torch.zeros([total_points_x, total_points_t])
for j in tqdm_notebook(range(total_points_t)):
un = u.clone()
for i in range(1,total_points_x-1):
rec[i,j] = u[i]
u[i] = un[i] - un[i] * dt/dx * (un[i]-un[i-1]) + viscosity * (dt/dx**2) * (un[i+1]- 2*un[i] + un[i-1])
if np.isnan(u[i]):
print(i, j, u[i-1])
break
return u, rec
x = torch.linspace(x_min, x_max, total_points_x)
u = torch.from_numpy(np.sin(np.pi*x.numpy()))
u_plot = u.clone()
plt.figure(figsize=(11, 7), dpi= 100)
plt.plot(x, u_plot, lw=2)
plt.xlabel('x')
plt.ylabel('u')
plt.title('Burgers Equation at t=0');
# Creating same amount of grid lattice as FDM
x = torch.linspace(x_min, x_max, total_points_x).view(-1,1)
t = torch.linspace(t_min, t_max, total_points_t).view(-1,1)
x.shape, t.shape
# Computing Finite Difference solution for u
print("Running Finite Difference Method...")
u_final, u_fem_2D = u_fem(x,t)
assert u_fem_2D.shape == torch.Size([total_points_x, total_points_t]),f"Expected [{total_points_x},{total_points_t}], got {u_fem_2D.shape}"
print("Completed successfully!")
Finite Difference Method - Solution¶
X, T = torch.meshgrid(x.squeeze(1),t.squeeze(1), indexing='ij')
X.shape, T.shape
plt.figure(figsize=(11, 7), dpi= 100)
plt.plot(x, u_plot, label="Initial Condition, t=0")
plt.plot(x, u_final, 'r',label="FDM, t=0.48")
plt.xlim([x_min,x_max])
plt.ylim([-1, 1])
plt.legend()
plt.xlabel('x')
plt.ylabel('u')
plt.title('Burgers Equation with Finite Difference Method');
plt.show()
print("FDM Solution Visualization")
plot3D(x,t,u_fem_2D)
Notice that as time approached 0.48, A shock wave id forming at the center. Due to limtations of the Finite difference method, aythong beyond t=0.48 will result in divergence. Shock wave is expected to form at t =0.5.
Physics Informed Neural Network¶
If we rearrange our PDE, we get:
$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} - \nu \frac{\partial ^2u}{\partial x^2}=0$$Neural Network¶
A Neural Network is a function approximator constructed as follows:
$$NN(X)=W_n\sigma_{n-1}(W_{n-1}\sigma_{n-2}(...(W_2(W_1X+b_1)+b_2)+..)+b_{n-1})+b_n$$Note: We usually train our NN by iteratively minimizing a loss function ($MSE$:mean squared error) in the training dataset(known data).
PINNs = Neural Network + PDE¶
We can use a neural network to approximate any function (Universal APproximation Theorem): $$NN(x,t)\approx u(x,t)$$
Since NN is a function, we can obtain its derivatives: $\frac{\partial NN}{\partial t},\frac{\partial^2 NN}{\partial x^2}$.(Automatic Diferentiation)
Assume:$$NN(t,x)\approx u(t,x)$$
Then:
$$\frac{\partial NN}{\partial t}+NN\frac{\partial NN}{\partial x}-\nu\frac{\partial^2 NN}{\partial x^2}\approx \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} - \nu \frac{\partial ^2u}{\partial x^2}=0$$And:
$$\frac{\partial NN}{\partial t}+NN\frac{\partial NN}{\partial x}-\nu\frac{\partial^2 NN}{\partial x^2}\approx 0$$We define this function as $f$:
$$f(t,x)=\left( \frac{\partial NN}{\partial t}+NN\frac{\partial NN}{\partial x}-\nu\frac{\partial^2 NN}{\partial x^2}\right)\rightarrow residue$$If $residue\rightarrow 0$ then our NN would be respecting the physical law. To minimize the $residue$ we will use L2 norm with matrix of zeros with same size as the $residue$
PINNs' Loss function¶
We evaluate our PDE in a certain number of "collocation points" ($N_f$) inside our domain $(x,t)$. Then we iteratively minimize a loss function related to $f$:
$$MSE_f=\frac{1}{N_f}\sum^{N_f}_{i=1}|f(t_f^i,x_f^i)|^2$$Usually, the training data set is a set of points from which we know the answer. In our case, we will use our boundary(BC) and initial conditions(IC).
Initial Condition:¶
$$y(x,0)=sin(\pi x)$$Boundary Conditions:¶
$$y(-1,t)=0$$$$y(1,t)=0$$Since we know the outcome, we select $N_u$ points from our BC and IC and used them to train our network.
$$MSE_{u}=\frac{1}{N_u}\sum^{N_u}_{i=1}|y(t_{u}^i,x_u^i)-NN(t_{u}^i,x_u^i)|^2$$Total Loss:¶
$$MSE=MSE_{u}+MSE_f$$#PINNs
# Creating same amount of grid lattice as FDM
x = torch.linspace(x_min, x_max, total_points_x).view(-1,1)
t = torch.linspace(t_min, t_max, total_points_t).view(-1,1)
x.shape, t.shape
X, T = torch.meshgrid(x.squeeze(1),t.squeeze(1), indexing='ij') #same as FDM
X.shape, T.shape
left_X = torch.hstack((X[:,0][:,None], T[:,0][:,None])) #horizontal stacking to create X, T dataset
left_U = torch.sin(np.pi*left_X[:,0]).unsqueeze(1) #initial condition is a sine wave
left_U.shape
plt.figure(figsize=(11, 7), dpi= 100)
plt.plot(x, left_U, lw=2)
plt.xlabel('x')
plt.ylabel('u')
plt.title('Burgers Equation at t=0');
# BC at x_min
bottom_X = torch.hstack((X[0,:][:,None],T[0,:][:,None]))
top_X = torch.hstack((X[-1,:][:,None],T[-1,:][:,None]))
bottom_U = torch.zeros(bottom_X.shape[0],1)
top_U = torch.zeros(top_X.shape[0],1)
bottom_X.shape
X_bc = torch.vstack([bottom_X, top_X])
U_bc = torch.vstack([bottom_U, top_U])
X_bc.shape
N_ic = 1000
N_bc = 1000 #Number of points on IC and BC
N_pde = 30000 #Number of points on PDE domain (Collocation Points)
#Now we will sample N_bc points at random
#from the X_train, U_train dataset
idx = np.random.choice(X_bc.shape[0],N_bc, replace=False)
X_bc_samples = X_bc[idx,:]
U_bc_samples = U_bc[idx,:]
idx = np.random.choice(left_X.shape[0],N_ic, replace=False)
X_ic_samples = left_X[idx,:]
U_ic_samples = left_U[idx,:]
#The boundary conditions will not change.
#Hence, these U values can be used as supervised labels during training
#For PDE collocation points, we will generate new X_train_pde dataset
#We do not know U(X,T) for these points
#Lets get the entire X,T dataset in a format suitable for Neural Network
#We will later use this for testing NN as well. So, lets call this x_test for convenience
x_test = torch.hstack((X.transpose(1,0).flatten()[:,None],
T.transpose(1,0).flatten()[:,None]))
#We need column major flattening to simlulte time-marching. Hence the transpose(1,0) or simply use .T
#we will use U generated from FEM as our u_test
#We will use u_test later in the process for calculating NN performance
u_test = u_fem_2D.transpose(1,0).flatten()[:,None]
x_test.shape
x_test.shape
#lower and upper bounds of x_test
lb = x_test[0]
ub = x_test[-1]
lb,ub
#Sampling (X,T) domain using LHS
lhs_samples = lhs(2,N_pde)
#2 since there are 2 variables in X_train, [x,t]
lhs_samples.shape
X_train_lhs = lb + (ub-lb)*lhs_samples
X_train_lhs.shape
plt.plot(pd.DataFrame(X_train_lhs)[0], "bo", markersize=.5)
plt.title("LHS Sampling of X [0,2]")
plt.show()
plt.plot(pd.DataFrame(X_train_lhs)[1], "ro", markersize=.5)
plt.title("LHS Sampling of T [0,0.48]")
plt.show()
X_train_final = torch.vstack((X_train_lhs, X_ic_samples, X_bc_samples))
X_train_final.shape
#Lets define a u_NN
class u_NN(nn.Module):
def __init__(self, layers_list):
super().__init__()
self.depth = len(layers_list)
self.loss_function = nn.MSELoss(reduction="mean")
self.activation = nn.Tanh() #This is important, ReLU wont work
self.linears = nn.ModuleList([nn.Linear(layers_list[i],layers_list[i+1]) for i in range(self.depth-1)])
for i in range(self.depth-1):
nn.init.xavier_normal_(self.linears[i].weight.data, gain=1.0) #xavier normalization of weights
nn.init.zeros_(self.linears[i].bias.data) #all biases set to zero
def Convert(self, x): #helper function
if torch.is_tensor(x) !=True:
x = torch.from_numpy(x)
return x.float().to(device)
def forward(self, x):
a = self.Convert(x)
for i in range(self.depth-2):
z = self.linears[i](a)
a = self.activation(z)
a = self.linears[-1](a)
return a
def loss_bc(self, x_bc, u_bc):
#This is similar to a Supervised Learning
l_bc = self.loss_function(self.forward(self.Convert(x_bc)), self.Convert(u_bc)) #L2 loss
return l_bc
def loss_ic(self, x_ic, u_ic):
#This is similar to a Supervised Learning
l_ic = self.loss_function(self.forward(self.Convert(x_ic)), self.Convert(u_ic)) #L2 loss
return l_ic
def loss_pde(self, x_pde):
# We will pass x_train_final here.
# Note that we do not have U_pde (labels) here to calculate loss. This is not Supervised Learning.
# Here we want to minimize the residues. So, we will first calculate the residue and then minimize it to be close to zero.
x_pde = self.Convert(x_pde)
x_pde_clone = x_pde.clone() ##VERY IMPORTANT
x_pde_clone.requires_grad = True #enable Auto Differentiation
NN = self.forward(x_pde_clone) #Generates predictions from u_NN
NNx_NNt = torch.autograd.grad(NN, x_pde_clone,self.Convert(torch.ones([x_pde_clone.shape[0],1])),retain_graph=True, create_graph=True)[0] #Jacobian of dx and dt
NNxx_NNtt = torch.autograd.grad(NNx_NNt,x_pde_clone, self.Convert(torch.ones(x_pde_clone.shape)), create_graph=True)[0] #Jacobian of dx2, dt2
NNxx = NNxx_NNtt[:,[0]] #Extract only dx2 terms
NNt = NNx_NNt[:,[1]] #Extract only dt terms
NNx = NNx_NNt[:,[0]] #Extract only dx terms
# {(du/dt) = viscosity * (d2u/dx2)} is the pde and the NN residue will be {du_NN/dt - viscosity*(d2u_NN)/dx2} which is == {NNt - viscosity*NNxx}
residue = NNt + self.forward(x_pde_clone)*(NNx) - (viscosity)*NNxx
# The residues need to be zero (or as low as possible). We'll create an arrazy of Zeros and minimize the residue
zeros = self.Convert(torch.zeros(residue.shape[0],1))
l_pde = self.loss_function(residue, zeros) #L2 Loss
return l_pde
def total_loss(self, x_ic, u_ic, x_bc, u_bc, x_pde): #Combine both loss
l_bc = self.loss_bc(x_bc, u_bc)
l_ic = self.loss_ic(x_ic, u_ic)
l_pde = self.loss_pde(x_pde)
return l_bc + l_pde + l_ic #this HAS to be a scalar value for auto differentiation to do its thing.
#Parameters for u_NN
EPOCHS = 100000
initial_lr = 0.001
layers_list = [2, 32, 128, 16, 128, 32, 1]
#batch_size = 32
# Instantiate a model
PINN = u_NN(layers_list).to(device)
print(PINN)
optimizer = torch.optim.Adam(PINN.parameters(), lr=initial_lr,amsgrad=False)
scheduler = torch.optim.lr_scheduler.ExponentialLR(optimizer, gamma=0.985)
history = pd.DataFrame(columns=["Epochs","Learning_Rate", "IC_Loss","BC_Loss","PDE_Loss","Total_Loss","Test_Loss"])
#****** Training ******#
print("Training Physics Informed Neural Network...")
Epoch = []
Learning_Rate = []
IC_Loss = []
BC_Loss = []
PDE_Loss = []
Total_Loss = []
Test_Loss = []
for i in tqdm_notebook(range(EPOCHS)):
if i==0:
print("Epoch \t Learning_Rate \t IC_Loss \t BC_Loss \t PDE_Loss \t Total_Loss \t Test_Loss")
l_ic = PINN.loss_ic(X_ic_samples,U_ic_samples)
l_bc = PINN.loss_bc(X_bc_samples,U_bc_samples)
l_pde = PINN.loss_pde(X_train_final)
loss = PINN.total_loss(X_ic_samples,U_ic_samples,X_bc_samples,U_bc_samples, X_train_final)
optimizer.zero_grad()
loss.backward()
optimizer.step()
if i%100 == 0: #print losses and step the exponential learning rate.
with torch.no_grad():
test_loss = PINN.loss_bc(x_test,u_test) #Here we are using loss_bc method as a helper function to calculate L2 loss
Epoch.append(i)
Learning_Rate.append(scheduler.get_last_lr()[0])
IC_Loss.append(l_ic.detach().cpu().numpy())
BC_Loss.append(l_bc.detach().cpu().numpy())
PDE_Loss.append(l_pde.detach().cpu().numpy())
Total_Loss.append(loss.detach().cpu().numpy())
Test_Loss.append(test_loss.detach().cpu().numpy())
if i%1000 ==0:
print(i,'\t',format(scheduler.get_last_lr()[0],".4E"),'\t',format(l_ic.detach().cpu().numpy(),".4E"),'\t',format(l_bc.detach().cpu().numpy(),".4E"),'\t',
format(l_pde.detach().cpu().numpy(),".4E"),'\t',format(loss.detach().cpu().numpy(),".4E"),'\t',format(test_loss.detach().cpu().numpy(),".4E"))
scheduler.step()
print("Completed!!")
Neural Network - Solution¶
# Predictions from the NN
u_NN_predict = PINN(x_test)
#Reshapping y1 to be used in plot3d()
u_NN_2D = u_NN_predict.reshape(shape=[total_points_t,total_points_x]).transpose(1,0).detach().cpu()
assert u_NN_2D.shape == torch.Size([total_points_x, total_points_t]),f"Expected [{total_points_x},{total_points_t}], got {u_NN_2D.shape}"
RMSE = torch.sqrt(torch.mean(torch.square(torch.subtract(u_NN_2D,u_fem_2D))))
print("The RMSE error between FDM and PINN is :",np.around(RMSE.item(),5))
file_name_loss = "LOSS_CURVES_RMSE_"+str(np.around(RMSE.item(),5))+"_Epochs_"+str(EPOCHS)+"_lr_"+str(initial_lr)+".png"
fig, ax1 = plt.subplots(figsize=(11, 7), dpi= 100)
ax2 = ax1.twinx()
ax1.plot(Epoch, IC_Loss, "b-",label = "IC Loss")
ax1.plot(Epoch, BC_Loss, "g-",label = "BC Loss")
ax1.plot(Epoch, PDE_Loss, "c-",label = "PDE Loss")
ax1.plot(Epoch, Total_Loss, "k-",label = "Total Loss")
ax1.plot(Epoch, Test_Loss, "m-",label = "Test Loss")
ax2.plot(Epoch,Learning_Rate, "ro",markersize=1,label = "Learning Rate")
ax1.set_xlabel('Epochs')
ax1.set_ylabel('Losses', color='k')
ax2.set_ylabel('Learning Rate', color='k')
ax2.legend(loc=7)
ax1.legend(loc=1)
plt.title(file_name_loss)
#plt.show()
plt.savefig("/content/drive/MyDrive/Colab Notebooks/PINNs/Burgers/"+file_name_loss)
3D plot of solution for U based on Neural Network (PINN)
plot3D(x,t,u_NN_2D)
3D plot of solution for U based on Finite Difference Method computed earlier.
plot3D(x,t,u_fem_2D)
3D plot of Error between Finite Difference Method and Physics Informed Neural Network model.
plot3D(x,t,(u_NN_2D - u_fem_2D)) #Error
file_name_result = "Result_RMSE_"+str(np.around(RMSE.item(),5))+"_Epochs_"+str(EPOCHS)+"_lr_"+str(initial_lr)+".png"
last_U_NN = u_NN_2D[:,-1].unsqueeze(1) #Extracting the last U values at t=0.48
fig, ax1 = plt.subplots(figsize=(11, 7), dpi= 100)
ax2 = ax1.twinx()
ax1.plot(x, u_plot, "g",label="Initial Condition at t=0")
ax1.plot(x, u_final, "r",label="FDM at t=0.48", )
ax1.plot(x, last_U_NN.detach().cpu().numpy(), "bo",label="PINN Predict at t=0.48", markersize=1.5)
ax1.legend()
plt.title(file_name_result)
plt.savefig("/content/drive/MyDrive/Colab Notebooks/PINNs/Burgers/"+file_name_result)
#Saving the model
file_name_model = "Result_RMSE_"+str(np.around(RMSE.item(),5))+"_Epochs_"+str(EPOCHS)+"_lr_"+str(initial_lr)+".pth"
torch.save(PINN, "/content/drive/MyDrive/Colab Notebooks/PINNs/Burgers/"+file_name_model)
Conclusion¶
The PINNs perform exceptionally well and are much more stable then Finite Difference Method especially at the shock region around t=0.48. There is some "error" between the FDM and PINNs at the shock region. Further assessment and investigation is necessary.
References¶
[1] Raissi, M., Perdikaris, P., & Karniadakis, G. E. (2017). Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations. arXiv preprint arXiv:1711.10561. http://arxiv.org/pdf/1711.10561v1
[2] https://www.iist.ac.in/sites/default/files/people/IN08026/Burgers_equation_viscous.pdf
[3] https://en.wikipedia.org/wiki/Burgers%27_equation
[4] Rackauckas Chris, Introduction to Scientific Machine Learning through Physics-Informed Neural Networks. https://book.sciml.ai/notes/03/
[5] https://github.com/jdtoscano94/Learning-Python-Physics-Informed-Machine-Learning-PINNs-DeepONets
[6] Kashefi, Ali and Tapan Mukerji. “Physics-informed PointNet: A deep learning solver for steady-state incompressible flows and thermal fields on multiple sets of irregular geometries.” J. Comput. Phys. 468 (2022): 111510.
#################### ANIMATION ########################
#Imports for animation and display within a jupyter notebook
from matplotlib import animation, rc
from IPython.display import HTML
#Generating the figure that will contain the animation
fig, ax = plt.subplots()
fig.set_dpi(100)
fig.set_size_inches(6, 5)
ax.set_xlim(( x_min, x_max))
ax.set_ylim((-1.2, 1.2))
PINN, = ax.plot([], [], "bo", markersize=3,label='Neural Network')
FDM, = ax.plot([], [], "r", lw=2,label='Finite Difference')
ax.legend();
plt.xlabel('x')
plt.ylabel('u')
plt.title('1D Burgers Equation Time Evolution from t=0 to t=0.48');
#Initialization function for funcanimation
def init():
FDM.set_data([], [])
PINN.set_data([], [])
return (PINN, FDM)
#Main animation function, each frame represents a time step in our calculation
def animate(j):
unn = u_NN_2D[:,j].unsqueeze(1).detach().cpu().numpy()
ufdm = u_fem_2D[:,j].unsqueeze(1).detach().cpu().numpy()
PINN.set_data(x, unn)
FDM.set_data(x, ufdm)
return (PINN, FDM,)
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=total_points_t, interval=100)
anim.save('/content/drive/MyDrive/Colab Notebooks/PINNs/Burgers/1dBurgers_06102022_1300.gif',animation.PillowWriter(fps=60))
#HTML(anim.to_jshtml())
#jupyter nbconvert --to html <filename.ipynb>