Navier Stokes Boundary Conditions

[svkarash]

Dear Dr. @lululxvi
I am dealing with a fluid dynamics problem involving a water tank, surface waves, and a set of partial differential equations (PDEs) along with boundary conditions. I am interested in using DeepXDE to solve this system, but I encountered issues with definition of boundary conditions and initial conditions.

Domain, Boundaries, and Initial Conditions
I defined a water tank with a horizontal axis (x) and a vertical axis (z).
The surface wave is represented by a Gaussian function.
Left, right, and bottom boundaries are solid, with specific conditions for velocity components ($u$ and $w$) and surface elevation ($\zeta$).

Issue
I got incorrect results, and my expectations are not met. Specifically, I expect the plot at $t=0$ to reflect the initial conditions and, at an arbitrary time, for $\zeta$ to approach zero at the left and right boundaries.

Equations

$$\frac{\partial \triangle \psi^{[1]}}{\partial t}-\nu \triangle \triangle\psi^{[1]}=0$$

in the bulk ($0\le z \le H$)

$$ \frac{\partial \psi^{[1]}}{\partial x} = 0 $$

$$ \frac{\partial \psi^{[1]}}{\partial z} = 0 $$

at the bottom $z=0$, and

$$\frac{\partial \zeta^{[1]}}{\partial t}-\frac{\partial \psi^{[1]}}{\partial x}=0;,$$

$$\frac{\partial^2 \psi^{[1]}}{\partial x^2}-\frac{\partial^2 \psi^{[1]}}{\partial z^2}=0 $$

$$ -\frac{\partial^2 \psi^{[1]}}{\partial t\partial z}+3\nu \frac{\partial^2 \psi^{[1]}}{\partial x^2\partial z} +\nu \frac{\partial^2 \psi^{[1]}}{\partial z^3}+g\frac{\partial \zeta^{[1]}}{\partial x}-\frac{\sigma }{\rho }\frac{\partial^3 \zeta^{[1]}}{\partial x^3} =0$$

on the surface $z=H$.
Note that $\psi=\psi(x,z,t)$ while $\zeta=\zeta(x,t)$

Code

L, H, T= 10, 1, 5
nu = 0.0011197
g = 9.81
sigma = 6.35e-2
rho = 1261

geom = dde.geometry.Rectangle([0, 0], [L, H]) 
timedomain = dde.geometry.TimeDomain(0, T)  
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

def pde(X, U):
    # psi:= component=0, zeta:=component=1
    # x:=0, z:=1, t:=2

    # psi derivatives
    psi, zeta = U[:, 0:1], U[:, 1:2]

    psi_x = dde.grad.jacobian(U, X, i=0, j=0)
    psi_z = dde.grad.jacobian(U, X, i=0, j=1)
    psi_t = dde.grad.jacobian(U, X, i=0, j=2)

    psi_xx = dde.grad.hessian(U, X, component=0, i=0, j=0)
    psi_zz = dde.grad.hessian(U, X, component=0, i=1, j=1)
    psi_zt = dde.grad.hessian(U, X, component=0, i=1, j=2)

    psi_xxt = dde.grad.jacobian(psi_xx, X, i=0, j=2)
    psi_zzt = dde.grad.jacobian(psi_zz, X, i=0, j=2)
    psi_xxz = dde.grad.jacobian(psi_xx, X, i=0, j=1)
    psi_zzz = dde.grad.jacobian(psi_zz, X, i=0, j=1)

    psi_xxxx = dde.grad.hessian(psi_xx, X,  i=0, j=0)
    psi_zzxx = dde.grad.hessian(psi_xx, X,  i=0, j=1)
    psi_zzzz = dde.grad.hessian(psi_zz, X,  i=0, j=1)

    # zeta derivaives
    zeta_x = dde.grad.jacobian(U, X, i=1, j=0)
    zeta_z = dde.grad.jacobian(U, X, i=1, j=1)
    zeta_t = dde.grad.jacobian(U, X, i=1, j=2)

    zeta_xx = dde.grad.hessian(U, X,component=1, i=0, j=0)

    zeta_xxx = dde.grad.jacobian(zeta_xx, X, i=0, j=0)

    # PDE from the bulk
    bulk_eq = psi_xxt + psi_zzt - nu * (psi_xxxx+2 * psi_zzxx + psi_zzzz)

    # Bottom boundary conditions (Neumann conditions for psi)
    bottom_bc_psi_x = psi_x  
    bottom_bc_psi_z = psi_z  

    # Surface boundary conditions
    surface_bc1 = zeta_t - psi_x  
    surface_bc2 = psi_xx - psi_zz  
    surface_bc3 = -psi_zt + 2 * nu * psi_xxz + nu * psi_zzz + g * zeta_x - sigma / rho * zeta_xxx  

    return bulk_eq 

def initial_condition_psi(x):
    return np.zeros_like(x[:, 0]) # 0.1*np.exp(-((x[:, 0] - 5) ** 2) / (2 * 0.1584 ** 2))

def initial_condition_zeta(x):
    return 0.1*np.exp(-((x[:, 0] - 5) ** 2) / (0.1584 ** 2))+H

ic_psi = dde.IC(geomtime, initial_condition_psi, lambda _, on_initial: on_initial, component=0)
ic_zeta = dde.IC(geomtime, initial_condition_zeta, lambda _, on_initial: on_initial, component=1)

# Define boundary condition functions
def left_boundary(x, on_boundary):
    return on_boundary and np.isclose(x[0], 0)

def right_boundary(x, on_boundary):
    return on_boundary and np.isclose(x[0], L)

def zero_condition(x):
    return np.zeros_like(x[:, 0])

left_bc_psi  = dde.DirichletBC(geomtime, zero_condition, left_boundary, component=0)
right_bc_psi = dde.DirichletBC(geomtime, zero_condition, right_boundary, component=0)
left_bc_zeta  = dde.DirichletBC(geomtime, zero_condition, left_boundary, component=1)
right_bc_zeta = dde.DirichletBC(geomtime, zero_condition, right_boundary, component=1)

def bottom_boundary(x, on_boundary):
    return on_boundary and np.isclose(x[1], 0)

bottom_bc_psi_x = dde.NeumannBC(geomtime, lambda x: 0, bottom_boundary, component=0)
bottom_bc_psi_z = dde.NeumannBC(geomtime, lambda x: 0, bottom_boundary, component=1)

def surface_boundary(x, on_boundary):
    return on_boundary and np.isclose(x[1], H)

surface_bc1 = dde.OperatorBC(geomtime, lambda x, y, _: dde.grad.jacobian(y, x, i=1, j=2) - dde.grad.jacobian(y, x, i=0, j=0), surface_boundary)

surface_bc2 = dde.OperatorBC(geomtime, lambda x, y, _: dde.grad.hessian(y[:, 0:1], x, i=0, j=0) - dde.grad.hessian(y[:, 0:1], x, i=1, j=1), surface_boundary)

def complex_surface_condition(X, U, _):
    psi_zt  = dde.grad.hessian(U, X, component=0, i=1, j=2)
    psi_xxz = dde.grad.jacobian(dde.grad.hessian(U, X, component=0, i=0, j=0), X, i=0, j=1)
    psi_zzz = dde.grad.jacobian(dde.grad.hessian(U, X, component=0, i=1, j=1), X, i=0, j=1)
    zeta_x   = dde.grad.jacobian(U, X, i=1, j=0)
    zeta_xxx = dde.grad.jacobian(dde.grad.hessian(U, X, component=1, i=0, j=0), X, i=0, j=0)
    return -psi_zt + 3 * nu * psi_xxz + nu * psi_zzz + g * zeta_x - sigma/rho * zeta_xxx

surface_bc3 = dde.OperatorBC(geomtime, complex_surface_condition, surface_boundary)

and the network is defined as below

data = dde.data.TimePDE(
    geomtime,
    pde,
    [left_bc_psi,right_bc_psi,left_bc_zeta,right_bc_zeta,
     bottom_bc_psi_x, bottom_bc_psi_z,
     surface_bc1, surface_bc2, surface_bc3,
     ic_psi, ic_zeta],
    num_domain=200,
    num_boundary=80,
    num_initial=80
)

net = dde.nn.FNN([3] + [50] * 3 + [2], "tanh", "Glorot uniform")
model = dde.Model(data, net)

I'm currently facing challenges with my code, and I'm unsure which part might be causing the incorrect results. I've implemented the domain, boundaries, and initial conditions in DeepXDE, but the outcomes don't align with my expectations.

I'm seeking help to identify and rectify any errors in the code. If possible, guidance on how to modify the code would be greatly appreciated.

Expected Results
What I expect to see for $\zeta$ is plotted below

Thanks in advance