Lid-driven cavity flow: Loss function converge really slow

[dang-trinh156]

Dear Prof. Lu Lu and DeepXDE users,

I'm trying to solve the classic lid-driven cavity flow benchmark (https://www.comsol.com/blogs/how-to-solve-a-classic-cfd-benchmark-the-lid-driven-cavity-problem/). It is a stationary problem where x,y is the input of PINN and outputs are u,v (velocity in x and y axis respectively), and p. In the boundary conditions, the top boundary, u = v = 1, and other boundaries are u = v = 0.

As this is a much simpler problem than other Navier Stokes examples in DeepXDE, I directly use the architecture of the Kovasznay flow example ([2] + 4*[50] + [3]) As far as I observe the loss terms, they converge relatively slowly, including the boundary conditions and Navier Stokes equations (momentum and continuity).

While taking a look at Kovasznay flow, all loss terms drop below 1e-4 values:

Given that I reuse the PDEs from the Kovaznay Flow example (https://github.com/lululxvi/deepxde/blob/master/examples/pinn_forward/Kovasznay_flow.py) and double-check my BCs.Have I defined boundary conditions incorrectly or it's just inherently my PINN is too simple to predict the flow?

Here is my code:

import deepxde as dde
import numpy as np
Re = 100
def pde(x,u):
    u_vel, v_vel, p = u[:, 0:1], u[:, 1:2], u[:, 2:]
    u_vel_x = dde.grad.jacobian(u, x, i=0, j=0)
    u_vel_y = dde.grad.jacobian(u, x, i=0, j=1)
    u_vel_xx = dde.grad.hessian(u, x, component=0, i=0, j=0)
    u_vel_yy = dde.grad.hessian(u, x, component=0, i=1, j=1)

    v_vel_x = dde.grad.jacobian(u, x, i=1, j=0)
    v_vel_y = dde.grad.jacobian(u, x, i=1, j=1)
    v_vel_xx = dde.grad.hessian(u, x, component=1, i=0, j=0)
    v_vel_yy = dde.grad.hessian(u, x, component=1, i=1, j=1)

    p_x = dde.grad.jacobian(u, x, i=2, j=0)
    p_y = dde.grad.jacobian(u, x, i=2, j=1)

    momentum_x = (
        u_vel * u_vel_x + v_vel * u_vel_y + p_x - 1 / Re * (u_vel_xx + u_vel_yy)
    )
    momentum_y = (
        u_vel * v_vel_x + v_vel * v_vel_y + p_y - 1 / Re * (v_vel_xx + v_vel_yy)
    )
    continuity = u_vel_x + v_vel_y

    return [momentum_x, momentum_y, continuity]

def BCs_u(x):
    return 0
def top_u(x):
    return 1

def top_BC(x, on_boundary):
    return on_boundary and dde.utils.isclose(x[1], 1)
def left_BC(x,on_boundary):
    return on_boundary and dde.utils.isclose(x[0], 0)
def right_BC(x,on_boundary):
    return on_boundary and dde.utils.isclose(x[0], 1)
def bottom_BC(x,on_boundary):
    return on_boundary and dde.utils.isclose(x[1], 0)
spatial_domain = dde.geometry.Rectangle(xmin=[0., 0.], xmax=[1., 1.])

boundary_condition_u_top = dde.icbc.DirichletBC(
    spatial_domain, top_u, top_BC, component=0)
boundary_condition_v_top = dde.icbc.DirichletBC(
    spatial_domain, top_u, top_BC, component=1)

boundary_condition_u_left = dde.icbc.DirichletBC(
    spatial_domain, BCs_u, left_BC, component=0)
boundary_condition_v_left = dde.icbc.DirichletBC(
    spatial_domain, BCs_u, left_BC, component=1)

boundary_condition_u_right = dde.icbc.DirichletBC(
    spatial_domain, BCs_u, right_BC, component=0)
boundary_condition_v_right = dde.icbc.DirichletBC(
    spatial_domain, BCs_u, right_BC, component=1)

boundary_condition_u_bottom = dde.icbc.DirichletBC(
    spatial_domain, BCs_u, bottom_BC, component=0)
boundary_condition_v_bottom = dde.icbc.DirichletBC(
    spatial_domain, BCs_u, bottom_BC, component=1)
data = dde.data.PDE(
    spatial_domain,
    pde,
    [boundary_condition_u_top, boundary_condition_v_top,
     boundary_condition_u_left,boundary_condition_v_left,
     boundary_condition_u_right,boundary_condition_v_right,
     boundary_condition_u_bottom,boundary_condition_v_bottom],
    num_domain=2601,
    num_boundary=400,
    num_test=100000,
)
net = dde.nn.FNN([2] + 4 * [50] + [3], "tanh", "Glorot normal")
model = dde.Model(data, net)

model.compile("adam", lr=1e-3)
model.train(iterations=30000)

#Visualization of PINN's predicted velocity

X = spatial_domain.uniform_points(10000,True)
Y = model.predict(X,operator = pde)
import matplotlib.pyplot as plt

plt.figure(figsize=(10, 6))
plt.scatter(X[:,0:1], X[:,1:], c=np.sqrt(Y[0]**2 + Y[1]**2), cmap='viridis', s=10)  # 'viridis' is the colormap, you can change it to other colormaps if desired
plt.colorbar(label='Velocity (v)')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('Spatial Domain with Velocity Field')
plt.grid(True)
plt.show()

Here is the visualization between the predicted velocity from the code above and Finite Element results from COMSOL Multiphysics. This shows a complete difference between the predicted velocity and F.E results

PINN's prediction

F.E's prediction (COMSOL)

Thank you for your answer!

[maxerial-benruscon]

Have you found a working architecture which compiles with the benchmark from ghiav et al? I started a thread on a similar issue...

[maxerial-benruscon]

Using LFBGS over adam is recommended for navier stokes, as it better optimizes on the domain