Dirichlet Boundary condition debugging (or visualization)

[KeunwooPark]

Hi.
Is there any tool for visualizing DirichletBC?

I'm trying to solve the Poisson equation with different 2D shapes and it doesn't work as expected.
I want to check if my boundary conditions are set as intended.
If there is a visualization tool for boundary conditions, my problem would be solved, but I couldn't find one.
Do you have any suggestions for debugging DirichletBC?

[lululxvi]

All the training points are in train.dat, and the test points with the network solution are in test.dat. The losses are displayed in screen. You can visualize the results using these files.

[KeunwooPark]

Thank you very much for your answer.

However, train.dat doesn't have enough information that I need.
I want to double-check if a DirichletBC that I set is what I actually intended.
When I plotted train.dat, it seems to contain internal points, too.

But I want to plot only data on the boundaries.
Moreover, I want to plot the boundary values together with the points.

Specifically, I want to make the following boundary condition.

This is what train.dat looks like.

How can I get only boundary points from train.dat?

[lululxvi]

You can check data.num_bcs which is a list of the number points used for each BC. In your case, you only have one BC, and thus the number should be equal to the number you set. For example, the number is 100, and then the first 100 points in train.datare the BC points.

You can also define the three BCs separately.

[KeunwooPark]

Thanks for your help!
Based on your answer, I could write a short code to visualize DirichletBCs.
I'm posting the code for anyone who needs it.

def plot_boundary_conditions(pde):
    boundary_conditions = pde.bcs
    X_train = np.array(pde.train_x)

    plt.figure()
    ax = plt.axes(projection=Axes3D.name)

    mat = []
    for bc in boundary_conditions:
        x = bc.filter(X_train)
        val = bc.func(x)
        m = np.hstack((x, val))
        mat.append(m)

    mat = np.vstack(mat)
    ax.plot3D(mat[:, 0], mat[:, 1], mat[:, 2], ".")
    plt.show()

[engsbk]

def plot_boundary_conditions(pde):
    boundary_conditions = pde.bcs
    X_train = np.array(pde.train_x)

    plt.figure()
    ax = plt.axes(projection=Axes3D.name)

    mat = []
    for bc in boundary_conditions:
        x = bc.filter(X_train)
        val = bc.func(x)
        m = np.hstack((x, val))
        mat.append(m)

    mat = np.vstack(mat)
    ax.plot3D(mat[:, 0], mat[:, 1], mat[:, 2], ".")
    plt.show()

Thanks for the valuable addition!

I tried to use your function but got this error:

----> 3     boundary_conditions = pde.bcs
AttributeError: 'function' object has no attribute 'bcs'

Please assist.

[KeunwooPark]

@engsbk Oh. I chose a very confusing variable name.
The below code should work.

def plot_boundary_conditions(pde_data):
    boundary_conditions = pde_data.bcs
    X_train = np.array(pde_data.train_x)

    plt.figure()
    ax = plt.axes(projection=Axes3D.name)

    mat = []
    for bc in boundary_conditions:
        x = bc.filter(X_train)
        val = bc.func(x)
        m = np.hstack((x, val))
        mat.append(m)

    mat = np.vstack(mat)
    ax.plot3D(mat[:, 0], mat[:, 1], mat[:, 2], ".")
    plt.show()

def main():
    def pde(x, y):
        dy_x = tf.gradients(y, x)[0]
        dy_x, dy_y = dy_x[:, 0:1], dy_x[:, 1:]
        dy_xx = tf.gradients(dy_x, x)[0][:, 0:1]
        dy_yy = tf.gradients(dy_y, x)[0][:, 1:]
        return -dy_xx - dy_yy

    def boundary_outer(x, on_boundary):
        norm = np.sqrt(x[0]**2 + x[1]**2)
        return on_boundary and (norm > 0.9)

    def boundary_inner(x, on_boundary):
        norm = np.sqrt(x[0]**2 + x[1]**2)
        return on_boundary and (norm < 0.9)

    def value_outer(x):
        num_data = x.shape[0]
        return np.zeros((num_data,1))

    def value_inner(x):
        num_data = x.shape[0]
        return np.ones((num_data,1))

    outer = dde.geometry.Rectangle([-1,-1], [1,1])
    inner = dde.geometry.Rectangle([-0.5,-0.5], [0.5, 0.5])
    geom = dde.geometry.CSGDifference(outer, inner)
    bc_outer = dde.DirichletBC(geom, value_outer, boundary_outer)
    bc_inner = dde.DirichletBC(geom, value_inner, boundary_inner)

    data = dde.data.PDE(geom, pde, [bc_outer, bc_inner], num_domain=6000, num_boundary=600, num_test=15000)
    net = dde.maps.FNN([2] + [50] * 4 + [1], "tanh", "Glorot uniform")
    model = dde.Model(data, net)

    plot_boundary_conditions(data)

    model.compile("adam", lr=0.001)
    losshistory, train_state = model.train(epochs=50000)
    model.compile("L-BFGS-B")
    losshistory, train_state = model.train()
    dde.saveplot(losshistory, train_state, issave=True, isplot=True)

if __name__ == "__main__":
    main()

I had to make a clear distinction between the pde function and pde data.

[kmache]

@KeunwooPark Hi,
thank you for your code,
I would like how can I generate my data and save it without train the model.
thanks

[KeunwooPark]

Hi @KarimMache
You can save dde.data.PDE as a file using pickle.

[kmache]

Hi @KeunwooPark,
thank you very,
I try it but get the following error
AttributeError: Can't pickle local object 'main..pde'
please may you share the piece of code to do so?
here is my code.
thanks
def main():
def pde(x, y):
dy_xx = dde.grad.hessian(y, x, i=0, j=0)
dy_yy = dde.grad.hessian(y, x, i=1, j=1)
return -dy_xx - dy_yy - 1

def boundary(_, on_boundary):
    return on_boundary
geom = dde.geometry.Polygon([[0, 0], [1, 0], [1, -1], [-1, -1], [-1, 1], [0, 1]])
bc = dde.DirichletBC(geom, lambda x: 0, boundary)
data = dde.data.PDE(geom, pde, bc, num_domain=1200, num_boundary=120, num_test=1500)
with open('data.pickle', 'wb') as f:
    pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)

if name == "main":
main()

[KeunwooPark]

@KarimMache
I guess you should pickle only the data part of the object.
Try to pickle train_x, train_y, test_x, test_y.
To reuse the data, load them and reassign the data to the PDE object.

[kmache]

Hi @KeunwooPark,
Thank you very much,
I am really new in PINNs. suppose we have 2D space and 1D time pde, I would like to save data with columns train_x, train_y, test_x, test_y, train_t, test_t, or just 'x', 'y', 't'
How can I generate the data and save it, please?

thank you a lot

[kmache]

@KeunwooPark,
may you share the piece of code please?

[KeunwooPark]

@KarimMache
I don't know much about time dependant data. But if you see my example code, it creates XY data by randomly sampling in a domain.

data = dde.data.PDE(geom, pde, [bc_outer, bc_inner], num_domain=6000, num_boundary=600, num_test=15000)
# geom: Domain infomation
# bc_* : boundary conditions
# num_domain: number of samples inside the domain for training
# num_boundary: number of samples on the boundary for training
# num_test: number of samples inside the domain for testing

For the time dependant PDE, Burgers might be helpful.

For the saving part, you can save each data matrix in a PDE object (train_x, train_y, etc.) in a pickle as you did in the previous comment.

[kmache]

@KeunwooPark,
Yeah I try it and It still me error
AttributeError: Can't pickle local object 'main..pde'
here is my code, I just take you code above I put the save part.
def main():
def pde(x, y):
dy_x = tf.gradients(y, x)[0]
dy_x, dy_y = dy_x[:, 0:1], dy_x[:, 1:]
dy_xx = tf.gradients(dy_x, x)[0][:, 0:1]
dy_yy = tf.gradients(dy_y, x)[0][:, 1:]
return -dy_xx - dy_yy

def boundary_outer(x, on_boundary):
    norm = np.sqrt(x[0]**2 + x[1]**2)
    return on_boundary and (norm > 0.9)

def boundary_inner(x, on_boundary):
    norm = np.sqrt(x[0]**2 + x[1]**2)
    return on_boundary and (norm < 0.9)

def value_outer(x):
    num_data = x.shape[0]
    return np.zeros((num_data,1))

def value_inner(x):
    num_data = x.shape[0]
    return np.ones((num_data,1))

outer = dde.geometry.Rectangle([-1,-1], [1,1])
inner = dde.geometry.Rectangle([-0.5,-0.5], [0.5, 0.5])
geom = dde.geometry.CSGDifference(outer, inner)
bc_outer = dde.DirichletBC(geom, value_outer, boundary_outer)
bc_inner = dde.DirichletBC(geom, value_inner, boundary_inner)

data = dde.data.PDE(geom, pde, [bc_outer, bc_inner], num_domain=6000, num_boundary=600, num_test=15000)
with open('data.pickle', 'wb') as f:
    pickle.dump(data, f, pickle.HIGHEST_PROTOCOL)

if name == "main":
main()

[KeunwooPark]

@KarimMache
I meant saving only matrix data of PDE data object.

data = dde.data.PDE(geom, pde, [bc_outer, bc_inner], num_domain=6000, num_boundary=600, num_test=15000)
with open("data.pkl", 'wb') as f:
    pickle.dump(data.train_x, f, pickle.HIGHEST_PROTOCOL)
    pickle.dump(data.train_y, f, pickle.HIGHEST_PROTOCOL)
    pickle.dump(data.test_x, f, pickle.HIGHEST_PROTOCOL)
    pickle.dump(data.test_y, f, pickle.HIGHEST_PROTOCOL)

When loading,

data = dde.data.PDE(geom, pde, [bc_outer, bc_inner])
with open("data.pkl", 'rb') as f:
        data.train_x = pickle.load(f)
        data.train_y = pickle.load(f)
        data.test_x = pickle.load(f)
        data.text_y = pickle.load(f)
        data.num_test = data.test_x.shape[0]