Navier Stokes Dynamic surface 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

[@AJAXJR24, @Ryszard2, @katayooneshkofti]

[jdellag]

What happens when you decrease your domain size and/or increase the amount of training points?

[svkarash]

What happens when you decrease your domain size and/or increase the amount of training points?

Thanks for your reply.
In both scenarios, the estimation of eta leading to unexpected results, and the shapes are not correct. fore example the result is

data = dde.data.TimePDE(
    geomtime,
    pde,
    [left_bc_psi,right_bc_psi,left_bc_eta,right_bc_eta,
     bottom_bc_psi_x, bottom_bc_psi_z,
     surface_bc1, surface_bc2, surface_bc3,
     ic_psi, ic_eta],
    num_domain=800,
    num_boundary=400,
    num_initial=400,
    num_test=400
)

model.compile("L-BFGS-B")
losshistory, train_state = model.train(iterations=1000)

Then

num_points_x = 400
num_points_t = 5

x = np.linspace(0, L, num_points_x)  # L is the domain length in x
t = np.linspace(0, T, num_points_t)  # T is the total time
X, T = np.meshgrid(x, t)
XT_flat = np.vstack((X.flatten(), T.flatten())).T
z_value = H  
Z_flat = np.full((XT_flat.shape[0],), z_value)
test_data = np.column_stack((XT_flat, Z_flat))

U_pred = model.predict(test_data)

the results are

psi_pred = U_pred[:, 0]  # First predicted value for each point
eta_pred = U_pred[:, 1]  

and plotted below. Clearly the result is not what I expected and physics predicts

[@engsbk, ]

[jdellag]

Do you have any memory limitations? I only ask because I try to use as many points as I can within my domain and boundary. For example, my resulting plots after training look something like this:

It may also be worth looking into using Anchors to specify specific points in your domain and/or on your boundary so that the PINN is forced to learn the solutions at these specific points.

[AJAXJR24]

Dear @svkarash
Well I am not sure about the things I am going to say but it's worth trying.

1. your Points and iterations seem not enough. Try to run for more iterations and more points. If your loss is ok and your problem converges that may not be the case.
2. In definition of boundary and initial conditions try to use x[:,0:1] for example np.zeros_like(x[:, 0:1]) because the shape of np.zeros_like(x[:, 0]) and aforementioned form are different.
3. review your code to see whether everything has been defined correctly especially in your pde because I assume the surface_bcs you defined in pde is not needed.

[svkarash]

It may also be worth looking into using Anchors to specify specific points in your domain and/or on your boundary so that the PINN is forced to learn the solutions at these specific points

Dear @jdellag
Since I am using the free plan of Google Colab, there are indeed memory limitations. This can affect the scale and complexity of the models I work with. However, in this particular case, I don't believe that memory is the root cause of the issue. The x-z plot should be zero on the left and right boundaries, and it's concerning that they are not.
As for your suggestion about using Anchors, I must admit that I am not familiar with this concept in the context of PINNs. I would greatly appreciate any insights or resources you could share on how to implement Anchors in the code.

[svkarash]

1. Try to run for more iterations and more points. If your loss is

Dear Reza [@AJAXJR24]
Thanks for your suggestions.
I'm planning to implement your suggestions. I think that the predicted eta plot would align with the initial condition at t=0, but unfortunately, the outcome is far from accurate.
Similarly, I had expected the values on the left and right to match the level H, yet this is not the case.

[jdellag]

Heres an example using the Laplace equation on a disk.

Say I wanted to have Anchor points evenly distributed around the disk, I could have a function that generates anchors like:

# Function to generate anchor points. Anchors can be used to improve the performance of the model by ensuring that 
# certain important points in the domain are always included in the training set.
def generate_points(n):
    """
    Generate anchor points in polar coordinates.
    
    Parameters:
    - n (int): Number of points to generate for both radius (r) and angle (theta).
    
    Returns:
    - np.array: Array of [r, theta] pairs.
    """
    # Define radius values ranging from 0 to 1
    r_values = np.linspace(0, 1.0, n)
    # Define angle values ranging from 0 to 2pi
    theta_values = np.linspace(0, 2*np.pi, n)
    points = []
    # Generate cartesian product of r_values and theta_values to get all combinations of [r, theta]
    for r in r_values:
        for theta in theta_values:
            points.append([r, theta])
    return np.array(points)


# Set the number of anchor points to be generated
n = 50
anchors = generate_points(n)

Then use these points when you define your dataset like so:

data = dde.data.PDE(
    geom,
    pde,
    bc,
    num_domain = 9500, 
    num_boundary = 1500,
    anchors = anchors
)

As far as limitations go, I can take care of most of my simulations using 1/10 of an A100 (i.e. under 8GB of memory using all 64bit floats for the ~12500 datapoints used above). While I have an allocation through ACCESS, if I look on RunPod I see that you can use a RTX A4000 for .34/hour.

[svkarash]

anchors = anchors

Dear @jdellag,
Thank you so much for sharing the link and code. I am very excited to read and apply it. hope it help me

[Nimava]

Dear @svkarash ,

I already tried to solve almost the same problem. Based on the help of @forxltk , may be you should define the free surface geometry first and then apply the BC on it. something like this:

fs_condition = tf.less(tf.abs(x[:, 1:2] - zeta), tol)
lmbd = tf.where(fs_condition, tf.ones_like(zeta), tf.zeros_like(zeta))

[svkarash]

Dear @svkarash ,

I already tried to solve almost the same problem. Based on the help of @forxltk , may be you should define the free surface geometry first and then apply the BC on it. something like this:

fs_condition = tf.less(tf.abs(x[:, 1:2] - zeta), tol)
lmbd = tf.where(fs_condition, tf.ones_like(zeta), tf.zeros_like(zeta))

Hi @Nimava,
First of all thanks for yor comment.
If I uderstand your point correctly, you're suggesting that I ought to initially specify the characteristics of the fluid and air, such as defining their presence in the domain, before attempting to solve the equation. This methodology aligns with the approach employed by CFD software like OpenFOAM.

I've implemented this for setting up initial conditions. Nevertheless, I'm inclined to think that this approach may not be suitable for Deepxde, as it seems necessary to define everything from the outset. This includes not only adjusting boundary conditions but also modifying the partial differential equations (PDEs).
I will try to implement your approach which is mentioned above.
I would also appreciate it if you could modify necessary section and comment them.

thanks in advance

[svkarash]

@Nimava,
I'd like to connect with you to discuss certain aspects of the problem and learn more from you.
How can I get in touch with you?

[Nimava]

Dear @svkarash,
Thanks, Here is my email but I am not sure that I can help you or not.
n.vaziri@gmail.com

[svkarash]

n.vaziri@gmail.com

thanks

[svkarash]

Dear @Nimava and experts
I changed the code as @Nimava advised. I only altered the surface equations in the PDE function and kept the surface boundaries the same. The results improved but are still incorrect. The first problem is the incorrect results at the left and right boundaries. Next, the initial wave at t=0 doesn't match the initial condition. Lastly, the wave's peak shifts to the left near x=2 instead of x=5.

initial condition

predicted $\zeta$

Please let me know how can I solve the issues

data = dde.data.TimePDE(
    geomtime,
    pde,
    [left_bc_psi,right_bc_psi,left_bc_eta,right_bc_eta,
     bottom_bc_psi_x, bottom_bc_psi_z,
     surface_bc1, surface_bc2, surface_bc3,
     ic_psi, ic_eta],
    num_domain=5090,
    num_boundary=2048,
    num_initial=2048,
    num_test=5090
)

layers = [3, [50] * 2, [50] * 2, 2]
net = dde.nn.PFNN(layers, "tanh", "Glorot uniform")

[Nimava]

May be these changing can improve the code:

1- You used Drichlet BC for all points of left and right walls but based on my experience in CFD, may be it is better that the uppest points of them (I mean the points that connect to free surface) should be without BC. In this condition, the free surface BCs are applied to these points.

2- I generally used Neumann Bc for these type of problems. You can check it too.

return on_boundary and np.isclose(x,1)

3- What is your initial wave amplitude? It should be very small to have linear waves.