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.