Physics-informed Neural Networks (PINN)
http://iailab.kaist.ac.kr/
Industrial AI Lab at KAIST
Table of Contents
1. Why Deep Learning Needs Physics?¶
(1) Advantages of Data-Driven Approaches
When sufficient data is available, data-driven models can achieve a high level of predictive accuracy without requiring extensive domain knowledge. These methods leverage large-scale datasets to learn complex patterns, making them effective for a wide range of applications.
(2) Challenges of Black-Box Data-Driven Methods
Despite their predictive capabilities, purely data-driven models, particularly black-box approaches, face several limitations:
They may produce physically inconsistent results.
They are prone to learning spurious correlations that appear valid only within the training and test datasets.
They often exhibit poor generalization, particularly when applied to out-of-sample prediction tasks.
The lack of interpretability hinders scientific understanding and trust in the model's predictions.
Understanding the underlying mechanisms governing a system is essential for scientific discovery and real-world applicability.
(3) Limitations of Purely Data-Driven Approaches
Limited Generalizability: Models struggle with out-of-distribution data, reducing reliability in real-world scenarios.
Lack of Interpretability: Predictions often lack physical meaning, making it difficult to extract scientifically relevant insights.
(4) Physics-informed Neural Networks (PINNs)
To address these challenges, Physics-informed Neural Networks (PINNs) integrate scientific domain knowledge into deep learning models.
Combine the strengths of data-driven methods with established physical laws, improving predictive accuracy and robustness.
Help mitigate issues related to imbalanced datasets and data scarcity by enforcing known physical constraints.
Enhance interpretability by ensuring outputs align with fundamental scientific principles, making them more reliable for real-world applications.
By embedding prior knowledge from physics and engineering, PINNs provide a more generalizable, interpretable, and physically consistent framework for deep learning, overcoming the limitations of purely data-driven methods.
Deep Neural Networks:
Each type of DNN architecture is tailored to handle specific data structures and tasks:
ANN (Artificial Neural Networks): Basic neural network structures that form the foundation of deep learning.
CNN (Convolutional Neural Networks): Designed primarily for tasks like image recognition and feature extraction.
RNN (Recurrent Neural Networks): Specialized for sequential data processing, such as time series or natural language tasks.
GNN (Graph Neural Networks): Utilized for processing graph-structured data, like social networks or molecular structures.
Generative Models: Models like GANs and VAEs for creating new data samples from learned distributions.
Knowledge Representation:
Knowledge representation refers to how domain-specific knowledge is encoded to be accessible and actionable for Deep Neural Networks. The choice of representation plays a critical role in enhancing learning efficiency and interpretability:
Differential Equations: Used for modeling continuous dynamical systems.
Algebraic Equations: Represent mathematical relationships and constraints.
Knowledge Graphs: Capture relational and semantic information between entities.
Simulation Results: Provide structured outputs from computational simulations.
Human Feedback: Incorporates expert insights and manual annotations.
Knowledge Integration:
Integrating external knowledge into DNNs is a critical step to ensure that models not only rely on large datasets but also leverage existing domain knowledge. This improves model accuracy, generalizability, and interpretability:
Feature Engineering: Tailoring input features to better align with domain-specific knowledge.
Designing: Modifying model architectures to incorporate external knowledge sources.
Regularizing: Adding constraints, priors, or penalties to the loss function ensures that the model adheres to known properties of the domain. This prevents overfitting and encourages the model to find solutions consistent with real-world knowledge.
1.2. Multilayer Feedforward Networks are Universal Approximators¶
The Universal Approximation Theorem (1989)
- Neural Networks are capable of approximating any Borel measurable function
- A foundational result in the theory of neural networks, stating that a feedforward neural network with at least one hidden layer and a finite number of neurons can approximate any continuous function to arbitrary accuracy, given sufficient neurons and an appropriate activation function.
1.3. Neural Networks for Solving Differential Equations¶
Neural networks have emerged as a powerful tool for solving differential equations, offering a flexible and efficient alternative to traditional numerical methods. By leveraging their ability to approximate complex functions, neural networks can model solutions to both ordinary differential equations (ODEs) and partial differential equations (PDEs). This approach is particularly advantageous for high-dimensional or nonlinear systems.
Two foundational works on the application of neural networks for solving differential equations are presented here.
- Neural Algorithm for Solving Differential Equations (1990)
- ANN for ODE and PDE (1998)
1.4. Journal of Computational Physics (2019)¶
The following paper has played a significant role in advancing the field of Physics-informed Neural Networks (PINNs). It is highly recommended for those interested in this emerging area of research:
This groundbreaking work explores how deep learning can be leveraged to solve complex forward and inverse problems in nonlinear partial differential equations, providing a foundation for further advancements in combining deep learning with domain-specific knowledge.
1.5. Nature Reviews Physics (2021)¶
The following paper explores several prevailing trends in embedding physics into machine learning, highlighting current capabilities, limitations, and diverse applications of physics-informed learning. This includes its use in solving both forward and inverse problems, uncovering hidden physical laws, and addressing high-dimensional challenges.
2. Physics-informed Neural Networks (PINNs) with PDE, IC, BC¶
Physics-informed Neural Networks (PINNs) integrate the laws of physics, represented by partial differential equations (PDEs), along with initial and boundary conditions into the training process of neural networks. This approach makes them powerful for solving forward and inverse problems, even with limited or noisy data.
Key Formulation
Consider a generic time-dependent PDE defined on a domain $ \Omega \times [0, T] $:
$$ \mathcal{N}[u(x, t)] = f(x, t), \quad (x, t) \in \Omega \times [0, T], $$
where:
- $ \mathcal{N} $ is a differential operator involving derivatives of $ u(x, t) $,
- $ u(x, t) $ is the unknown solution to be approximated,
- $ f(x, t) $ is a source or forcing term,
- $ \Omega $ is the spatial domain, and $[0, T]$ is the temporal interval.
Initial and Boundary Conditions
To uniquely solve the PDE, we include:
(1) Initial Condition: Specifies the solution at the start of the time domain:
$$ u(x, 0) = g(x), \quad x \in \Omega. $$
(2) Boundary Conditions: Specifies the behavior of the solution at the spatial domain boundaries. For example:
- Dirichlet Boundary Condition:
$$ u(x_b, t) = h_1(x_b, t), \quad x_b \in \partial \Omega, \, t \in [0, T]. $$
- Neumann Boundary Condition:
$$ \frac{\partial u}{\partial n}(x_b, t) = h_2(x_b, t), \quad x_b \in \partial \Omega, \, t \in [0, T]. $$
PINN Loss Function
PINNs approximate the solution $ u(x, t) $ using a neural network $ u_\theta(x, t) $, parameterized by $ \theta $. The training involves minimizing a total loss function composed of the following components:
(1) Physics Equation (PDE):
$$ \mathcal{L}_{\text{physics}} = \frac{1}{N_f} \sum_{i=1}^{N_f} \left| \mathcal{N}[u_\theta(x_i, t_i)] - f(x_i, t_i) \right|^2, $$
where $ \{(x_i, t_i)\}_{i=1}^{N_f} $ are collocation points sampled from the domain $ \Omega \times [0, T] $.
(2) Initial Condition:
$$ \mathcal{L}_{\text{initial}} = \frac{1}{N_I} \sum_{i=1}^{N_I} \left| u_\theta(x_i, 0) - g(x_i) \right|^2, $$
where $ \{x_i\}_{i=1}^{N_I} $ are points sampled from $ \Omega $ at $ t = 0 $.
(3) Boundary Conditions:
For Dirichlet boundary conditions:
$$ \mathcal{L}_{\text{Dirichlet}} = \frac{1}{N_B} \sum_{i=1}^{N_B} \left| u_\theta(x_b, t_i) - h_1(x_b, t_i) \right|^2, $$
where $ \{(x_b, t_i)\}_{i=1}^{N_B} $ are points sampled on the boundary $ \partial \Omega $.
For Neumann boundary conditions:
$$ \mathcal{L}_{\text{Neumann}} = \frac{1}{N_B} \sum_{i=1}^{N_B} \left| \frac{\partial u_\theta}{\partial n}(x_b, t_i) - h_2(x_b, t_i) \right|^2. $$
(4) Total Loss Function:
The total loss function combines all components:
$$ \mathcal{L}_{\text{total}} = \mathcal{L}_{\text{physics}} + \lambda_I \mathcal{L}_{\text{initial}} + \lambda_B \mathcal{L}_{\text{boundary}}, $$
where $ \lambda_I $ and $ \lambda_B $ are weights to balance the contributions of initial and boundary losses.
By encoding physical laws and constraints directly into the learning process, PINNs offer a robust and data-efficient framework for solving challenging real-world problems.
3. Inverse Problems using Observed Data¶
We investigate the potential of Physics-informed Neural Networks (PINNs) as a neural network-based approach for solving forward problems. However, this method does not offer significant advantages over conventional computational techniques.
In contrast, the application of PINNs for inverse problems presents a promising opportunity for advancing scientific computing.
Inverse problems involve determining unknown parameters, functions, or underlying physical laws based on limited information. Given the typically ill-posed nature of inverse problems, they often suffer from insufficient information, making it challenging to obtain unique and stable solutions. To address this issue, a combination of observed data and physical laws is essential.
PINNs provide a robust framework for solving inverse problems by integrating both data-driven insights and physics-based constraints derived from governing equations. This hybrid approach enhances model stability and accuracy, making PINNs a powerful tool for tackling complex inverse problems in scientific computing.
Let's examine how PINNs can be used to estimate unknown parameters, serving as an example of solving inverse problems.
Formulation of an Inverse Problem with PINNs
Given a system governed by a partial differential equation (PDE):
$$ \mathcal{N}[u(x, t; \lambda)] = f(x, t), \quad (x, t) \in \Omega \times [0, T], $$
where:
- $ \mathcal{N} $ is a differential operator acting on the unknown solution $ u(x,t) $,
- $ \lambda $ is an unknown parameter or function to be inferred,
- $ f(x, t) $ is a known forcing term,
- $ \Omega $ is the spatial domain, and $ [0, T] $ is the temporal interval.
The goal is to estimate $ \lambda $, given a set of observed data points:
$$ \{ (x_i, t_i, u_{\text{obs}}(x_i, t_i) ) \}_{i=1}^{N_d}. $$
PINNs solve this problem by learning both the solution $ u_\theta(x,t) $ and the unknown parameter $ \lambda $ simultaneously.
Loss Function for the Inverse Problem
To estimate $ \lambda $, we define a composite loss function that minimizes:
(1) Physics Residual Loss: Ensures the PDE is satisfied at collocation points:
$$ \mathcal{L}_{\text{physics}} = \frac{1}{N_f} \sum_{i=1}^{N_f} \left| \mathcal{N}[u_\theta(x_i, t_i; \lambda)] - f(x_i, t_i) \right|^2. $$
(2) Initial Condition Loss: Enforces the initial condition of the system:
$$ \mathcal{L}_{\text{initial}} = \frac{1}{N_I} \sum_{i=1}^{N_I} \left| u_\theta(x_i, 0) - g(x_i) \right|^2. $$
(3) Boundary Condition Loss: Enforces boundary conditions:
$$ \mathcal{L}_{\text{boundary}} = \frac{1}{N_B} \sum_{i=1}^{N_B} \left| u_\theta(x_b, t_i) - h_1(x_b, t_i) \right|^2. $$
(4) Data Loss: Ensures the predicted solution matches the observed data:
$$ \mathcal{L}_{\text{data}} = \frac{1}{N_d} \sum_{i=1}^{N_d} \left| u_\theta(x_i, t_i) - u_{\text{obs}}(x_i, t_i) \right|^2. $$
(5) Total Loss Function: Combining all components:
$$ \mathcal{L}_{\text{total}} = \mathcal{L}_{\text{physics}} + \lambda_I \mathcal{L}_{\text{initial}} + \lambda_B \mathcal{L}_{\text{boundary}} + \lambda_D \mathcal{L}_{\text{data}}. $$
where $ \lambda_I, \lambda_B, \lambda_D $ are weights that balance the contributions of each loss term.
By combining observed data with known physics, PINNs provide a powerful framework for solving inverse problems across various scientific and engineering domains.
4. PINN Implementation¶
We have discussed the fundamental theory behind Physics-informed Neural Networks (PINNs); now, let's explore how PINNs are implemented in practice. This will involve constructing the neural network, defining the loss function incorporating physical constraints, and training the model to solve forward problems.
4.1. Lab 1: Simple Example¶
Here, we will begin with the simplest example of an ordinary differential equation (ODE), which will be implemented from scratch using Keras.
- Let's look at the ODE:
$$\frac{du}{dt} = \cos2\pi t$$
- Initial condition:
$$u(0) = 1$$
- The exact solution:
$$u(t) = \frac{1}{2\pi}\sin2\pi t + 1$$
Make a neural network and loss functions like below:
Import Library
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
tf.compat.v1.enable_eager_execution()
import random
import warnings
warnings.filterwarnings("ignore")
Define Network and Hyper-parameters
NN = tf.keras.models.Sequential([
tf.keras.layers.Input((1,)),
tf.keras.layers.Dense(units = 32, activation = 'tanh'),
tf.keras.layers.Dense(units = 32, activation = 'tanh'),
tf.keras.layers.Dense(units = 1)
])
NN.summary()
optm = tf.keras.optimizers.Adam(learning_rate = 0.001)
Define ODE System
- ODE loss:
$$L_{\text{ODE}} = \frac{1}{N_f}\sum_{i=1}^{N_f} \left(\frac{d \text{u}(t_i)}{dt}- \cos2\pi t_i\right)^2$$
- Initial condition loss:
$$L_{IC} = \frac{1}{N_I}\sum_{i=1}^{N_I} \left({\text{u}(0)}- 1\right)^2$$
- Total loss:
$$L_{\text{Total}} = L_{\text{ODE}} + L_{\text{IC}}$$
def ode_system(t, net):
t = t.reshape(-1,1)
t = tf.constant(t, dtype = tf.float32)
t_0 = tf.zeros((1,1))
one = tf.ones((1,1))
with tf.GradientTape() as tape:
tape.watch(t)
u = net(t) # Calculate the ODE loss
u_t = tape.gradient(u, t) # Calculate the Initial Condition loss
ode_loss = u_t - tf.math.cos(2*np.pi*t)
IC_loss = net(t_0) - one
square_loss = tf.square(ode_loss) + tf.square(IC_loss)
total_loss = tf.reduce_mean(square_loss)
return total_loss
Train & Prediction
train_loss_record = []
for itr in range(3000):
train_t = (np.random.rand(20)*2).reshape(-1, 1)
with tf.GradientTape() as tape:
train_loss = ode_system(train_t, NN)
train_loss_record.append(train_loss)
grad_w = tape.gradient(train_loss, NN.trainable_variables)
optm.apply_gradients(zip(grad_w, NN.trainable_variables))
if itr % 1000 == 0:
print('Epoch: {} Loss: {:.4f}'.format(itr, train_loss.numpy()))
plt.figure(figsize = (6, 4))
plt.plot(train_loss_record)
plt.show()
test_t = np.linspace(0, 2, 100)
test_t = tf.constant(test_t, dtype = tf.float32)
pred_u = NN(test_t).numpy()
true_u = np.sin(2*np.pi*test_t)/(2*np.pi) + 1
plt.figure(figsize = (6, 4))
plt.plot(test_t, true_u, 'k', label = 'True', alpha = 0.3)
plt.plot(test_t, pred_u, '--r', label = 'Prediction', linewidth = 3)
plt.legend()
plt.xlabel('t')
plt.ylabel('u')
plt.show()
4. Lab 2: Solve Lab 1 Again using DeepXDE¶
As Physics-informed Neural Networks (PINNs) have gained increasing attention, DeepXDE, a framework specifically designed for PINNs, has been developed by Prof. Lu Lu. Here, we will solve the same example presented above using DeepXDE.
4.1. DeepXDE¶
- Using DeepXDE libray
- DeepXDE is a useful library for solving forward and inverse partial differential equations (PDEs) via physics-informed neural network (PINN)
!pip install deepxde
Change DeepXDE backends
- DeepXDE supports backends for TensorFlow 1.x, TensorFlow 2.x, and PyTorch. We have to change it into TensorFlow 2.x
from deepxde.backend.set_default_backend import set_default_backend
set_default_backend("tensorflow")
4.2. Problem Setup Again¶
- ODE
$$\frac{du}{dt} = \cos2\pi t$$
- Boundary condition:
$$u(0) = 1$$
- The exact solution:
$$u(t) = \frac{1}{2\pi}\sin2\pi t + 1$$
Import Library
import tensorflow as tf
import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
import math as m
Define ODE System
pi = tf.constant(m.pi)
def ode_system(t, u):
du_t = dde.grad.jacobian(u, t)
return du_t - tf.math.cos(2*pi*t)
Define Initial Condition
def boundary(t, on_initial):
return on_initial and np.isclose(t[0], 0)
Define Geometry & Implement Initial Condition
geom = dde.geometry.TimeDomain(0, 2)
ic = dde.IC(geom, lambda t: 1, boundary)
# Reference solution to compute the error
def true_solution(t):
return np.sin(2*np.pi*t)/(2*np.pi) + 1
data = dde.data.PDE(geom,
ode_system,
ic,
num_domain = 30,
num_boundary = 2,
solution = true_solution,
num_test = 100)
Define Network and Hyper-parameters
layer_size = [1] + [32] + [32] + [1]
activation = "tanh"
initializer = "Glorot uniform"
NN = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, NN)
model.compile("adam", lr = 0.001)
Train & Prediction
losshistory, train_state = model.train(epochs = 3000)
dde.saveplot(losshistory, train_state, issave = False, isplot = True)
5. Lab 3: Euler Beam (Solid Mechanics)¶
Having explored the implementation of Physics-informed Neural Networks (PINNs) using both Keras and DeepXDE, we now extend our study to a more practical and familiar problem. Specifically, we will solve the Euler beam bending problem, a fundamental topic in solid mechanics. For this implementation, we will use PyTorch to implement the solution.
5.1. Problem Setup¶
1D problem
- Problem properties (nondimensionalized)
$$E = 1 \operatorname{pa}, \quad I = 1 \operatorname{kg\cdot m^2}, \quad q = 1 \operatorname{N/m}, \quad l = 1 \operatorname{m}$$
- Partial differential equations & boundary conditions
$${\partial^4 y \over \partial x^4} + q = 0, \qquad \text{where} \quad x \in [0,1]$$
- One Dirichlet boundary condition on the left boundary:
$$y(0) = 0$$
- One Neumann boundary condition on the left boundary:
$$y'(0) = 0$$
- Two boundary conditions on the right boundary:
$$y''(1) = 0, \quad y'''(1) = 0$$
- The exact solution is
$$y(x) = -{1 \over 24}x^4 + {1 \over 6}x^3 - {1 \over 4}x^2$$
Make a neural network and loss functions like below:
import numpy as np
import torch
import matplotlib.pyplot as plt
import os
import warnings
warnings.filterwarnings("ignore")
device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
print(device)
Define Collocation Points
x_pts = np.linspace(0, 1, 25)
x0 = x_pts[0] # Left boundary of y = 0
x1 = x_pts[-1] # Right boundary of y = 1
plt.figure(figsize = (6, 3))
plt.scatter(x_pts, np.zeros_like(x_pts), s = 10)
plt.scatter(x0, 0, label = 'x = 0', s = 10, c = 'red')
plt.scatter(x1, 0, label = 'x = 1', s = 10, c = 'red')
plt.xlabel('x')
plt.legend()
plt.title('Collocation Points with Boundaries')
plt.show()
Define Network and Hyper-parameters
import torch.nn as nn
from torch.nn import Linear, Tanh
class PINN(nn.Module):
def __init__(self):
super(PINN, self).__init__()
self.net = nn.Sequential(
Linear(1, 32),
Tanh(),
Linear(32, 32),
Tanh(),
Linear(32, 32),
Tanh(),
Linear(32, 32),
Tanh(),
Linear(32, 1),
Tanh()
)
def forward(self, x):
return self.net(x)
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-3)
loss_fn = nn.MSELoss()
Define Functions to Compute Physics Loss
def derivative(y, t):
return torch.autograd.grad(y, t, grad_outputs = torch.ones_like(y).to(device), create_graph = True)[0]
def requires_grad(x):
return torch.tensor(x, dtype = torch.float32, requires_grad = True).to(device)
Define PDE System
$$\frac{\partial^4 y}{\partial x^4} + q = 0, \qquad \text{where} \quad x \in [0,1]$$
def PDE(model, x):
u = model(x)
u_x = derivative(u, x)
u_xx = derivative(u_x, x)
u_xxx = derivative(u_xx, x)
u_xxxx = derivative(u_xxx, x)
pde = u_xxxx + 1
return pde
Define Boundary Condition
- One Dirichlet boundary condition on the left boundary:
$$y(0) = 0$$
- One Neumann boundary condition on the left boundary:
$$y'(0) = 0$$
- Two boundary conditions on the right boundary:
$$y''(1) = 0, \quad y'''(1) = 0$$
def BC_left(model, x_0):
# x = 0
u_0 = model(x_0)
bc1 = u_0
u_x_0 = derivative(u_0, x_0)
bc2 = u_x_0
return bc1, bc2
def BC_right(model, x_f):
# x = 1
u_f = model(x_f)
u_x_f = derivative(u_f, x_f)
u_xx_f = derivative(u_x_f, x_f)
bc3 = u_xx_f
u_xxx_f = derivative(u_xx_f, x_f)
bc4 = u_xxx_f
return bc3, bc4
Train
x = requires_grad(x_pts.reshape(-1, 1))
x0 = requires_grad(x0.reshape(-1, 1))
x1 = requires_grad(x1.reshape(-1, 1))
for epoch in range(1001):
pde_term = PDE(model, x)
bc1, bc2 = BC_left(model, x0)
bc3, bc4 = BC_right(model, x1)
loss_pde = loss_fn(pde_term, torch.zeros_like(pde_term).to(device))
loss_bc1 = loss_fn(bc1, torch.zeros_like(bc1).to(device))
loss_bc2 = loss_fn(bc2, torch.zeros_like(bc2).to(device))
loss_bc3 = loss_fn(bc3, torch.zeros_like(bc3).to(device))
loss_bc4 = loss_fn(bc4, torch.zeros_like(bc4).to(device))
loss_bc = loss_bc1 + loss_bc2 + loss_bc3 + loss_bc4
loss = loss_pde + loss_bc
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 500 == 0:
print('Epoch: {} Loss: {:.4f} PDELoss: {:.4f} BCLoss: {:.4f}'.format(epoch,
loss.item(),
loss_pde.item(),
loss_bc.item()))
Define Exact Solution
def exact_solution(x):
return -(x ** 4) / 24 + x ** 3 / 6 - x ** 2 / 4
Result
x_ = np.linspace(0, 1, 100)
x_ = torch.tensor(x_.reshape(-1, 1), dtype=torch.float32)
with torch.no_grad():
model.cpu()
model.eval()
u_ = model(x_)
plt.figure(figsize = (6, 4))
plt.plot(x_, exact_solution(x_), c = 'k', label = 'Exact Solution', alpha = 0.5)
plt.plot(x_, u_.detach().numpy(), c = 'r', label = 'Predicted u', linestyle = '--', linewidth = 3)
plt.legend(fontsize = 10)
plt.xticks(fontsize = 10)
plt.yticks(fontsize = 10)
plt.xlabel('x', fontsize = 10)
plt.ylabel('Displacement, y', fontsize = 10)
plt.show()
6. Lab 4: Thin Plate (Solid Mechanics)¶
So far, we have demonstrated the implementation of PINNs using Keras, DeepXDE, and PyTorch. By now, we hope you have gained a solid understanding of how to apply PINNs to solve differential equations. However, all the examples covered thus far have been one-dimensional problems.
Next, we extend our approach to a two-dimensional setting by addressing the thin plate displacement prediction problem. This will illustrate how PINNs adapt to higher-dimensional problems and highlight the additional complexities introduced when solving 2D partial differential equations.
6.1. Problem Setup¶
We will solve thin plate equations to find displacement and stress distribution of thin plate in 2D
- Based on Kirchhoff-Love plate theory, three hypotheses were used
- straight lines normal to the mid-surface remain straight after deformation
- straight lines normal to the mid-surface remain normal to the mid-surface after deformation
- the thickness of the plate does not change during a deformation
A non-uniform stretching force is applied to square elastic plate
Only one quarter of the plate is considered since the geometry and in-plane forces are symmetric (yellow domain)
- Problem properties
$$ E = 50 \operatorname{Mpa}, \quad \nu = 0.3, \quad \omega = 20 \operatorname{mm}, \quad h = 1 \operatorname{mm}, \quad f = 1 \operatorname{Mpa} $$
- Governing equations (Föppl–von Kármán equations) for the isotropic elastic plate:
$$ \begin{align*} &{E \over 1 - \nu^2}\left({\partial^2 u \over \partial x^2} + {1 - \nu \over 2}{\partial^2 u \over \partial y^2} + {1 + \nu \over 2}{\partial^2 v \over \partial x \partial y} \right) = 0\\\\ &{E \over 1 - \nu^2}\left({\partial^2 v \over \partial y^2} + {1 - \nu \over 2}{\partial^2 v \over \partial x^2} + {1 + \nu \over 2}{\partial^2 x \over \partial x \partial y} \right) = 0 \end{align*} $$
- Two Dirichlet boundary conditions at $x = 0,\, y = 0\; (B.C.①, B.C.②)$:
$$ v(x,y) = 0 \qquad \text{at} \quad y = 0\\\\ u(x,y) = 0 \qquad \text{at} \quad x = 0 $$
- Two free boundary conditions at $y = \omega / 2\; (B.C.③)$:
$$ \sigma_{yy} = 0,\quad \sigma_{yx} = 0 $$
- Free boundary condition and in-plane force boundary condition at $x = \omega / 2\; (B.C.④)$:
$$\sigma_{xx} = P \cdot h,\quad \sigma_{xy} = 0$$
Make a neural network and loss funcitons like below:
Numerical Solution (= Exact Solution)
Numerical solution of this problem is illustrated in below figures
- $x, y$ direction displacement $u$, $v$ and
- Stress $\sigma_{xx}$ and $\sigma_{yy}$
Solve this problem using PINN and then compare with a numerical solution
import torch
import matplotlib.pyplot as plt
import numpy as np
import os
import math
device = torch.device('cuda' if torch.cuda.is_available() else 'cpu')
print(device)
# properties
E = 50.0
nu = 0.3
a = 10.0
b = 10.0
h = 1.0
f = 1.0
Define Collocation Points
Nx = 50 # Number of samples
Ny = 50 # Number of samples
x = torch.linspace(0, b, Nx) # Input data for x (Nx x 1)
y = torch.linspace(0, b, Ny) # Input data for y (Ny x 1)
xy = torch.meshgrid(x, y, indexing = 'ij')
Y_domain = torch.cat([xy[0].reshape(-1, 1), xy[1].reshape(-1, 1)], dim = 1)
Y_left = Y_domain[Y_domain[:, 0] == 0]
Y_right = Y_domain[Y_domain[:, 0] == b]
Y_top = Y_domain[Y_domain[:, 1] == b]
Y_bottom = Y_domain[Y_domain[:, 1] == 0]
plt.figure(figsize = (5, 5))
plt.scatter(Y_domain[:, 0], Y_domain[:, 1], s = 4, c = 'gray', alpha = 0.3)
plt.scatter(Y_bottom[:, 0], Y_bottom[:, 1], s = 4, c = 'red')
plt.scatter(Y_top[:, 0], Y_top[:, 1], s = 4, c = 'blue')
plt.scatter(Y_left[:, 0], Y_left[:, 1], s = 4, c = 'green')
plt.scatter(Y_right[:, 0], Y_right[:, 1], s = 4, c = 'purple')
plt.title('Collocation Points')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(True)
plt.show()
Define Network and Hyper-parameters
import torch.nn as nn
from torch.nn import Linear, Tanh
class PINN(nn.Module):
def __init__(self):
super(PINN, self).__init__()
self.net = nn.Sequential(
Linear(2, 64),
Tanh(),
Linear(64, 64),
Tanh(),
Linear(64, 64),
Tanh(),
Linear(64, 64),
Tanh(),
Linear(64, 64),
Tanh(),
Linear(64, 64),
Tanh(),
Linear(64, 2)
)
self._initialize_weights()
def _initialize_weights(self):
for m in self.modules():
if isinstance(m, nn.Linear):
nn.init.xavier_normal_(m.weight)
nn.init.zeros_(m.bias)
def forward(self, x):
x = x.float()
output = self.net(x)
u, v = output[:, 0:1], output[:, 1:2]
return u, v
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-3)
loss_fn = nn.MSELoss()
Define Functions to Compute Physics Loss
def derivative(y, t):
df = torch.autograd.grad(y, t, grad_outputs = torch.ones_like(y).to(device), create_graph = True)[0]
df_x = df[:, 0:1]
df_y = df[:, 1:2]
return df_x, df_y
def requires_grad(x):
return torch.tensor(x, dtype = torch.float32, requires_grad = True).to(device)
Define PDE with Boundary Conditions
Governing equations (Föppl-von Kármán equations) for the isotropic elastic plate:
$$ \begin{align*} &{E \over 1 - \nu^2}\left({\partial^2 u \over \partial x^2} + {1 - \nu \over 2}{\partial^2 u \over \partial y^2} + {1 + \nu \over 2}{\partial^2 v \over \partial x \partial y} \right) = 0\\\\ &{E \over 1 - \nu^2}\left({\partial^2 v \over \partial y^2} + {1 - \nu \over 2}{\partial^2 v \over \partial x^2} + {1 + \nu \over 2}{\partial^2 u \over \partial x \partial y} \right) = 0 \end{align*} $$
def PDE(U, Y_domain):
w, v = model(Y_domain)
dw_x, dw_y = derivative(w, Y_domain)
dw_xx, dw_xy = derivative(dw_x, Y_domain)
_, dw_yy = derivative(dw_y, Y_domain)
dv_x, dv_y = derivative(v, Y_domain)
dv_xx, dv_xy = derivative(dv_x, Y_domain)
_, dv_yy = derivative(dv_y, Y_domain)
force_eq_x = (dw_xx + 0.5 * (1 - nu) * dw_yy + 0.5 * (1 + nu) * dv_xy) * E / (1 - nu**2)
force_eq_y = (dv_yy + 0.5 * (1 - nu) * dv_xx + 0.5 * (1 + nu) * dw_xy) * E / (1 - nu**2)
return force_eq_x.float(), force_eq_y.float()
Two Dirichlet boundary conditions at $x = 0,\, y = 0\; (B.C.① \text{ (left)}, B.C.② \text{ (bottom)})$:
$$ \begin{align*} v(x,y) &= 0 \qquad \text{at} \quad y = 0 \text{: (bottom)} \\ \\ u(x,y) &= 0 \qquad \text{at} \quad x = 0 \text{: (left)} \end{align*} $$
def BC_left(model, Y_left):
# Left edge
w_left, _ = model(Y_left)
return w_left
def BC_bottom(model, Y_bottom):
# Bottom edge
_, v_bottom = model(Y_bottom)
return v_bottom
Two free boundary conditions at $y = \omega / 2\; (B.C.③ \text{ (top)})$:
$$ \sigma_{yy} = 0,\quad \sigma_{yx} = 0 $$
Free boundary condition and in-plane force boundary condition at $x = \omega / 2\; (B.C.④ \text{ (right)})$:
$$ \begin{align*} \sigma_{xx} &= P \cdot h \\ \sigma_{xy} &= 0 \end{align*} $$
Stress
$$ \begin{align*} \sigma_{xx} &= \frac{E}{1 - \nu^2} \left( \frac{\partial u}{\partial x} + \nu \frac{\partial v}{\partial y} \right) \\ \sigma_{yy} &= \frac{E}{1 - \nu^2} \left( \frac{\partial v}{\partial y} + \nu \frac{\partial u}{\partial x} \right)\\ \sigma_{xy} &= \sigma_{yx} = \frac{E}{2(1 + \nu)} \left( \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right) \end{align*} $$
def BC_top(U_top, Y_top):
# Top edge
w, v = model(Y_top)
dw_x, dw_y = derivative(w, Y_top)
dv_x, dv_y = derivative(v, Y_top)
sig_xy = (dv_x + dw_y) * E / (1 + nu) / 2
sig_yy = (dv_y + nu * dw_x) * E / (1 - nu**2)
return sig_xy.float(), sig_yy.float()
def BC_right(U_right, Y_right):
# Right edge
w, v = model(Y_right)
dw_x, dw_y = derivative(w, Y_right)
dv_x, dv_y = derivative(v, Y_right)
sig_xx = (dw_x + nu * dv_y) * E / (1 - nu**2)
sig_xy = (dv_x + dw_y) * E / (1 + nu) / 2
sig_ex = f * h * torch.cos(math.pi / (2 * b) * Y_right[:, 1:2])
sig_xx_ex = sig_xx.float() - sig_ex.float()
return sig_xx_ex.float(), sig_xy.float()
Train
Training a PINN can be computationally expensive and time-consuming. To facilitate experimentation and reproducibility, we provide a trained model along with its corresponding parameters. This allows you to directly evaluate and analyze the results without waiting for extended training times. However, if needed, you can still train the model to explore its behavior further.
## Requires grad
Y_domain = requires_grad(Y_domain)
Y_bottom = requires_grad(Y_bottom)
Y_top = requires_grad(Y_top)
Y_left = requires_grad(Y_left)
Y_right = requires_grad(Y_right)
epochs = 10001
for epoch in range(epochs):
## PDE loss
force_eq_x, force_eq_y = PDE(model, Y_domain)
loss_PDE_x = loss_fn(force_eq_x, torch.zeros_like(force_eq_x))
loss_PDE_y = loss_fn(force_eq_y, torch.zeros_like(force_eq_y))
loss_PDE = loss_PDE_x + loss_PDE_y
## BC loss
bottom_v = BC_bottom(model, Y_bottom)
left_w = BC_left(model, Y_left)
top_sigyy, top_sigyx = BC_top(model, Y_top)
right_sig_xx_ex, right_sigxy = BC_right(model, Y_right)
loss_BC_bottom = loss_fn(bottom_v, torch.zeros_like(bottom_v))
loss_BC_left = loss_fn(left_w, torch.zeros_like(left_w))
loss_BC_top_1 = loss_fn(top_sigyy, torch.zeros_like(top_sigyy))
loss_BC_top_2 = loss_fn(top_sigyx, torch.zeros_like(top_sigyx))
loss_BC_right_1 = loss_fn(right_sig_xx_ex, torch.zeros_like(right_sig_xx_ex))
loss_BC_right_2 = loss_fn(right_sigxy, torch.zeros_like(right_sigxy))
loss_BC_right = loss_BC_right_1 + loss_BC_right_2
loss_BC_top = loss_BC_top_1 + loss_BC_top_2
loss_BC = loss_BC_bottom + loss_BC_left + loss_BC_top + loss_BC_right
loss = loss_PDE + loss_BC
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print('Epoch: {} Loss: {:.4f} PDELoss: {:.4f} BCLoss: {:.4f}'.\
format(epoch, loss.item(), loss_PDE.item(), loss_BC.item()))
Load the Model Parameters
The trained model will be used to avoid the long training time.
from google.colab import drive
drive.mount('/content/drive')
# GPU available
if torch.cuda.is_available():
model = torch.load('/content/drive/MyDrive/DL/DL_data/lab2_adam_model.pt')
Plot Displacement and Stress
def check_stress(model, Y):
w, v = model(Y)
w_x, _ = derivative(w, Y)
_, v_y = derivative(v, Y)
sig_xx = E / (1 - nu ** 2) * (w_x + nu * v_y)
sig_yy = E / (1 - nu ** 2) * (v_y + nu * w_x)
return sig_xx, sig_yy
def plot_results(data, color_legend, titles, w):
plt.figure(figsize = (8, 3))
for idx in range(2):
plt.subplot(1, 2, idx + 1)
plt.scatter(test[:, 0], test[:, 1],
c = data[idx].detach().cpu().numpy(),
cmap = 'rainbow',
s = 5)
plt.clim(color_legend[idx])
plt.title(titles[idx])
plt.axis('square')
plt.xlim((0, b))
plt.ylim((0, b))
plt.colorbar()
plt.tight_layout()
plt.show()
Nx = 100 # Number of samples
Ny = 100 # Number of samples
x_test = torch.linspace(0, b, Nx) # Input data for x (Nx x 1)
y_test = torch.linspace(0, b, Ny) # Input data for y (Ny x 1)
Y_test = torch.meshgrid(x_test, y_test)
test = torch.cat([Y_test[0].reshape(-1, 1), Y_test[1].reshape(-1, 1)], dim = 1)
test_tensor = requires_grad(test)
with torch.no_grad():
model.eval()
result = model(test_tensor)
sigma_xx, sigma_yy = check_stress(model, test_tensor)
sigma = [sigma_xx, sigma_yy]
titles = ['x-displacement ($u$)', 'y-displacement ($v$)',
'$\sigma_{xx}$', '$\sigma_{yy}$']
color_legend = [[0, 0.182], [-0.06, 0.011], [-0.0022, 1.0], [-0.15, 0.45]]
plot_results(result, color_legend[:2], titles[:2], b)
plot_results(sigma, color_legend[2:], titles[2:], b)
We generate two sets of plots - one for the Finite Element Method (FEM) and another for PINN - to quantitatively compare their performance. Additionally, to enable a more qualitative comparison, we plot the solution values along the line $ y = x $.
Plate_data = np.load('/content/drive/MyDrive/DL/DL_data/lab2_Plate_data.npy')
loc = Plate_data[:, 0:2]
w = Plate_data[:, 2:3]
v = Plate_data[:, 3:4]
stress = Plate_data[:, 4:6]
def RESULT(test_model, diag):
diag = requires_grad(diag)
pde_w, pde_v = test_model(diag)
pde_disp = np.hstack([pde_w.cpu().detach().numpy(), pde_v.cpu().detach().numpy()])
sig_x, sig_y = check_stress(test_model, diag)
pde_sig = np.hstack([sig_x.cpu().detach().numpy(), sig_y.cpu().detach().numpy()])
return pde_disp, pde_sig
diag_ = [i for i in range(w.shape[0]) if loc[i, 0] + loc[i, 1] == 10]
diag_x = np.linspace(0, 10, 101).reshape(-1, 1)
diag_y = -diag_x + 10
diag = np.concatenate((diag_x, diag_y), 1)
model_Adam = torch.load('/content/drive/MyDrive/DL/DL_data/lab2_adam_model.pt')
results = {
"FEM": [w[diag_], v[diag_], stress[diag_, 0], stress[diag_, 1]],
"PINN": RESULT(model_Adam, diag),
}
for key in ["PINN"]:
disp, sig = results[key]
results[key] = [disp[:, 0], disp[:, 1], sig[:, 0], sig[:, 1]]
titles = ['x-displacement ($u$)', 'y-displacement ($v$)', '$\sigma_{xx}$', '$\sigma_{yy}$']
line_styles = {'FEM': 'k', 'PINN': '--'}
plt.figure(figsize = (7, 6))
for idx, title in enumerate(titles):
plt.subplot(2, 2, idx + 1)
for label, result in results.items():
plt.plot(diag[:, 0], result[idx], line_styles[label], linewidth = 3, label = label)
plt.xlabel('x')
plt.ylabel(title)
plt.legend()
plt.tight_layout()
plt.show()
6.3. PINN + Data¶
The results for y-displacement and y-direction stress deviate slightly from the ground truth. To improve the prediction performance, we could further optimize the hyperparameters or increase the number of collocation points. However, instead of focusing on these conventional optimization techniques, we will introduce a different approach to enhance accuracy.
Assume that observed data is available, such as displacement values obtained from physical sensors placed on the thin plate or extracted from CAE results. Incorporating such observed data into the Physics-informed Neural Network (PINN) framework can significantly improve predictive accuracy. In the following section, we will demonstrate how to integrate observed data into a PINN and examine how it contributes to enhanced model performance.
Adding data constraints for regularization
$$ \begin{align*} \mathcal{L} &= \omega_{\text{data}} \mathcal{L}_{\text{data}} + \omega_{\text{PDE}} \mathcal{L}_{\text{PDE}} + \omega_{\text{BC}} \mathcal{L}_{\text{BC}}\\\\\\ \mathcal{L}_{\text{data}} &= \frac{1}{N_{\text{data}}} \sum \left(u - \hat u \right)^2 \\ \end{align*} $$
Define Observed Data
x_values = [0, 2, 5, 10]
tr_sample = np.where(np.isin(loc[:, 0], x_values))[0]
x_sample = loc[tr_sample, :]
w_sample = w[tr_sample, :]
v_sample = v[tr_sample, :]
stress_sample = stress[tr_sample, :]
data = np.concatenate([x_sample, w_sample, v_sample, stress_sample], axis = 1)
plt.figure(figsize = (5, 5))
plt.scatter(x_sample[:, 0], x_sample[:, 1], s = 1)
plt.xlim([-0.5, 10.5])
plt.xlabel('x')
plt.ylabel('y')
plt.show()
Define Data-Driven Loss Term
def DATA(model, loc):
w, v = model(loc)
w_x, _ = derivative(w, loc)
_, v_y = derivative(v, loc)
sigma_xx = E / (1 - nu ** 2) * (w_x + nu * v_y)
sigma_yy = E / (1 - nu ** 2) * (v_y + nu * w_x)
sigma = torch.cat((sigma_xx, sigma_yy), dim = 1)
return w, v, sigma
Train
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-3)
loss_fn = nn.MSELoss()
epochs = 10001
loc = requires_grad(data[:, 0:2])
gt_w = requires_grad(data[:, 2:3])
gt_v = requires_grad(data[:, 3:4])
gt_sigma = requires_grad(data[:, 4:6])
for epoch in range(epochs):
# PDE
force_eq_x, force_eq_y = PDE(model, Y_domain)
loss_PDE_x = loss_fn(force_eq_x, torch.zeros_like(force_eq_x))
loss_PDE_y = loss_fn(force_eq_y, torch.zeros_like(force_eq_y))
loss_PDE = loss_PDE_x + loss_PDE_y
# BC
bottom_v = BC_bottom(model, Y_bottom)
left_w = BC_left(model, Y_left)
top_sigyy, top_sigyx = BC_top(model, Y_top)
right_sigxx_ex, right_sigxy = BC_right(model, Y_right)
loss_BC_bottom = loss_fn(bottom_v, torch.zeros_like(bottom_v))
loss_BC_left = loss_fn(left_w, torch.zeros_like(left_w))
loss_BC_top_1 = loss_fn(top_sigyy, torch.zeros_like(top_sigyy))
loss_BC_top_2 = loss_fn(top_sigyx, torch.zeros_like(top_sigyx))
loss_BC_top = loss_BC_top_1 + loss_BC_top_2
loss_BC_right_1 = loss_fn(right_sigxx_ex, torch.zeros_like(right_sigxx_ex))
loss_BC_right_2 = loss_fn(right_sigxy, torch.zeros_like(right_sigxy))
loss_BC_right = loss_BC_right_1 + loss_BC_right_2
loss_BC = loss_BC_bottom + loss_BC_left + loss_BC_top + loss_BC_right
# DATA
data_w, data_v, data_sigma = DATA(model, loc)
loss_Data_w = loss_fn(data_w, gt_w)
loss_Data_v = loss_fn(data_v, gt_v)
loss_Data_sigma = loss_fn(data_sigma, gt_sigma)
loss_Data = loss_Data_w + loss_Data_v + loss_Data_sigma
# Total
loss = loss_PDE + loss_BC + loss_Data
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print('Epoch: {} Loss: {:.4f} PDELoss: {:.4f} BCLoss: {:.4f} DATALoss: {:.4f}'.\
format(epoch, loss.item(), loss_PDE.item(), loss_BC.item(), loss_Data.item()))
with torch.no_grad():
model = torch.load('/content/drive/MyDrive/DL/DL_data/lab2_data_model.pt')
result = model(test_tensor)
sigma_xx, sigma_yy = check_stress(model, test_tensor)
sigma = [sigma_xx, sigma_yy]
plot_results(result, color_legend[:2], titles[:2], w)
plot_results(sigma, color_legend[2:], titles[2:], w)
Comparison with FEM data along the line $y=x$
loc = Plate_data[:, 0:2]
w = Plate_data[:, 2:3]
v = Plate_data[:, 3:4]
stress = Plate_data[:, 4:6]
def RESULT(test_model, diag):
diag = requires_grad(diag)
pde_w, pde_v = test_model(diag)
pde_disp = np.hstack([pde_w.cpu().detach().numpy(), pde_v.cpu().detach().numpy()])
sig_x, sig_y = check_stress(test_model, diag)
pde_sig = np.hstack([sig_x.cpu().detach().numpy(), sig_y.cpu().detach().numpy()])
return pde_disp, pde_sig
diag_ = [i for i in range(w.shape[0]) if loc[i, 0] + loc[i, 1] == 10]
diag_x = np.linspace(0, 10, 101).reshape(-1, 1)
diag_y = -diag_x + 10
diag = np.concatenate((diag_x, diag_y), 1)
model_Adam = torch.load('/content/drive/MyDrive/DL/DL_data/lab2_adam_model.pt')
model_Data = torch.load('/content/drive/MyDrive/DL/DL_data/lab2_data_model.pt')
results = {
"FEM": [w[diag_], v[diag_], stress[diag_, 0], stress[diag_, 1]],
"PINN": RESULT(model_Adam, diag),
"PINN + data": RESULT(model_Data, diag)
}
for key in ["PINN", "PINN + data"]:
disp, sig = results[key]
results[key] = [disp[:, 0], disp[:, 1], sig[:, 0], sig[:, 1]]
titles = ['x-displacement ($u$)', 'y-displacement ($v$)', '$\sigma_{xx}$', '$\sigma_{yy}$']
line_styles = {'FEM': 'k', 'PINN': '--', 'PINN + data': ':'}
plt.figure(figsize = (7, 6))
for idx, title in enumerate(titles):
plt.subplot(2, 2, idx + 1)
for label, result in results.items():
plt.plot(diag[:, 0], result[idx], line_styles[label], linewidth = 3, label = label)
plt.xlabel('x')
plt.ylabel(title)
plt.xlim((0, 10))
plt.legend()
plt.tight_layout()
plt.show()
From this example, we can learn two key lessons:
How to incorporate observational data into the PINN model.
When available, leveraging data enhances the performance of PINNs, making the approach more effective and generally recommended.
7. Inverse Problem: Unknown Parameter Estimation¶
We explore the capabilities of PINNs as a neural network-driven approach for solving forward problems. However, as previously discussed, this method does not demonstrate clear advantages over traditional computational techniques. On the other hand, the use of PINNs for inverse problem-solving presents a promising avenue for advancing scientific computing.
Here, we focus on the inverse problem of not only predicting the flow field around a cylinder but also identifying unknown parameters from observed data. By integrating observational data with governing physical equations, PINNs provide a powerful framework for solving these problems while maintaining consistency with known scientific principles.
7.1. Problem Setup¶
Now, we will solve an inverse problem using PINN
- Unknown properties
$$\rho = ? \quad \mu = ?$$
- Partial differential equations & boundary conditions
$$ \color{red}\rho\left(u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + \frac{1}{\color{red}\rho}\frac{\partial p}{\partial x} \right) - \color{red}\mu\left( \frac{\partial^2u}{\partial^2 x} + \frac{\partial^2u}{\partial^2 y}\right) = 0$$
$$ \color{red}\rho\left(u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + \frac{1}{\color{red}\rho}\frac{\partial p}{\partial y} \right) - \color{red}\mu\left( \frac{\partial^2v}{\partial^2 x} + \frac{\partial^2v}{\partial^2 y}\right) = 0$$
- Inlet boundary condition ($x=-0.5$):
$$u|_{x = -0.5} = 1$$
- Outlet boundary condition ($x=1.5$):
$$v|_{x = 1.5} = 0$$$$p|_{x = 1.5} = 0$$
- Non-slip boundary condition ($y=\pm 0.5$):
$$u|_{y = \pm 0.5} = 0$$$$v|_{y = \pm 0.5} = 0$$
Make a neural network and loss functions like below :
import torch
import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore")
# true values
C1true = 1.0
C2true = 0.01
device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
print(device)
Y_domain= np.load('/content/drive/MyDrive/DL/DL_data/lab3_XY_domain.npy')
gt_domain = np.load('/content/drive/MyDrive/DL/DL_data/lab3_gt_domain.npy')
print(Y_domain.shape)
print(gt_domain.shape)
Define Collocation Points
'Boundary Conditions'
bc_top_x = np.linspace(-0.5, 1.5, 200).reshape(-1, 1)
bc_top_y = 0.5 * np.ones_like(bc_top_x).reshape(-1, 1)
bc_bottom_x = np.linspace(-0.5, 1.5, 200).reshape(-1, 1)
bc_bottom_y = -0.5 * np.ones_like(bc_bottom_x).reshape(-1, 1)
bc_inlet_y = np.linspace(-0.5, 0.5, 100).reshape(-1, 1)
bc_inlet_x = -0.5 * np.ones_like(bc_inlet_y).reshape(-1, 1)
bc_outlet_y = np.linspace(-0.5, 0.5, 100).reshape(-1, 1)
bc_outlet_x = 1.5 * np.ones_like(bc_outlet_y).reshape(-1, 1)
Y_top = np.concatenate((bc_top_x, bc_top_y), 1)
Y_bottom = np.concatenate((bc_bottom_x, bc_bottom_y), 1)
Y_inlet = np.concatenate((bc_inlet_x, bc_inlet_y), 1)
Y_outlet = np.concatenate((bc_outlet_x, bc_outlet_y), 1)
radius = 0.05
theta = np.linspace(0, 2 * np.pi, 200)
bc_cylinder_x = (0 + radius * np.cos(theta)).reshape(-1, 1)
bc_cylinder_y = (0 + radius * np.sin(theta)).reshape(-1, 1)
Y_cylinder = np.concatenate((bc_cylinder_x, bc_cylinder_y), 1)
Y_wall = np.concatenate((Y_top, Y_bottom, Y_cylinder))
plt.scatter(Y_domain[:, 0], Y_domain[:, 1], s = 1, label = 'domain')
plt.scatter(Y_wall[:, 0], Y_wall[:, 1], s = 5, label = 'wall')
plt.scatter(Y_inlet[:, 0], Y_inlet[:, 1], s = 5, label = 'inlet')
plt.scatter(Y_outlet[:, 0], Y_outlet[:, 1], s = 5, label = 'outlet')
plt.axis('scaled')
plt.legend()
plt.show()
Define Observation Points
This CFD data is where the density $\rho$ is 1 and the viscosity $\mu$ is 0.01.
Of course, these two values are assumed to be unknown throughout this example.
Also, assume that 49 sensor locations are randomly selected.
Assume that x-velocity, y-velocity, and pressure are available at these locations.
idx_obs = np.random.choice(len(Y_domain), int(0.0025 * len(Y_domain)), replace = False)
Y_obs = Y_domain[idx_obs]
gt_obs = gt_domain[idx_obs]
print('Y obs: {}'.format(Y_obs.shape))
print('GT obs: {}'.format(gt_obs.shape))
plt.scatter(Y_domain[:, 0], Y_domain[:, 1], c = gt_domain[:, 0], s = 1, cmap = 'jet')
plt.scatter(Y_obs[:, 0], Y_obs[:, 1], marker = 'x', linewidths = 1, s = 20, c = 'k')
plt.axis('scaled')
plt.show()
Define PINN Network
import torch.nn as nn
from torch.nn import Linear, Tanh
class PINN(nn.Module):
def __init__(self):
super(PINN, self).__init__()
self.net = nn.Sequential(
Linear(2, 256),
Tanh(),
Linear(256, 256),
Tanh(),
Linear(256, 256),
Tanh(),
Linear(256, 256),
Tanh(),
Linear(256, 3),
)
self._initialize_weights()
def _initialize_weights(self):
for m in self.modules():
if isinstance(m, nn.Linear):
nn.init.xavier_normal_(m.weight)
nn.init.zeros_(m.bias)
def forward(self, x):
x = x.float()
output = self.net(x)
v1, v2, p = output[:, 0:1], output[:, 1:2], output[:, 2:3]
return v1, v2, p
from IPython.display import clear_output
def PLOT(model, Y_domain, rho_list, vis_list):
pred = model(Y_domain)
plt.figure(figsize = (9, 2))
plt.subplot(1, 3, 1)
plt.title('CFD')
plt.scatter(Y_domain[:, 0].detach().cpu().numpy(),
Y_domain[:, 1].detach().cpu().numpy(),
c = gt_domain[:, 0],
s = 0.1, cmap = 'jet')
plt.axis('scaled')
plt.subplot(1, 3, 2)
plt.title('Prediction')
plt.scatter(Y_domain[:, 0].detach().cpu().numpy(),
Y_domain[:, 1].detach().cpu().numpy(),
c = pred[0].detach().cpu().numpy(),
s = 0.1, cmap = 'jet')
plt.axis('scaled')
plt.subplot(1, 3, 3)
plt.title('Error')
plt.scatter(Y_domain[:, 0].detach().cpu().numpy(),
Y_domain[:, 1].detach().cpu().numpy(),
c = np.abs(pred[0].detach().cpu().numpy() - gt_domain[:, 0:1]),
s = 0.1, cmap = 'jet')
plt.axis('scaled')
plt.tight_layout()
plt.show()
plt.figure(figsize = (6, 2))
plt.subplot(1, 2, 1)
plt.title('Rho Prediction')
plt.plot(rho_list)
plt.subplot(1, 2, 2)
plt.title('Vis Prediction')
plt.plot(vis_list)
plt.tight_layout()
plt.show()
clear_output(wait = True)
Define Functions to Compute Physics Loss
def derivative(y, t):
df = torch.autograd.grad(y, t, grad_outputs = torch.ones_like(y).to(device), create_graph = True)[0]
df_x = df[:, 0:1]
df_y = df[:, 1:2]
return df_x, df_y
def requires_grad(x):
return torch.tensor(x, dtype = torch.float32, requires_grad = True).to(device)
Define PDE
2D Navier-Stokes Equations & boundary conditions (for steady state)
$$ \begin{align*} \rho \left(u{\partial u \over \partial x} + v{\partial u \over \partial y} + {1 \over \rho}{\partial p \over \partial x}\right) - \mu \ \left({\partial^2 u \over {\partial x^2}} + {\partial^2 u \over {\partial y^2}}\right) &= 0\\\\ \rho \left(u{\partial v \over \partial x} + v{\partial v \over \partial y} + {1 \over \rho}{\partial p \over \partial y}\right) - \mu \ \left({\partial^2 v \over {\partial x^2}} + {\partial^2 v \over {\partial y^2}}\right) &= 0\\\\ {\partial u \over \partial x} + {\partial v \over \partial y} &= 0 \end{align*} $$
def PDE(model, Y_domain):
v1, v2, p = model(Y_domain)
dv1_x, dv1_y = derivative(v1, Y_domain)
dv1_xx, _ = derivative(dv1_x, Y_domain)
_, dv1_yy = derivative(dv1_y, Y_domain)
dv2_x, dv2_y = derivative(v2, Y_domain)
dv2_xx, _ = derivative(dv2_x, Y_domain)
_, dv2_yy = derivative(dv2_y, Y_domain)
dp_x, dp_y = derivative(p, Y_domain)
pde_v1 = rho * (v1 * dv1_x + v2 * dv1_y) + dp_x - vis * (dv1_xx + dv1_yy)
pde_v2 = rho * (v1 * dv2_x + v2 * dv2_y) + dp_y - vis * (dv2_xx + dv2_yy)
pde_cont = dv1_x + dv2_y
return pde_v1, pde_v2, pde_cont
Boundary Conditions
- Two Dirichlet boundary conditions on the plate boundary (no-slip condition),
$$u(x,y) = 0, \quad v(x,y) = 0 \qquad \text{at} \quad y = \frac{D}{2} \ \; \text{or} \; -\frac{D}{2}$$
- Two Dirichlet boundary conditions at the cylinder boundary
$$u(\text{cylinder}) = 0, \quad v(\text{cylinder}) = 0$$
- Two Dirichlet boundary conditions at the inlet boundary
$$u(-1,y) = u_{\text{in}}, \quad v(-1,y) = 0$$
- Two Dirichlet boundary conditions at the outlet boundary
$$p(1,y) = 0, \quad v(1,y) = 0$$
def BC_wall(model, Y_wall):
v1_wall, v2_wall, _ = model(Y_wall)
return v1_wall, v2_wall
def BC_inlet(model, Y_inlet):
v1_inlet, v2_inlet, _ = model(Y_inlet)
v1_inlet = v1_inlet - torch.ones_like(v1_inlet).to(device)
return v1_inlet, v2_inlet
def BC_outlet(model, Y_outlet):
_, v2_outlet, p_outlet = model(Y_outlet)
return v2_outlet, p_outlet
Parameters to be Optimized
All we need to do is to define the unknown parameters, specifically $\rho$ (density) and $v$ (viscosity), as trainable variables within the PINN.
# Initialize \rho and \mu as NumPy arrays with the value 1.
rho = torch.ones(1).to(device).requires_grad_(True).float()
vis = torch.ones(1).to(device).requires_grad_(True).float()
model = PINN().to(device)
optimizer = torch.optim.Adam([{'params': model.parameters(), 'lr': 1e-4},
{'params': rho, 'lr': 1e-4},
{'params': vis, 'lr': 1e-4}])
loss_fn = nn.MSELoss()
Train
$$ \begin{align*} \mathcal{L} &= \omega_{\text{data}} \mathcal{L}_{\text{data}} + \omega_{\text{PDE}} \mathcal{L}_{\text{PDE}} + \omega_{\text{BC}} \mathcal{L}_{\text{BC}}\\\\\\ \mathcal{L}_{\text{data}} &= \frac{1}{N_{\text{data}}} \sum \left( \left(u - \hat u \right)^2 + \left(v - \hat v \right)^2 + \left(p - \hat p \right)^2 \right) \\ \end{align*} $$
epoch = 0
best_loss = np.inf
loss_list = []
rho_list, vis_list = [], []
## Requires grad
Y_domain = requires_grad(Y_domain)
Y_wall = requires_grad(Y_wall)
Y_inlet = requires_grad(Y_inlet)
Y_outlet = requires_grad(Y_outlet)
Y_obs = requires_grad(Y_obs)
gt_obs = requires_grad(gt_obs)
while epoch < 50001:
## PDE loss
PDE_v1, PDE_v2, PDE_cont = PDE(model, Y_domain)
loss_PDE_v1 = loss_fn(PDE_v1, torch.zeros_like(PDE_v1))
loss_PDE_v2 = loss_fn(PDE_v2, torch.zeros_like(PDE_v2))
loss_PDE_cont = loss_fn(PDE_cont, torch.zeros_like(PDE_cont))
loss_pde = loss_PDE_v1 + loss_PDE_v2 + loss_PDE_cont
## BC loss
v1_wall, v2_wall = BC_wall(model, Y_wall)
v1_inlet, v2_inlet = BC_inlet(model, Y_inlet)
v2_outlet, p_outlet = BC_outlet(model, Y_outlet)
loss_BC_wall_u = loss_fn(v1_wall, torch.zeros_like(v1_wall))
loss_BC_wall_v = loss_fn(v2_wall, torch.zeros_like(v2_wall))
loss_BC_inlet_u = loss_fn(v1_inlet, torch.zeros_like(v1_inlet))
loss_BC_inlet_v = loss_fn(v2_inlet, torch.zeros_like(v2_inlet))
loss_BC_outlet_v = loss_fn(v2_outlet, torch.zeros_like(v2_outlet))
loss_BC_outlet_p = loss_fn(p_outlet, torch.zeros_like(p_outlet))
loss_bc = loss_BC_wall_u + loss_BC_wall_v + loss_BC_inlet_u + loss_BC_inlet_v + loss_BC_outlet_v + loss_BC_outlet_p
u_obs, v_obs, p_obs = model(Y_obs)
loss_data_u = loss_fn(u_obs, gt_obs[:, 0:1])
loss_data_v = loss_fn(v_obs, gt_obs[:, 1:2])
loss_data_p = loss_fn(p_obs, gt_obs[:, 2:3])
loss_data = loss_data_u + loss_data_v + loss_data_p
loss = loss_pde + loss_bc + 10 * loss_data
optimizer.zero_grad()
loss.backward()
optimizer.step()
loss_list.append(loss.item())
rho_list.append(rho.item())
vis_list.append(vis.item())
if epoch % 1000 == 0:
print('EPOCH : %6d/%6d | Loss_PDE : %5.4f | Loss_BC : %5.4f | Loss_DATA : %5.4f | RHO : %.4f | VIS : %.4f' \
%(epoch, 50000, loss_pde.item(), loss_bc.item(), loss_data.item(), rho.item(), vis.item()))
PLOT(model, Y_domain, rho_list, vis_list)
epoch += 1
rho_list = np.load('/content/drive/MyDrive/DL/DL_data/rho_list.npy', allow_pickle = True).tolist()
vis_list = np.load('/content/drive/MyDrive/DL/DL_data/vis_list.npy', allow_pickle = True).tolist()
with torch.no_grad():
print('Predicted RHO: {:.4f}, Predicted VIS: {:.4f}\n\n'.format(rho_list[-1], vis_list[-1]))
model = torch.load('/content/drive/MyDrive/DL/DL_data/lab3.pt')
PLOT(model, Y_domain, rho_list, vis_list)
8. Inverse Problem: Unknown Boundary Condition Estimation¶
Another common type of inverse problem involves unknown boundary conditions. Here, we will consider a 2D heat transfer problem where one of the boundary conditions is not known. By leveraging PINNs, we can infer the missing boundary condition using observed temperature data while ensuring consistency with the governing heat equation.
8.1 Problem Setup¶
Now, we will predict the unknown boundary condition using PINN in solving the 2D heat transfer equation.
- Problem properties
$$ \kappa = 1 $$
- Partial differential equations & boundary conditions
$$ \kappa \left(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} \right) = 0$$
- Left boundary condition ($x=-1$):
$$T|_{x = -1} = 0.5$$
- Right boundary condition ($x=1$):
$$T|_{x = 1} = 1$$
- Bottom boundary condition ($y = -1$):
$$T|_{y = -1} = 0$$
- The top boundary condition $(y = 1)$ is unknown.
$$T|_{y = 1} = \; ?$$
Make a neural network and loss functions like below:
import torch
import matplotlib.pyplot as plt
import numpy as np
import warnings
warnings.filterwarnings("ignore")
device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
print(device)
gt_domain = np.load('/content/drive/MyDrive/DL/DL_data/lab4_gt_domain.npy')
Define Observation Points
Assume that nine thermometers are positioned at specific locations as shown below. These thermometers provide observed temperature data, which will be used to infer the unknown boundary condition in the 2D heat transfer problem.
plt.figure(figsize = (5, 4))
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
x, y = np.meshgrid(x, y)
idx_xs = np.array([15, 50, 85])
idx_ys = np.array([15, 50, 85])
Y_obs = []
gt_obs = []
for idx_x in idx_xs:
for idx_y in idx_ys:
Y_obs.append(np.concatenate((x[idx_x, idx_y].reshape(-1, 1),
y[idx_x, idx_y].reshape(-1, 1)), 1)[0, :])
gt_obs.append(gt_domain[idx_x, idx_y])
Y_obs = np.array(Y_obs)
gt_obs = np.array(gt_obs).reshape(-1, 1)
plt.pcolormesh(x, y, gt_domain, cmap = 'jet')
plt.colorbar()
plt.scatter(Y_obs[:, 0], Y_obs[:, 1], s = 30, c = 'r', cmap = 'jet')
plt.title('FDM')
plt.xlabel('x')
plt.ylabel('y')
plt.axis("square")
plt.show()
Define Collocation Points
# 'Domain Points'
x = np.linspace(-1, 1, 100)
y = np.linspace(-1, 1, 100)
x, y = np.meshgrid(x, y)
Y_domain = np.concatenate((x.reshape(-1, 1), y.reshape(-1, 1)), 1)
# 'Boundary Conditions'
bc_top_x = np.linspace(-1, 1, 100).reshape(-1, 1)
bc_top_y = np.ones_like(bc_top_x)
bc_bottom_x = np.linspace(-1, 1, 100).reshape(-1, 1)
bc_bottom_y = -1 * np.ones_like(bc_bottom_x)
bc_right_y = np.linspace(-1, 1, 100).reshape(-1, 1)
bc_right_x = np.ones_like(bc_right_y)
bc_left_y = np.linspace(-1, 1, 100).reshape(-1, 1)
bc_left_x = -1 * np.ones_like(bc_left_y)
Y_top = np.concatenate((bc_top_x, bc_top_y), 1)
Y_bottom = np.concatenate((bc_bottom_x, bc_bottom_y), 1)
Y_right = np.concatenate((bc_right_x, bc_right_y), 1)
Y_left = np.concatenate((bc_left_x, bc_left_y), 1)
plt.scatter(Y_domain[:, 0], Y_domain[:, 1], s = 1, alpha = 0.5)
plt.scatter(Y_top[:, 0], Y_top[:, 1], s = 5)
plt.scatter(Y_bottom[:, 0], Y_bottom[:, 1], s = 5)
plt.scatter(Y_right[:, 0], Y_right[:, 1], s = 5)
plt.scatter(Y_left[:, 0], Y_left[:, 1], s = 5)
plt.title('Collocation Points')
plt.axis('scaled')
plt.show()
Define PINN Network
import torch.nn as nn
from torch.nn import Linear, Tanh
class PINN(nn.Module):
def __init__(self):
super(PINN, self).__init__()
self.net = nn.Sequential(
Linear(2, 256),
Tanh(),
Linear(256, 256),
Tanh(),
Linear(256, 256),
Tanh(),
Linear(256, 256),
Tanh(),
Linear(256, 1),
)
self._initialize_weights()
def _initialize_weights(self):
for m in self.modules():
if isinstance(m, nn.Linear):
nn.init.xavier_normal_(m.weight)
nn.init.zeros_(m.bias)
def forward(self, x):
x = x.float()
return self.net(x)
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=1e-4)
loss_fn = nn.MSELoss()
from IPython.display import clear_output
def PLOT(model_test, Y_domain):
pred = model_test(Y_domain)
plt.figure(figsize = (9, 2.3))
plt.subplot(1, 3, 1)
plt.title('CFD')
plt.scatter(Y_domain[:, 0].detach().cpu().numpy(),
Y_domain[:, 1].detach().cpu().numpy(),
c = gt_domain.detach().cpu().numpy().reshape(-1, 1),
vmin = 0, vmax = 1, s = 1, cmap = 'jet')
plt.colorbar()
plt.axis('scaled')
plt.subplot(1, 3, 2)
plt.title('Prediction')
plt.scatter(Y_domain[:, 0].detach().cpu().numpy(),
Y_domain[:, 1].detach().cpu().numpy(),
c = pred.detach().cpu().numpy(),
vmin = 0, vmax = 1, s = 1, cmap = 'jet')
plt.colorbar()
plt.axis('scaled')
plt.subplot(1, 3, 3)
plt.title('Error')
plt.scatter(Y_domain[:, 0].detach().cpu().numpy(),
Y_domain[:, 1].detach().cpu().numpy(),
c = np.abs(pred.detach().cpu().numpy() - np.array(gt_domain.detach().cpu().numpy().reshape(-1, 1))),
vmin = 0, vmax = 0.2, s = 1, cmap = 'jet')
plt.colorbar()
plt.axis('scaled')
plt.tight_layout()
plt.show()
clear_output(wait=True)
Define Functions to Compute Physics Loss
def derivative(y, t):
df = torch.autograd.grad(y, t, grad_outputs = torch.ones_like(y).to(device), create_graph = True)[0]
df_x = df[:, 0:1]
df_y = df[:, 1:2]
return df_x, df_y
def requires_grad(x):
return torch.tensor(x, dtype = torch.float32, requires_grad = True).to(device)
Define Governing Equation and Three Known Boundary Conditions
def PDE(model, Y_domain):
T = model(Y_domain)
dT_x, dT_y = derivative(T, Y_domain)
dT_xx, _ = derivative(dT_x, Y_domain)
_, dT_yy = derivative(dT_y, Y_domain)
pde = 1 * (dT_xx + dT_yy)
return pde
def BC_bottom(model, Y_bottom):
T_bottom = model(Y_bottom)
return T_bottom
def BC_right(model, Y_right):
T_right = model(Y_right)
T_right = T_right - torch.ones_like(T_right).to(device)
return T_right
def BC_left(model, Y_left):
T_left = model(Y_left)
T_left = T_left - 0.5 * torch.ones_like(T_left).to(device)
return T_left
Train
## Requires grad
Y_domain = requires_grad(Y_domain)
Y_obs = requires_grad(Y_obs)
Y_bottom = requires_grad(Y_bottom)
Y_right = requires_grad(Y_right)
Y_left = requires_grad(Y_left)
Y_top = requires_grad(Y_top)
gt_domain = requires_grad(gt_domain)
gt_obs = requires_grad(gt_obs)
best_loss = np.inf
for epoch in range(20001):
## PDE loss
pde_term = PDE(model, Y_domain)
loss_pde = loss_fn(pde_term, torch.zeros_like(pde_term))
## BC loss
T_bottom = BC_bottom(model, Y_bottom)
T_right = BC_right(model, Y_right)
T_left = BC_left(model, Y_left)
loss_BC_bottom = loss_fn(T_bottom, torch.zeros_like(T_bottom))
loss_BC_right = loss_fn(T_right, torch.zeros_like(T_right))
loss_BC_left = loss_fn(T_left, torch.zeros_like(T_left))
loss_bc = loss_BC_bottom + loss_BC_right + loss_BC_left
loss = loss_pde + loss_bc
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print('EPOCH : %6d/%6d| Loss : %5.4f | Loss_PDE : %5.4f | Loss_BC : %5.4f' \
%(epoch, epochs-1, loss.item(), loss_pde.item(), loss_bc.item()))
PLOT(model, Y_domain)
Result
As shown, the prediction deviates significantly from the CFD ground truth due to the absence of one boundary condition. It is important to note that most conventional computational methods fail to produce meaningful results unless all boundary conditions are fully specified.
In contrast, PINNs can still generate a solution by minimizing the loss function, even though the inferred result may not be entirely accurate. This behavior is somewhat analogous to the case of multiple solutions in a system of linear equations, where an underdetermined system can still yield a set of possible solutions that satisfy the given constraints.
Load the model parameters
with torch.no_grad():
model = torch.load('/content/drive/MyDrive/DL/DL_data/lab4_wo_data.pt')
PLOT(model, Y_domain)
8.3. Using Observed Data for Predicting Unknown Boundary Condition¶
We incorporate additional observed data from nine sensors to provide PINN with more information, enhancing its ability to infer the missing boundary condition. This scenario is analogous to finding the best approximate solution in a system of linear equations, where an overdetermined system - having more equations than unknowns - allows for an optimized solution that best satisfies the given constraints.
Train with data
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-4)
loss_fn = nn.MSELoss()
best_loss = np.inf
for epoch in range(20001):
## PDE loss
pde_term = PDE(model, Y_domain)
loss_pde = loss_fn(pde_term, torch.zeros_like(pde_term))
## BC loss
T_bottom = BC_bottom(model, Y_bottom)
T_right = BC_right(model, Y_right)
T_left = BC_left(model, Y_left)
loss_BC_bottom = loss_fn(T_bottom, torch.zeros_like(T_bottom))
loss_BC_right = loss_fn(T_right, torch.zeros_like(T_right))
loss_BC_left = loss_fn(T_left, torch.zeros_like(T_left))
loss_bc = loss_BC_bottom + loss_BC_right + loss_BC_left
## Data loss
T_obs = model(Y_obs)
loss_data = loss_fn(T_obs, gt_obs)
loss = loss_pde + loss_bc + loss_data
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print('EPOCH : %6d/%6d| Loss : %5.4f | Loss_PDE : %5.4f | Loss_BC : %5.4f | Loss_DATA : %5.4f' \
%(epoch, epochs-1, loss.item(), loss_pde.item(), loss_bc.item(), loss_data.item()))
PLOT(model, Y_domain)
Result
The results demonstrate improved predictions, although they still do not perfectly match the ground truth. One approach to enhance accuracy is to place additional sensors to provide more observed data. However, beyond simply increasing the number of sensors, other methods could further improve the solution.
Load the model parameters
with torch.no_grad():
model = torch.load('/content/drive/MyDrive/DL/DL_data/lab4_w_data.pt')
PLOT(model, Y_domain)
8.4 Using Observed Data and Prior Knowledge¶
To make this example more interesting, let's assume that we have prior knowledge that the temperature at the missing boundary is unknown but constant. This additional information can be incorporated into the PINN framework as an additional constraint, improving its ability to infer a more accurate solution.
Train again using both data and prior knowledge
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr = 1e-4)
loss_fn = nn.MSELoss()
for epoch in range(20001):
## PDE loss
pde_term = PDE(model, Y_domain)
loss_pde = loss_fn(pde_term, torch.zeros_like(pde_term))
## BC loss
T_bottom = BC_bottom(model, Y_bottom)
T_right = BC_right(model, Y_right)
T_left = BC_left(model, Y_left)
loss_BC_bottom = loss_fn(T_bottom, torch.zeros_like(T_bottom))
loss_BC_right = loss_fn(T_right, torch.zeros_like(T_right))
loss_BC_left = loss_fn(T_left, torch.zeros_like(T_left))
loss_bc = loss_BC_bottom + loss_BC_right + loss_BC_left
## Data loss
T_obs = model(Y_obs)
loss_data = loss_fn(T_obs, gt_obs)
## Prior Knowledge
T_top = model(Y_top)
mean_value = torch.mean(T_top)
loss_prior = loss_fn(T_top, mean_value)
loss = loss_pde + loss_bc + loss_data + loss_prior
optimizer.zero_grad()
loss.backward()
optimizer.step()
if epoch % 1000 == 0:
print('EPOCH : %6d/%6d| Loss : %5.4f | Loss_PDE : %5.4f | Loss_BC : %5.4f | Loss_DATA : %5.4f | Loss_PRIOR : %5.4f' \
%(epoch, epochs - 1, loss.item(), loss_pde.item(), loss_bc.item(), loss_data.item(), loss_prior.item()))
PLOT(model, Y_domain)
with torch.no_grad():
model = torch.load('/content/drive/MyDrive/DL/DL_data/lab4_w_prior.pt')
PLOT(model, Y_domain)
Result
It is clearly observed that the final prediction is the most accurate among the three different scenarios. The key takeaways from this example are as follows:
Understanding the PINN Framework - We explored how PINNs can integrate physics-based constraints, observed data, and prior knowledge, demonstrating their flexibility in solving inverse problems.
Strategic Problem-Solving in Engineering - When addressing complex engineering challenges, start by assessing the availability of data. Then, consider leveraging additional forms of knowledge — such as physical laws, boundary conditions, or empirical insights — beyond just raw data to enhance predictive accuracy and reliability.
8.5 Comparison of Each Model¶
To evaluate the model's performance more precisely, we plot the predicted temperature along the top boundary. This visualization allows for a direct comparison between the predicted and ground truth values, providing deeper insights into how well the model captures the missing boundary condition.
with torch.no_grad():
model_wo_data = torch.load('/content/drive/MyDrive/DL/DL_data/lab4_wo_data.pt')
model_w_data = torch.load('/content/drive/MyDrive/DL/DL_data/lab4_w_data.pt')
model_w_prior = torch.load('/content/drive/MyDrive/DL/DL_data/lab4_w_prior.pt')
y_top = Y_top[:,0].detach().cpu().numpy()
GT = np.zeros_like(y_top)
T_wo_data = model_wo_data(Y_top)
T_w_data = model_w_data(Y_top)
T_w_prior = model_w_prior(Y_top)
plt.figure(figsize = (6, 4))
plt.plot(y_top, GT, c = 'k', linestyle = '--', linewidth = 2, label = 'GT')
plt.plot(y_top, T_wo_data.detach().cpu().numpy(), c = 'b', linestyle = '--', linewidth = 2, label = 'w/o data')
plt.plot(y_top, T_w_data.detach().cpu().numpy(), c = 'g', linestyle = '--', linewidth = 2, label = 'w/ data')
plt.plot(y_top, T_w_prior.detach().cpu().numpy(), c = 'r', linewidth = 2, label = 'w/ prior')
plt.legend()
plt.grid(True)
plt.show()
%%javascript
$.getScript('https://kmahelona.github.io/ipython_notebook_goodies/ipython_notebook_toc.js')