Solving Boussinesq equations

[ruthwik96]

Hello,
I'm trying to solve the 2-D incompressible boussinesq equations.

Domain:

Initial solution plots:

Temperature plot:

But the final results of temperature at t=0 do not match with the initial condition of the temperature though used more num_initial points.

The domain is closed and has a no-slip boundary condition at the walls. I am not sure if i have defined the boundary conditions properly. A reference paper is attached to which has same initial conditions but with the no-penetration boundary conditions along with a zero vorticity boundary condition. But, my case is simple where I have used only no-slip boundary condition.

Loss values:

Reference Paper:
https://math.okstate.edu/people/jiahong/DoeringWuZhaoZhengPhysicaD2018.pdf

Code:

# This code is to solve 2-D Boussinesq Navier-Stokes equations
# It was written by Edom Belayneh, Xiaobai(Chris) Chen, Ruthwik Nadam, and Xiaoming Zheng
# of Central Michigan University, Mount Pleasant, MI 48859
# July 31, 2024

import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
import tensorflow as tf
import os
os.environ["TF_USE_LEGACY_KERAS"]="1

D = 1
L = 1
time = 10
nyu = 1

#Printing values for reference
print("D",D)
print("L",L)
print("time",time)

print("Kinematic Viscocity", nyu)

#Define the PDE
def pde(X, Y):

  u_vel = Y[:,0:1]
  v_vel = Y[:,1:2]
  Temp = Y[:,3:4]

  # Y first column is u, second col is v and third column is p
  # X first column is x, second col is y and third column is t
  # time terms
  du_dt = dde.grad.jacobian(Y,X, i=0, j=2)
  dv_dt = dde.grad.jacobian(Y,X, i=1, j=2)
  dT_dt = dde.grad.jacobian(Y,X, i=3, j=2)

  # space terms
  du_dx = dde.grad.jacobian(Y,X, i=0, j=0)
  du_dy = dde.grad.jacobian(Y,X, i=0, j=1)
  dv_dx = dde.grad.jacobian(Y,X, i=1, j=0)
  dv_dy = dde.grad.jacobian(Y,X, i=1, j=1)
  dp_dx = dde.grad.jacobian(Y,X, i=2, j=0)
  dp_dy = dde.grad.jacobian(Y,X, i=2, j=1)
  dT_dx = dde.grad.jacobian(Y,X, i=3, j=0)
  dT_dy = dde.grad.jacobian(Y,X, i=3, j=1)

  # diffusion terms
  ddu_dxx = dde.grad.hessian(Y,X, i=0, j=0)
  ddu_dyy = dde.grad.hessian(Y,X, i=0, j=1)
  ddv_dxx = dde.grad.hessian(Y,X, i=1, j=0)
  ddv_dyy = dde.grad.hessian(Y,X, i=1, j=1)

  #PDE
  pde_u = du_dt + (u_vel*du_dx) + (v_vel*du_dy) - (nyu * (ddu_dxx + ddu_dyy)) + dp_dx
  pde_v = dv_dt + (u_vel*dv_dx) + (v_vel*dv_dy) - (nyu * (ddv_dxx + ddv_dyy)) + dp_dy - Temp
  pde_T = dT_dt + (u_vel*dT_dx) + (v_vel*dT_dy)
  pde_cont = du_dx + dv_dy

  return [pde_u, pde_v, pde_T, pde_cont]

#Define the Sapce and Time Domain Geometry
geom = dde.geometry.Rectangle(xmin=[-L, -D], xmax=[L, D])
timedomain = dde.geometry.TimeDomain(0,time)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

#No slip wall boundary conditions
bc_wall_u = dde.icbc.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary, component=0)
bc_wall_v = dde.icbc.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary, component=1)

#Initial conditions at t=0
ic_u = dde.icbc.IC(geomtime, lambda x: np.sin(2 * np.pi * x[1]) * np.sin(np.pi * x[0])**2, lambda _, on_initial: on_initial, component=0)
ic_v = dde.icbc.IC(geomtime, lambda x: -np.sin(2 * np.pi * x[0]) * np.sin(np.pi * x[1])**2, lambda _, on_initial: on_initial, component=1)
ic_T = dde.icbc.IC(geomtime, lambda x: np.sin(3 * x[0]) * np.cos(2 * x[1]) + (x[0] - 0.5)**3 + 1/(x[1] + 10), lambda _, on_initial: on_initial, component=3)

train_distribution = "LHS"

data = dde.data.TimePDE(geomtime, pde, [bc_wall_u, bc_wall_v, ic_u, ic_v, ic_T], train_distribution=train_distribution,
                       num_domain= 4000, num_boundary=3000, num_initial=3000, num_test=20000)

print(data.train_x.shape)
print(data.train_x_all.shape)
print(data.train_x_bc.shape)

#Plotting 3D domain
fig = plt.figure(figsize=(12, 6))
ax1 = fig.add_subplot(121, projection='3d')
ax1.scatter(data.train_x[:,0], data.train_x[:,1], data.train_x[:,2],color='blue', label='Domain Points', s=0.5)
ax1.set_title('Points in Spatial Domain')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('t')
ax1.legend()
plt.savefig(f"3D_XYT_plot_{train_distribution}.png")
plt.show()

#Plotting the initial conditions to verify
x = data.bc_points()[6000:15000,0]
y = data.bc_points()[6000:15000,1]
t = data.bc_points()[6000:15000,2]

X, Y = np.meshgrid(x, y)

def ic_u_func(x, y):
    return np.sin(2 * np.pi * y) * (np.sin(np.pi * x))**2

def ic_v_func(x, y):
    return -np.sin(2 * np.pi * x) * (np.sin(np.pi * y))**2

def ic_T_func(x, y):
    return np.sin(3 * x) * np.cos(2 * y) + ((x - 0.5)**3) + (1 / (y + 10))

U = ic_u_func(X, Y)
V = ic_v_func(X, Y)
T = ic_T_func(X, Y)

# Plotting the initial conditions using scatter plots
fig, axs = plt.subplots(1, 3, figsize=(18, 5))

sc1 = axs[0].scatter(X.ravel(), Y.ravel(), c=U.ravel(), s=2, cmap='jet')
fig.colorbar(sc1, ax=axs[0])
axs[0].set_title('Initial Condition $u$')
axs[0].set_xlabel('x')
axs[0].set_ylabel('y')

sc2 = axs[1].scatter(X.ravel(), Y.ravel(), c=V.ravel(), s=2, cmap='jet')
fig.colorbar(sc2, ax=axs[1])
axs[1].set_title('Initial Condition $v$')
axs[1].set_xlabel('x')
axs[1].set_ylabel('y')

sc3 = axs[2].scatter(X.ravel(), Y.ravel(), c=T.ravel(), s=2, cmap='jet')
fig.colorbar(sc3, ax=axs[2])
axs[2].set_title('Initial Condition $T$')
axs[2].set_xlabel('x')
axs[2].set_ylabel('y')

plt.tight_layout()
plt.savefig(f"initial_conditions_{train_distribution}.png")
plt.show()

# Creating the neural network with 3 inputs(x,y,t), 2 hidden layers with 80 neurons and 4 outputs (u,v,p,T).
# Activation function is tanh
net = dde.nn.FNN([3] + [80] * 4 + [4], "tanh", "Glorot normal")

#Creating the model
model = dde.Model(data, net)

#Compile the model with Adam Optimizer
model.compile("adam", lr=1e-3)

# learning rates: 1e-2, 1e-3, 8e-4, 8e-5

#Train the model
losshistory, train_state = model.train(iterations=20000)

#Compile the model again with LBFGS Optimizer after ADAM optimiser is done
dde.optimizers.config.set_LBFGS_options(maxiter=3000)
model.compile("L-BFGS")
losshistory, train_state = model.train()

#Train vs test loss plot
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

#Plotting the temperature plots
x = np.linspace(-L, L, 201)
y = np.linspace(-D, D, 201)
t = np.linspace(0, 1000, 201)
xx, yy, tt = np.meshgrid(x, y, t)

samples = np.vstack([xx.ravel(), yy.ravel(), tt.ravel()]).T
result = model.predict(samples)

# Extract temperature predictions and reshape
Temp_pred = result[:, 3].reshape(201, 201, 201)

# Define the time indices for plotting (0, 2, 4, 6, 8, 10)
time_indices = [0, 40, 80, 120, 160, 200]

# Plot the temperature scatter plots
fig, axes = plt.subplots(2, 3, figsize=(18, 12))
axes = axes.ravel()  # Flatten the axes array for easy iteration

for idx, time_index in enumerate(time_indices):
    ax = axes[idx]
    sc = ax.scatter(xx[:, :, time_index].ravel(), yy[:, :, time_index].ravel(), c=Temp_pred[:, :, time_index].ravel(), s=5, cmap="jet")
    ax.set_title(f't = {t[time_index]:.2f}')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    fig.colorbar(sc, ax=ax, label='Temperature')

plt.tight_layout()
plt.savefig("Temperature_Scatter_Plots.png")
plt.show()

[lululxvi]

• Why your network has 4 outputs? I only see u, v and T in your code.
• Can you try to solve the PDE in the time domain [0, 1]?
• Does the network predict IC correctly?
• Try the hard constraint for BC.

[ruthwik96]

Hello Lu Lu,

Thank you so much for responding and apologies for the late reply as I was busy with my academic semester. Let me answer your questions here,

1. Why your network has 4 outputs? I only see u, v and T in your code.
   The unknown parameters in the PDEs are u,v,p and T. As you can see in the X and Y momentum equations, the solution of u and v are dependent are p. So, we needed the value of p to be in the network and update as the solution of u and v. When we tried providing only 3 outputs, it gave out an error. Typically, we do not give the initial value of the pressure because the Navier-Stokes equations are well-defined with only the initial value of the velocity. In other words, the initial value of the pressure is not needed at all. For a realistic problem where the exact solution is unknown, we can simply give the zero initial value of the pressure, that is, setting the initial value of pressure as 0. But it still gives incorrect results.

2. Can you try to solve the PDE in the time domain [0, 1]?
   Yes, all my initial test runs were for the time domain [0,1]. The loss values are still high for all of them and give incorrect results.

3. Does the network predict IC correctly?
   Unfortunately, it does not. As you can see in the images, the initial condition of T plotting after generating the data points, and the plot at t=0 for T do not match. Still trying to figure out why.

4. Try the hard constraint for BC.
   After looking into your feedback, I did try the hard constraint but it gave out more fancy plots (incorrect). It still does not predict the initial condition well.

Based on my study, I think it probably is because of the high non linearity in the PDEs which troubles the network to predict the initial condition properly. We are now trying to optimize the network by:

• Increasing the network complexity (more neurons or layers)
• Reducing the learning rate and train for more iterations.
• Increasing the number of domain and initial points to improve constraint satisfaction

If an initial condition for pressure is mandatory for the network then, please let me know. And if I missed out any major or minor thing, please correct me as well.

Thank you so much for your time and understanding.