Proper usage of model.predict for TimePDEs and boundary conditions

[yklcs]

I am trying to solve the Fermi-Pasta-Ulam problem using TimePDE.
The equation I'm trying to solve is defined on a 1D interval x and has a time variable: y(x, t).
I tried using Model.predict to predict results for a custom timeframe. However, the returned value is 1D, so I attempted to reshape it, as I want to be able to generate a y-x-t plot. Also, I am trying to apply the periodic boundary condition y(0) = y(x) and zero boundary condition y(0) = y(x) = 0, but my usage of DirichletBC and PeriodicBC doesn't seem to work.

import numpy as np
import tensorflow as tf
import matplotlib.pyplot as plt
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.collections import PolyCollection
from matplotlib.animation import FuncAnimation
import deepxde as dde

SIGMA = 0.1


def main():
    def pde(x, u):
        du_x = tf.gradients(u, x)[0]
        du_x, du_t = du_x[:, 0:1], du_x[:, 1:2]
        du_xx = tf.gradients(du_x, x)[0][:, 0:1]
        du_xxx = tf.gradients(du_xx, x)[0][:, 0:1]
        # return du_t + 3 * u ** 2 * du_x + du_xxx
        return (
            du_t
            + u * du_x
            + SIGMA ** 2 * du_xxx
            + 2 * SIGMA ** 2 * du_x * du_xx
            + SIGMA ** 2 * u * du_xxx
        )

    def func(x):
        return (
            4 * (1 / np.cosh(2 * (x[:, 0:1] + 8))) ** 2
            + 1 * (1 / np.cosh(x[:, 0:1] + 1)) ** 2
        )

    def boundary_l(x, on_boundary):
        return on_boundary and np.isclose(x[0], -10)

    def boundary_r(x, on_boundary):
        return on_boundary and np.isclose(x[0], 10)

    geom = dde.geometry.Interval(-10, 10)
    timedomain = dde.geometry.TimeDomain(0, 100)
    geomtime = dde.geometry.GeometryXTime(geom, timedomain)

    # bc = dde.DirichletBC(
    #     geomtime, lambda x: np.zeros((len(x), 1)), lambda _, on_boundary: on_boundary
    # )

    bcl = dde.PeriodicBC(geomtime, 0, boundary_l)
    bcr = dde.PeriodicBC(geomtime, 0, boundary_r)

    ic = dde.IC(geomtime, func, lambda _, on_initial: on_initial)

    num_domain, num_boundary, num_initial = 40, 30, 20
    data = dde.data.TimePDE(
        geomtime,
        1,
        pde,
        [bcl, bcr, ic],
        num_domain,
        num_boundary,
        num_initial,
        train_distribution="uniform",
        func=func,
        num_test=1000,
    )

    layer_size = [2] + [32] * 3 + [1]
    activation = "tanh"
    initializer = "Glorot uniform"
    net = dde.maps.FNN(layer_size, activation, initializer)

    model = dde.Model(data, net)

    model.compile("adam", lr=0.001, metrics=["l2 relative error"])
    losshistory, train_state = model.train(epochs=10000)
    # np.savetxt("train_state.dat", train_state)
    # dde.saveplot(losshistory, train_state, issave=True, isplot=True)

    t = geomtime.uniform_points(10000, True)
    res = model.predict(t)[:10000]
    res = res.reshape(250, 40)
    x = np.linspace(-10, 10, 40)
    y = np.linspace(0, 100, 250)

    fig = plt.figure()
    ax = plt.axes(xlim=(-10, 10), ylim=(-5, 5))
    line, = ax.plot([], [], lw=3)

    def init():
        line.set_data([], [])
        return (line,)

    def animate(i):
        print(i)
        line.set_data(x, res[i])
        return (line,)

    anim = FuncAnimation(fig, animate, init_func=init, frames=250, interval=30)

    anim.save("wow.mp4", writer="ffmpeg")

if __name__ == "__main__":
    main()

My questions:

1. Shouldn't the result be 2D, being defined on the interval and having a time variable?
2. How can I properly apply the aforementioned BCs?

Thank you so much for your earlier response and the work on this package!

[lululxvi]

I am a little confused about your problem:

• "periodic boundary condition y(0) = y(x)" should be y(-10, t) = y(10, t)? Your code is right bcl = dde.PeriodicBC(geomtime, 0, boundary_l), and you only need either bcl or bcr.
• "zero boundary condition y(0) = y(x) = 0" should be y(-10, t) = y(10,t) = 0? Your code is right bc = dde.DirichletBC(geomtime, lambda x: np.zeros((len(x), 1)), lambda _, on_boundary: on_boundary)

Do you want to use one BC only or both BCs? If you use zero boundary condition on left and right boundary, then it is automatically periodic BC, right? (The PDE is third order, and thus you need 3 BCs?) Could you show some details why you think the BC doesn't work?

For the prediction, you predict the whole field directly, and then reshape, but this may be hard to debug. An easy way is to predict one timeframe by one timeframe. e.g., to predict t=0

x = np.linspace(-10, 10, 40)
t = np.zeros_like(x)
X = np.vstack((x, t)).T
res = np.ravel(model.predict(X))

You can put this code in your animate by changing t according to i.

[engsbk]

I am a little confused about your problem:

• "periodic boundary condition y(0) = y(x)" should be y(-10, t) = y(10, t)? Your code is right bcl = dde.PeriodicBC(geomtime, 0, boundary_l), and you only need either bcl or bcr.
• "zero boundary condition y(0) = y(x) = 0" should be y(-10, t) = y(10,t) = 0? Your code is right bc = dde.DirichletBC(geomtime, lambda x: np.zeros((len(x), 1)), lambda _, on_boundary: on_boundary)

Do you want to use one BC only or both BCs? If you use zero boundary condition on left and right boundary, then it is automatically periodic BC, right? (The PDE is third order, and thus you need 3 BCs?) Could you show some details why you think the BC doesn't work?

For the prediction, you predict the whole field directly, and then reshape, but this may be hard to debug. An easy way is to predict one timeframe by one timeframe. e.g., to predict t=0

x = np.linspace(-10, 10, 40)
t = np.zeros_like(x)
X = np.vstack((x, t)).T
res = np.ravel(model.predict(X))

You can put this code in your animate by changing t according to i.

This is a great help to display and understand results via animation!
What is the equivalent lines for a 2D domain?
I tried the following code:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import matplotlib
matplotlib.use('Qt5Agg')
%matplotlib inline


x = np.linspace(-1, 1, 256)
y = np.linspace(-1, 1, 256)
#t = np.linspace(0, 0.99, 100)
t = np.zeros_like(x)

X = np.vstack((x,y,t)).T
res = np.ravel(model.predict(X))

#Initially
fig, ax = plt.subplots(figsize=(5, 3))
ax.set_aspect('equal')

xx, yy= np.meshgrid(x,y)
X_star = np.vstack((xx.flatten()[:,None], yy.flatten()[:,None])) 
print ("X_star:", X_star.shape)
print ("res:", res.shape)
res2D = griddata(X, res.flatten(), (xx, yy), method='linear')

cax = ax.pcolormesh(x, y, res2D, vmin=-1, vmax=1, cmap='Blues')
fig.colorbar(cax)

print ("res:", res.shape)

def animate(i):
    t = np.zeros_like(x)+(i/100)
    X = np.vstack((x,y,t)).T
    res = np.ravel(model.predict(X))
    xx, yy = np.meshgrid(x,y)
    X_star = np.vstack((xx.flatten()[:,None], yy.flatten()[:,None]))
    res2D = griddata(X, res.flatten(), (xx, yy), method='linear')
    cax.set_array(res2D)
    ax.set_title('t = ' + str(i/100))
    
anim = FuncAnimation(
    fig, animate, interval=100, frames=len(t)-1)
 
plt.draw()
plt.show()

anim.save('2D_wave' + '.gif', writer='imagemagick')

but I got the following error:

---------------------------------------------------------------------------
QhullError                                Traceback (most recent call last)
<ipython-input-32-279f7796f1ed> in <module>
     23 print ("X_star:", X_star.shape)
     24 print ("res:", res.shape)
---> 25 res2D = griddata(X, res.flatten(), (xx, yy), method='linear')
     26 
     27 cax = ax.pcolormesh(x, y, res2D, vmin=-1, vmax=1, cmap='Blues')

C:\ProgramData\Anaconda3\lib\site-packages\scipy\interpolate\ndgriddata.py in griddata(points, values, xi, method, fill_value, rescale)
    216         return ip(xi)
    217     elif method == 'linear':
--> 218         ip = LinearNDInterpolator(points, values, fill_value=fill_value,
    219                                   rescale=rescale)
    220         return ip(xi)

interpnd.pyx in scipy.interpolate.interpnd.LinearNDInterpolator.__init__()

qhull.pyx in scipy.spatial.qhull.Delaunay.__init__()

qhull.pyx in scipy.spatial.qhull._Qhull.__init__()

QhullError: QH6154 Qhull precision error: Initial simplex is flat (facet 1 is coplanar with the interior point)

While executing:  | qhull d Q12 Qt Qz Qbb Qc
Options selected for Qhull 2019.1.r 2019/06/21:
  run-id 400320134  delaunay  Q12-allow-wide  Qtriangulate  Qz-infinity-point
  Qbbound-last  Qcoplanar-keep  _pre-merge  _zero-centrum  Qinterior-keep
  Pgood  _max-width  2  Error-roundoff 2e-15  _one-merge 1.8e-14
  Visible-distance 1.2e-14  U-max-coplanar 1.2e-14  Width-outside 2.4e-14
  _wide-facet 7.3e-14  _maxoutside 2.4e-14

precision problems (corrected unless 'Q0' or an error)
      3 degenerate hyperplanes recomputed with gaussian elimination
      3 nearly singular or axis-parallel hyperplanes
      3 zero divisors during back substitute
    465 zero divisors during gaussian elimination

The input to qhull appears to be less than 4 dimensional, or a
computation has overflowed.

Qhull could not construct a clearly convex simplex from points:
- p1(v5): -0.99 -0.99     0  0.89
- p256(v4):     0     0     0     1
- p127(v3): -0.0039 -0.0039     0     0
- p255(v2):     1     1     0  0.91
- p0(v1):    -1    -1     0  0.91

The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet.  The maximum round off error for
computing distances is 2e-15.  The center point, facets and distances
to the center point are as follows:

center point  -0.1992  -0.1992        0   0.7426

facet p256 p127 p255 p0 distance=    0
facet p1 p127 p255 p0 distance=    0
facet p1 p256 p255 p0 distance=    0
facet p1 p256 p127 p0 distance=    0
facet p1 p256 p127 p255 distance=    0

These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates.  Trial points
are first selected from points that maximize a coordinate.

The min and max coordinates for each dimension are:
  0:        -1         1  difference=    2
  1:        -1         1  difference=    2
  2:         0         0  difference=    0
  3:         0         1  difference=    1

If the input should be full dimensional, you have several options that
may determine an initial simplex:
  - use 'QJ'  to joggle the input and make it full dimensional
  - use 'QbB' to scale the points to the unit cube
  - use 'QR0' to randomly rotate the input for different maximum points
  - use 'Qs'  to search all points for the initial simplex
  - use 'En'  to specify a maximum roundoff error less than 2e-15.
  - trace execution with 'T3' to see the determinant for each point.

If the input is lower dimensional:
  - use 'QJ' to joggle the input and make it full dimensional
  - use 'Qbk:0Bk:0' to delete coordinate k from the input.  You should
    pick the coordinate with the least range.  The hull will have the
    correct topology.
  - determine the flat containing the points, rotate the points
    into a coordinate plane, and delete the other coordinates.
  - add one or more points to make the input full dimensional.

Please note that 'cubic' interpolation is not applicable for 2D with time as I tried it here.

[lululxvi]

Have no experience on 2d case. If I do this, I would save all the figures first, which is easy, and then use other software to combine figures into video.