I am trying to solve coupled three PDEs using a DeepXDE PINN model. The Three PDEs are coupled and dependent on each other. I would like to train this model on big time intervals because I would like to make prediction at the end at time = 1E7 (s). Here is my code:
import os
os.environ['DDE_BACKEND'] = 'tensorflow'
import numpy as np
import deepxde as dde
from tensorflow import keras
import tensorflow as tf
import time
import matplotlib.pyplot as plt
# Constants for the first and second PDEs
R = 8.67
h0 = 0.0076
# Constants for the third PDE
Af = 5.8446E5
Ac = 5.8264E8
AlphaP = 0.27
Eaf = 75400
Eac = 1.038E5
MRFTO = 0.05
alpha = 0.1 # 3.1536E6 m^2/year = 0.1 m^2/s
D0 = 1E-11 # 3.1536E-1 m^s/year = 1E-11 m^2/s
def pde_system(X, y, c0):
T, Pb, CAb = y[:, 0:1], y[:, 1:2], y[:, 2:3]
T = tf.nn.softplus(T)
Pb = tf.nn.softplus(Pb)
#CAb = tf.nn.softplus(CAb)
# Using dde.grad.jacobian to calculate first derivatives
dT = dde.grad.jacobian(y, X, i=0)
dPb = dde.grad.jacobian(y, X, i=1)
dCAb = dde.grad.jacobian(y, X, i=2)
dT_x, dT_y, dT_t = dT[:, 0:1], dT[:, 1:2], dT[:, 2:3]
dPb_x, dPb_y, dPb_t = dPb[:, 0:1], dPb[:, 1:2], dPb[:, 2:3]
dCAb_t = dCAb[:, 2:3]
# Using dde.grad.hessian to calculate second derivatives
dT_xx = dde.grad.hessian(y, X, component=0, i=0, j=0)
dT_yy = dde.grad.hessian(y, X, component=0, i=1, j=1)
dPb_xx = dde.grad.hessian(y, X, component=1, i=0, j=0)
dPb_yy = dde.grad.hessian(y, X, component=1, i=1, j=1)
h = h0 * (1 + 0.0215 * (T - 273.15 - 30))
# Calculate kf and kc with the original values
kf = Af * (Pb / 101325) ** AlphaP * tf.exp(-Eaf / R / T)
kc = Ac * (Pb / 101325) ** AlphaP * tf.exp(-Eac / R / T)
return [dT_t - alpha * (dT_xx + dT_yy), dPb_t - (D0 * (dPb_xx + dPb_yy) - (T * R * c0 / h) * dCAb_t), dCAb_t - (MRFTO * kf * tf.exp(-kf * X[:, 2:3]) + kc)]
# Boundary conditions for T
def func_bc_T(x):
top_bottom_bc = 303.15
right_left_bc = 274.15
# Scale the BCs to [0, 1]
top_bottom_bc_scaled = scale(top_bottom_bc, 274.15, 303.15, 0, 1)
right_left_bc_scaled = scale(right_left_bc, 274.15, 303.15, 0, 1)
return np.where((x[:, 1:2] == 0) | (x[:, 1:2] == 1), top_bottom_bc_scaled, right_left_bc_scaled)
# Initial condition for T
def func_ic_T(x):
ic = 295.65
# Scale the IC to [0, 1]
ic_scaled = scale(ic, 274.15, 303.15, 0, 1)
return np.full_like(x[:, 0:1], ic_scaled)
# Boundary conditions for Pb
def func_bc_Pb(x):
top_bottom_bc = 1E7
right_left_bc = 1E6
# Scale the BCs to [0, 1]
top_bottom_bc_scaled = scale(top_bottom_bc, 1E6, 1E7, 0, 1)
right_left_bc_scaled = scale(right_left_bc, 1E6, 1E7, 0, 1)
return np.where((x[:, 1:2] == 0) | (x[:, 1:2] == 1), top_bottom_bc_scaled, right_left_bc_scaled)
# Initial condition for Pb
def func_ic_Pb(x):
ic = 1E7
# Scale the IC to [0, 1]
ic_scaled = scale(ic, 1E6, 1E7, 0, 1)
return np.full_like(x[:, 0:1], ic_scaled)
# Boundary conditions for CAb
def func_bc_CAb(x):
top_bottom_bc = 1.0
right_left_bc = 0.942
top_bottom_bc_scaled = scale(top_bottom_bc, 0.942, 1.0, 0, 1)
right_left_bc_scaled = scale(right_left_bc, 0.942, 1.0, 0, 1)
return np.where((x[:, 1:2] == 0) | (x[:, 1:2] == 1), top_bottom_bc_scaled, right_left_bc_scaled)
# Initial condition for CAb
def func_ic_CAb(x):
ic = 0.942
ic_scaled = scale(ic, 0.942, 1.0, 0, 1)
return np.full_like(x[:, 0:1], ic_scaled)
# Scaling function
def scale(x, old_min, old_max, new_min, new_max):
return ((x - old_min) / (old_max - old_min)) * (new_max - new_min) + new_min
# Inverse scaling function
def inverse_scale(x, old_min, old_max, new_min, new_max):
return ((x - new_min) / (new_max - new_min)) * (old_max - old_min) + old_min
geom = dde.geometry.geometry_2d.Rectangle([0,0], [1e7,1e7])
timedomain = dde.geometry.TimeDomain(0, 1E7)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)
bc_T = dde.DirichletBC(geomtime, func_bc_T, lambda _, on_boundary: on_boundary, component=0)
ic_T = dde.IC(geomtime, func_ic_T, lambda _, on_initial: on_initial, component=0)
bc_Pb = dde.DirichletBC(geomtime, func_bc_Pb, lambda _, on_boundary: on_boundary, component=1)
ic_Pb = dde.IC(geomtime, func_ic_Pb, lambda _, on_initial: on_initial, component=1)
bc_CAb = dde.DirichletBC(geomtime, func_bc_CAb, lambda _, on_boundary: on_boundary, component=2)
ic_CAb = dde.IC(geomtime, func_ic_CAb, lambda _, on_initial: on_initial, component=2)
data = dde.data.TimePDE(geomtime, lambda X, y: pde_system(X, y, c0), [bc_T, ic_T, bc_Pb, ic_Pb, bc_CAb, ic_CAb], num_domain=1000, num_boundary=1000, num_initial=1000, num_test=1000)
layer_size = [3] + [30] * 3 + [3] # Three output nodes for T, Pb, and CAb
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)
model.compile("adam", lr=0.001)
# Initialize an empty list to store the predictions
predictions = []
t0 = time.time()
t1 = time.time()
# Generate a grid of points in the x, y, and t dimensions
x = np.linspace(0, 1, 100)
y = np.linspace(0, 1, 100)
t = np.array([1E7])
X, Y, T = np.meshgrid(x, y, t)
# Reshape the grid of points into a 2D array
grid = np.vstack((X.flatten(), Y.flatten(), T.flatten())).T
# Curriculum learning parameters
initial_c0 = 1 # Start with a lower value
final_c0 = 371 # Final value
num_steps = 1000 # Number of steps for curriculum learning
# Instead of using model.train(epochs=50000), use a loop
for step in range(num_steps): # 200 steps * 1000 epochs/step = 200,000 epochs
# Update c0 value
c0 = initial_c0 + step * (final_c0 - initial_c0) / num_steps
print(f"Step {step + 1}/{num_steps}, c0: {c0}")
# Train for 1000 epochs
losshistory, train_state = model.train(1000)
# Predict the temperature and Pb at each point in the grid
predictions = model.predict(grid)
# Separate the predictions into temperature, Pb, and CAb
temp, Pb, CAb = predictions[:, 0:1], predictions[:, 1:2], predictions[:, 2:3]
# Reshape the predicted temperatures, Pb, and CAb into the shape of the grid
temp = temp.reshape(X.shape)
Pb = Pb.reshape(X.shape)
CAb = CAb.reshape(X.shape)
# Unscale the temperatures, Pb, and CAb
temp_unscaled = inverse_scale(temp, 274.15, 303.15, 0, 1)
Pb_unscaled = inverse_scale(Pb, 1E6, 1E7, 0, 1)
CAb_unscaled = inverse_scale(CAb, 0.942, 1.0, 0, 1)
# Since we only have one time slice, we select that for the heatmaps
temp_slice_unscaled = temp_unscaled[:, :, 0]
Pb_slice_unscaled = Pb_unscaled[:, :, 0]
CAb_slice_unscaled = CAb_unscaled[:, :, 0]
# Plot the heatmap for temperature
plt.figure(figsize=(8, 6))
plt.imshow(temp_slice_unscaled, cmap='turbo', origin='lower', extent=[0, 1, 0, 1])
plt.colorbar(label='Temperature (K)')
plt.title('Temperature Distribution at t = 1')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# Plot the heatmap for Pb
plt.figure(figsize=(8, 6))
plt.imshow(Pb_slice_unscaled, cmap='turbo', origin='lower', extent=[0, 1, 0, 1])
plt.colorbar(label='Pb')
plt.title('Pb Distribution at t = 1')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
# Plot the heatmap for CAb
plt.figure(figsize=(8, 6))
plt.imshow(CAb_slice_unscaled, cmap='turbo', origin='lower', extent=[0, 1, 0, 1])
plt.colorbar(label='CAb')
plt.title('CAb Distribution at t = 1')
plt.xlabel('x')
plt.ylabel('y')
plt.show()
t2 = time.time()
print("Training time: ", t2 - t1)
plt.figure()
plt.semilogy(losshistory.loss_train, label="Train loss")
plt.legend()
plt.show()
save_folder = 'C:\Users\mkhadijeh\Desktop\Codes to test Algorithms' # replace with your directory
plt.savefig(save_folder +'loss_history')
plt.close()
So if I simply made the time for the training between 0 and 1E7, the model will not work! and the predicition results are inaccurate
timedomain = dde.geometry.TimeDomain(0, 1E7)
if I made it:
timedomain = dde.geometry.TimeDomain(0, 1)
the model will work but the prediction will be only accurate for a range of inputs “t” between 0 and 1.
How to deal with this problem?
i.e., what strategies can be applied to to train PINNs with big time intervals
I tried scaling but when scaling the time, it is like I am training the model between 0 and 1 so it did not work!
I am expecting a solution to deal with PINN for big time intervals!