How to choose suitable operator for boundary conditions

[chenr86]

Hello @lululxvi

How does PointSetBC set the boundary condition? What's the limitation of this operator? Can I use PointSetBC if there are some symmetry boundaries?
And what's the difference between PointSetBC and applying exact operators (i.e. Dirichlet, Neumann, ...)for each boundary?

Thank you,
Ru

[lululxvi]

Do you mean OperatorBC or PointSetBC?

[chenr86]

我想知道的是有对称边界的情况下PointSetBC是否还适用。我的理解是如果用了PointSetBC之后，就不用指明具体每个边界是什么类型，可以很方便的设定边界条件和初始条件。那他有什么限制吗，所有的边界条件类型都适用？
I want to know the difference of PointSetBC and other operators, including DirichletBC, NeumannBC, OperatorBC and IC. It seems that once PointSetBC is applied, detailed boundary condition can be ignored(i.e. I don't need to specify which boundary is inlet, wall, outlet, etc. in my code). So I am wondering the applicable condition of PointSetBC.
Besides, if there is a wall with roughness z=0.01m, how can I set the boundary condition?

[chenr86]

Hello @lululxvi

How does PointSetBC set the boundary condition? What's the limitation of this operator? Can I use PointSetBC if there are some symmetry boundaries? And what's the difference between PointSetBC and applying exact operators (i.e. Dirichlet, Neumann, ...)for each boundary?

Thank you, Ru

This problem comes from the example of flow around cylinder. I am now using my CFD result for training data, but the results seems incorrect. Here, PointSetBC is applied for the boundary conditions. All the sample points come from anchors, the num_domain, num_boundary and num_initial are set to be 0.

This is the plot of training point distribution.

For each time step, results using PINNs and from CFD are compared. You can see great difference between them.

[forxltk]

But the Navier-Stokes example of incompressible flow around cylinder is an inverse problem. That is why I ask you the question. For inverse problem, the predicted u should not be difference from the reference solution, since you have some observe u. So I think there must be some mistake in your assignment of data points.

[chenr86]

I show you my code here. The code is modified from the NS inverse problem, now it can solve NS forward problem. I have tried:

• reduce training space from the whole rectangle removing circle to a smaller space after the cylinder
• consider fully developed flow for t = 10~20

But the result does not improve. Any suggestion is welcome.

import deepxde as dde
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
import re
import pandas as pd
import matplotlib as mpl

dde.config.set_random_seed(48)
dde.config.set_default_float('float64')

def load_training_data(num):
    """
    num    : 从读取的数据中抽取的样本点数        
    """
    Time_Interval = 70       # 计算时间取值范围
    InitVel_Interval = 0.1   # 初速度取值范围

    Time_Segment     = 1     # 时间间隔
    InitVel_Segment  = 0.1    # 初速度取值间隔
    data = pd.DataFrame()
    for i in range(1, int(InitVel_Interval/InitVel_Segment+1)):
        for j in range(60, int(Time_Interval/Time_Segment+1)):
            # 按列标题读取指定路径下的数据
            path = str("../data/cylinder-Leo/fac/case1/RESU/v_0")+str(int(i))+str("post/v_0")+str(int(i))+str("-")+str(int(j))+str(".csv")
            feature_name = ['Points_0','Points_1','Points_2','Velocity_0','Velocity_1','Pressure']
            df = pd.read_csv(path, usecols = feature_name)

            # 舍弃不在水平面上的数据；将动压转化成总压；添加时间及初始速度列
            df.drop(df[(df.Points_2 > 0)].index, inplace = True)           
            df = df.drop('Points_2',axis=1)
            df['Pressure'] = df['Pressure'].apply(lambda x: x + 101325)
            df.insert(loc=0, column='Initial_Velocity',value=InitVel_Segment*i)
            df.insert(loc=0, column='Time',value=Time_Segment*j)

            data = pd.concat([data, df])
    data = data.reset_index(drop=True)
    data = np.array(data)
    # data.to_csv("../data/cylinder-Leo/fac/test.csv", sep=',')
    
    # define the training domain: T = [0, 90]
    data1 = data[:, :][data[:, 0] <= 90]               # T_max
    data2 = data1[:, :][data1[:, 2] <= 0.2]           # x_max
    data3 = data2[:, :][data2[:, 2] >= 0.01]           # x_min
    data4 = data3[:, :][data3[:, 3] <= 0.07]           # y_max
    data_domain = data4[:, :][data4[:, 3] >= -0.07]    # y_min
    # choose number of training points: num =7000
    idx = np.random.choice(data_domain.shape[0], num, replace=False)
    t_train = data_domain[idx, 0:1]
    uInit_train = data_domain[idx, 1:2]
    x_train = data_domain[idx, 2:3]
    y_train = data_domain[idx, 3:4]
    p_train = data_domain[idx, 4:5]
    u_train = data_domain[idx, 5:6]
    v_train = data_domain[idx, 6:7]
    # return [t_train, x_train, y_train, p_train, u_train, v_train, uInit_train]
    return [t_train, x_train, y_train, p_train, u_train, v_train]

 # Parameters to set the problem
C1 = 1
C2 = 0.01 

 # Define Navier Stokes Equations (Time-dependent PDEs)
def Navier_Stokes_Equation(x, y):
    u = y[:, 0:1]
    v = y[:, 1:2]
    p = y[:, 2:3]
    du_x = dde.grad.jacobian(y, x, i=0, j=0)
    du_y = dde.grad.jacobian(y, x, i=0, j=1)
    du_t = dde.grad.jacobian(y, x, i=0, j=2)
    dv_x = dde.grad.jacobian(y, x, i=1, j=0)
    dv_y = dde.grad.jacobian(y, x, i=1, j=1)
    dv_t = dde.grad.jacobian(y, x, i=1, j=2)
    dp_x = dde.grad.jacobian(y, x, i=2, j=0)
    dp_y = dde.grad.jacobian(y, x, i=2, j=1)
    du_xx = dde.grad.hessian(y, x, component=0, i=0, j=0)
    du_yy = dde.grad.hessian(y, x, component=0, i=1, j=1)
    dv_xx = dde.grad.hessian(y, x, component=1, i=0, j=0)
    dv_yy = dde.grad.hessian(y, x, component=1, i=1, j=1)
    continuity = du_x + dv_y
    x_momentum = du_t + C1 * (u * du_x + v * du_y) + dp_x - C2 * (du_xx + du_yy)
    y_momentum = dv_t + C1 * (u * dv_x + v * dv_y) + dp_y - C2 * (dv_xx + dv_yy)
    return [continuity, x_momentum, y_momentum]

# Define Spatio-temporal domain
# Rectangular
# x_min, x_max = -0.125, 0.575
# y_min, y_max = -0.205, 0.205
x_min, x_max = 0.01, 0.2
y_min, y_max = -0.07, 0.07
uInit_min, uInit_max = 0, 0.5
Time_Interval = 10
# 计算域外轮廓: X × Y × u= [-0.125, 0.575] × [−0.205, 0.205] 
outline_domain = dde.geometry.Rectangle([x_min, y_min], [x_max, y_max])
# 圆柱区域
center = [0, 0]
radius = 0.01
cylinder_domain = dde.geometry.Disk(center, radius)

# 实际计算域：长方形区域-圆柱区域
space_domain = dde.geometry.CSGDifference(outline_domain, cylinder_domain)

# Time domain: T = [0, 90]
time_domain = dde.geometry.TimeDomain(0, Time_Interval)
# Spatio-temporal domain
geomtime = dde.geometry.GeometryXTime(space_domain, time_domain)

[ob_t, ob_x, ob_y, ob_p, ob_u, ob_v] = load_training_data(num=7000)
ob_xyt = np.hstack((ob_x, ob_y, ob_t))
plt.figure(figsize = (16, 9))
plt.scatter(ob_x, ob_y)
plt.axis("equal")

observe_u = dde.icbc.PointSetBC(ob_xyt, ob_u, component=0)
observe_v = dde.icbc.PointSetBC(ob_xyt, ob_v, component=1)
bcs = [observe_u, observe_v]

# Training datasets and Loss
data = dde.data.TimePDE(
    geomtime,
    Navier_Stokes_Equation,
    bcs,
    num_domain=700,
    num_boundary=200,
    num_initial=100,
    anchors=ob_xyt,
)

# Neural Network setup
layer_size = [3] + [50] * 9 + [3]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.FNN(layer_size, activation, initializer)
model = dde.Model(data, net)

# # callbacks for storing results
# fnamevar = "variables.dat"
# variable = dde.callbacks.VariableValue([C1, C2], period=100, filename=fnamevar)

# Compile, train and save model
model.compile("adam", lr=5e-4, loss_weights = [1, 1, 1, 2, 2])
loss_history, train_state = model.train(epochs=10000, display_every=1000, disregard_previous_best=True
)
dde.saveplot(loss_history, train_state, issave=True, isplot=True)

# model.compile("adam", lr=1e-4)
# loss_history, train_state = model.train(epochs=10000, display_every=1000, disregard_previous_best=True)
model.compile(optimizer = 'L-BFGS-B', loss_weights = [1, 1, 1, 2, 2])
model.train_step.optimizer_kwargs = {'options': {'maxcor': 50, 
                                                   'ftol': 1.0 * np.finfo(float).eps, 
                                                   'maxfun':  50000, 
                                                   'maxiter': 50000, 
                                                   'maxls': 50}}
losshistory, train_state = model.train(display_every = 100)
dde.saveplot(loss_history, train_state, issave=True, isplot=True)
# # model.save(save_path = "./NS_inverse_model/model")
# f = model.predict(ob_input, operator=Navier_Stokes_Equation)
# print("Mean residual:", np.mean(np.absolute(f)))

# Plot the velocity distribution of the flow field:
for t in range(0, 11):
    [ob_t, ob_x, ob_y, ob_p, ob_u, ob_v] = load_training_data(num=5000)
    xyt_pred = np.hstack((ob_x, ob_y, t * np.ones((len(ob_x), 1))))
    uvp_pred = model.predict(xyt_pred)
    x_pred, y_pred, t_pred = xyt_pred[:, 0], xyt_pred[:, 1], xyt_pred[:, 2]
    u_pred, v_pred, p_pred = uvp_pred[:, 0], uvp_pred[:, 1], uvp_pred[:, 2]
    x_true = ob_x[ob_t == t]
    y_true = ob_y[ob_t == t]
    u_true = ob_u[ob_t == t]


    fig, ax = plt.subplots(2, 1, figsize = (16, 9))
    cntr0 = ax[0].tricontourf(x_pred, y_pred, u_pred, levels=np.linspace(-0.1, 0.15, 80), cmap="rainbow")
    cb0 = plt.colorbar(cntr0, ax=ax[0])
    cntr1 = ax[1].tricontourf(x_true, y_true, u_true, levels=np.linspace(-0.1, 0.15, 80), cmap="rainbow")
    cb1 = plt.colorbar(cntr1, ax=ax[1])
    ax[0].set_title("u-PINN " + "(t=" + str(t) + ")", fontsize=16)
    ax[0].axis("scaled")
    ax[0].set_xlabel("X", fontsize=12, family="Arial")
    ax[0].set_ylabel("Y", fontsize=12, family="Arial")

    ax[1].set_title("u-Reference solution " + "(t=" + str(t) + ")", fontsize=16)
    ax[1].axis("scaled")
    ax[1].set_xlabel("X", fontsize=12, family="Arial")
    ax[1].set_ylabel("Y", fontsize=12, family="Arial")
    
    # 标注圆柱所在位置
    cylinder = plt.Circle((0,0), 0.01, color = 'w')
    ax[0].add_patch(mpl.patches.Circle((0,0), 0.01, color='w', fill =False))
    ax[1].add_patch(mpl.patches.Circle((0,0), 0.01, color='w', fill =False))
    
    fig.tight_layout()
    plt.show()

[forxltk]

You mixed the usage of PointSet BC in forward problem and inverse problem. See demos and issues that how to solve a forward problem with suitable boundary condition, for example #80, #634. The latter one, in particular, has to solve a problem basically similar to yours.

[chenr86]

Do you mean PointSetBC can only be used in inverse problem? The examples you gave did not use PointSetBC. I don't know what's the difference between boundary conditions in forward problem and inverse problem.

[forxltk]

No, I mean there is no need for you to use PointSetBC for this problem.

[chenr86]

Now the boundary conditions are defined as follows:

# # Boundary geometry:入口inlet, 出口outlet, 上下对称边界symmetry, 圆柱壁面wall
def boundary_inlet(x,on_boundary):
    return on_boundary and np.isclose(x[0], x_min)

def boundary_outlet(x,on_boundary):
    return on_boundary and np.isclose(x[0], x_max)

def boundary_symm(x,on_boundary):
    return on_boundary and (np.isclose(x[1], y_min) or np.isclose(x[1], y_max))

def boundary_wall(x,on_boundary):
    return on_boundary and (x[0]**2+x[1]**2 < 0.015**2)

def boundary_init(_, on_initial):
    return on_initial

# Boundary values definition
def fun_farfield(x, y, _):
    return dde.grad.jacobian(y, x, i=0, j=1)

def funOne(x):
    return 0.1

def funZero(x):
    return 0.0

# Boundary condition
bc_inlet_u = dde.icbc.DirichletBC(geomtime, funOne, boundary_inlet, component=0)
bc_inlet_v = dde.icbc.DirichletBC(geomtime, funZero, boundary_inlet, component=1)

bc_wall_u = dde.icbc.DirichletBC(geomtime, funZero, boundary_wall, component=0)
bc_wall_v = dde.icbc.DirichletBC(geomtime, funZero, boundary_wall, component=1)

bc_symm_u = dde.icbc.DirichletBC(geomtime, funOne, boundary_symm, component=0)
bc_symm_v = dde.icbc.DirichletBC(geomtime, funZero, boundary_symm, component=1)

bc_outlet_p = dde.DirichletBC(geomtime, funZero, boundary_outlet, component=2)

ic_u = dde.icbc.DirichletBC(geomtime, funOne, boundary_init, component=0)
ic_v = dde.icbc.DirichletBC(geomtime, funZero, boundary_init, component=1)

bcs = [bc_inlet_u, bc_inlet_v, bc_wall_u, bc_wall_v, bc_symm_u, bc_symm_v, bc_outlet_p, ic_u, ic_v]

Still, the result is not reasonable. The velocity of flow decreases before meeting the cylinder. The initial condition may be wrong defined. But I don't know how to modify it. I just copy from other examples. In my code, VP form of NS equation is used, Re is about 1000 with a cylinder of r=0.005, u_inlet=0.1. Training points are not from CFD result any more. Instead, they are randomly generated with more points around cylinder.
Hope someone can help.

[lululxvi]

@chenr86 It seems that you try to solve a forward problem, but also assume knowing the true solution of u, right?

[chenr86]

This is my previous idea. In the first version of code, the true solution is used for training. But for the latest version, the true solution is only used for plotting to compare with PINNs' result.

[lululxvi]

I see. You may start with a simpler setup, e.g., smaller Re, larger circle, etc