lululxvi / deepxde

A library for scientific machine learning and physics-informed learning
https://deepxde.readthedocs.io
GNU Lesser General Public License v2.1
2.64k stars 739 forks source link

Error from trying to build a 2D diffusion equation #159

Closed Enowyz closed 3 years ago

Enowyz commented 3 years ago

Hi, thanks for this awesome library DeepXDE.

I'm trying to build a 2D diffusion equation based on your 1D sample and some other DeepXDE APIs.

Screenshot_2020-11-10 Diffusion2D_coding_NEW - Jupyter Notebook

My code:

def main():
    def pde(x, y):
        D = tf.Variable(2.0)
        #-----------
        dy_x = tf.gradients(y, x)[0]
        dy_x = dy_x[:, 0:1]
        dy_y = dy_x[:, 1:]
        #-----------
        dy_xx = tf.gradients(dy_x, x)[0][:, 0:1]
        dy_yy = tf.gradients(dy_y, x)[0][:, 1:]  **# a question here: why this part is dy_y gradient descent by x?**
        return D * (dy_xx + dy_yy)

    def func(x):
        return np.zeros([len(x), 1])

    geom = dde.geometry.Rectangle([0,0],[10,10])
    timedomain = dde.geometry.TimeDomain(0, 1)
    geomtime = dde.geometry.GeometryXTime(geom, timedomain)

    def boundary_top(x, on_boundary):
        return on_boundary and np.isclose(x[1], 10) 

    def boundary_bottom(x, on_boundary):
        return on_boundary and np.isclose(x[1], 0) # y coordinate x[1] is 0

    bc_1 = dde.DirichletBC(geomtime, func, boundary_top)     
    bc_2 = dde.DirichletBC(geomtime, func, boundary_bottom)     

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

    data = dde.data.TimePDE(
        geomtime, 
        pde,
        [bc_1, bc_2, ic],
        num_domain=40,
        num_boundary=20,
        num_initial=10,
        #train_distribution='sobol', 
        #anchors=None,
        solution=func,
        num_test=10000,
    )

    layer_size = [2] + [50] * 3 + [1]  # The network input is 2D, the first layer should have two neurons
    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)

    dde.saveplot(losshistory, train_state, issave=True, isplot=True)

My error: ValueError: Cannot feed value of shape (92, 3) for Tensor 'Placeholder_23:0', which has shape '(?, 2)'

I think this might be something wrong with my matrix, I need to reshape it as 92x3? But I don't know where I went wrong.


Updated_1:

I changed my input layer, error is gone. But why instead 2, I need 3 neurons for input layer?


Updated_2:

I tried the solution you wrote in 1D sample.

def func(x):
      return np.sin(np.pi * x[:, 0:1]) * np.exp(-x[:, 1:])

Got an error: RuntimeError: DirichletBC should output 1D values. Use argument 'component' for different components.

What should I do with this component? bc_1 = dde.DirichletBC(geomtime, func, boundary_top, component= ) bc_2 = dde.DirichletBC(geomtime, func, boundary_bottom, component= )


Updated_3: I tried this 'NeumannBC'

bc_1 = dde.NeumannBC(geomtime, func, boundary_top)
bc_2 = dde.NeumannBC(geomtime, func, boundary_bottom)

New shape problem: ValueError: Cannot feed value of shape (92, 2) for Tensor 'Placeholder_74:0', which has shape '(?, 1)'

How does this tensor shape work? I thought it is 2 dimensions, so I have x and y as input, so it should be 2 as input layer.


Updated_4: Should I modify output layer?

lululxvi commented 3 years ago
Enowyz commented 3 years ago

@lululxvi

I changed my code to the way expressing my PDE. 图片


def main():
    def pde(x, y, p):
        D = tf.Variable(2.0)
        #-----------
        dp_x = tf.gradients(p, x)[0][:, 0:1]
        dp_y = tf.gradients(p, y)[0][:, 0:1]
        #-----------
        dp_t = dp_x[:, 1:]
        #-----------
        dp_xx = tf.gradients(dp_x, x)[0][:, 0:1]
        dp_yy = tf.gradients(dp_y, y)[0][:, 0:1] 
        return D * (dp_xx + dp_yy)

    def func(x,y):  
        return ( np.sin(np.pi * x[:, 0:1]) * np.exp(-x[:, 1:]) , 
                np.cos(np.pi * y[:, 0:1]) * np.exp(-x[:, 1:]) )

    geom = dde.geometry.Rectangle([-100,-100],[100,100])
    timedomain = dde.geometry.TimeDomain(0, 1)
    geomtime = dde.geometry.GeometryXTime(geom, timedomain)

    def boundary_top(x, on_boundary):
        return on_boundary and np.isclose(x[1], 100) 

    def boundary_bottom(x, on_boundary):
        return on_boundary and np.isclose(x[1], -100) # y coordinate x[1] is 0

    bc_1 = dde.RobinBC(geomtime, func, boundary_top, component=0)     
    bc_2 = dde.RobinBC(geomtime, func, boundary_bottom, component=0)  

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

    data = dde.data.TimePDE(
        geomtime, 
        pde,
        [bc_1, bc_2, ic],
        num_domain=160,
        num_boundary=20,
        num_initial=10,
        solution=func,
        num_test=1000,
    )

    layer_size = [3] + [50] * 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)

    dde.saveplot(losshistory, train_state, issave=True, isplot=True)

I got an error : TypeError: func() missing 1 required positional argument: 'y'

My question is:

  1. What is this error? How can I fix it?
  2. Could you please explain a little bit more about how you generated func as initial solution using sin and cos?
  3. Is my PDE generated correctly?
  4. For the boundary condition, I see DirichletBC can only deal with 1D parameter. If I want to use it in 2D, I think I should modify component. If so, how should I change it? Because I tried other integers, not working. If not, what should I do for 2D boundary issue then?

Thank you and looking forward to your response.

lululxvi commented 3 years ago

There are two errors:

  1. pde(x, y, p) should be pde(X, p). Here X = (x,y,t) is all inputs, i.e., an array of size Nx3. see the Burgers example and diffusion_1d example
  2. Similarly, func(x,y) should be func(x).
engsbk commented 3 years ago

In 2D Poisson example dy_yy = dde.grad.hessian(y, x, i=1, j=1), but if I'm also interested in time as well, how to define dy_tt? dy_xx = dde.grad.hessian(y, x, i=0, j=0) dy_yy = dde.grad.hessian(y, x, i=1, j=1) dy_tt = ?? (I'm not sure how to define this)

lululxvi commented 3 years ago

dy_tt = dde.grad.hessian(y, x, i=2, j=2)

engsbk commented 3 years ago

Thanks, Lu!

I got this error: i=2 is not valid.

when I'm using: du_tt = dde.grad.hessian(u, x, i=2, j=2)

I suspect the error cane from the meshgrid, so I used a new definition:

def ndmesh(*args):
   args = map(np.asarray,args)
   return np.broadcast_arrays(*[x[(slice(None),)+(None,)*i] for i, x in enumerate(args)])

and called it:

xx, yy, tt = ndmesh(x,y,t)
X = np.vstack((np.ravel(xx), np.ravel(yy), np.ravel(tt))).T
u = exact.flatten()[:, None]

but the error is still occurring..

!Update! It was the first layer of the network. Thanks.

engsbk commented 3 years ago

How can we specify initial conditions and boundary conditions for both x and y in a 2D equation? Can you give an example?

lululxvi commented 3 years ago