Periodic boundary condition issue

[tangqi]

I think there may be a bug in the periodic boundary condition. See this test code

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import matplotlib.pyplot as plt
import numpy as np
import tensorflow as tf
import sys
import deepxde as dde

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-1

    def boundary(x, on_boundary):
        return on_boundary

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

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

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

    def boundary_t(x, on_boundary):
        return on_boundary and np.isclose(x[0], 1)

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

    geom = dde.geometry.Rectangle([0,0],[1,1])
    bc = dde.DirichletBC(geom, func_bdy, boundary)

    bc1 = dde.DirichletBC(geom, func_bdy, boundary_l)
    bc2 = dde.DirichletBC(geom, func_bdy, boundary_r)
    bc3 = dde.PeriodicBC(geom, 0, boundary_t)
    #bc4 = dde.PeriodicBC(geom, func_bdy, boundary_b)

    data = dde.data.PDE(
        geom, 1, pde, [bc1, bc2, bc3], num_domain=1200, num_boundary=120, num_test=1500
    )
    metrics_=None
    net = dde.maps.FNN([2] + [50] * 4 + [1], "tanh", "Glorot uniform")
    model = dde.Model(data, net)

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

if __name__ == "__main__":
    main()

It is easy to work out that the exact solution is a quadratic function in x. But the results are not.

By looking at a similar issue in https://github.com/lululxvi/deepxde/issues/10, I want to point out the periodic boundary condition is not just u(x, 0)=u(x, 1). In fact, it should be any derivative at both sides are identical. (The formal condition is u(x, y-1)=u(x, y) for any x,y.) There are a few ways to force it. One way has been described in the PINNs paper by Raissi et al.

It is possible that I have a mistake in my code, since I am new to your package. Let me know your suggestion.

Sorry for too many issues, your package has been entertaining us a lot. Thanks a lot for a great package.

Thanks, Qi

[lululxvi]

The exact requirement of periodic BC depends on the problem. For example, let us consider the 1D example u(x) = |sin(\pi x)|, which satisfies u(x-1) = u(x) for any x, so it is periodic. However, the derivatives of u(x) do not exist. Even if we consider right derivative of u(x) at x = 0 and left derivative of u(x) at x = 1, they are not equal. Therefore, from the condition u(x-1) = u(x), we can only have u(0) = u(1). This is why in DeepXDE periodic BC means u(0) = u(1).

For PDE problems, although in many cases, the periodic BC means u and du/dx. But I am not sure whether there is a formal definition. Let me know if you have a reference.

I can allow periodic BC in DeepXDE support u_x(0) = u_x(1), or even higher order derivatives. But I am not sure whether high order derivative is necessary.

[tangqi]

Thanks for the discussion. That absolute value is a fair example. But since derivative is not defined at x=1, it normally will not be a solution to PDE, at least in a strong sense. Anyway, I think u(0)=u(1) is what I thought you have implemented. Thanks for confirming that.

The problem I tried is a well posed problem. But the current implementation will not give the right answer. I think it needs to enforce one more constraint of u’(0)=u’(1).

So thinking about it a bit more, I think one can use your general bc interface to add this constraint. I assume one could enforce two different boundary conditions at the same boundary, as long as the problem is well posed? For instance, for this Poisson case, one could enforce your periodicbc and add another bc of the derivative constraint at the same boundary. Will that be acceptable in the code?

This seems an advantage of your package. I also agree with you it is not a good idea to enforce higher derivatives, as the number of constraints needed depends on pde. BTW, typical periodic means u(x-L)=u(x), certainly it does not say anything about its derivatives as you pointed out.

[lululxvi]

• For your problem, three BCs (bc1, bc2, bc3) may be not enough. Do you have a reference solution?
• Yes, DeepXDE accepts more than one BC at the same boundary.
• But u'(0)=u'(1) is a little tricky to implement. I will implement another periodic BC for this so that users can use it.

[tangqi]

So the problem is essentially just -u_xx = 1 with zero boundary condition on x. So the exact solution should be u=-0.5*x^2+0.5x

BTW, I thought that's how the periodic boundary condition enforced in the previous work of PINNs.

This is not high priority task for us but I am happy we are doing things too wrong. Thanks again for your help!

[tangqi]

OK. One more trouble case with periodic boundary condition. The following code has issue to solve the problem (exact solution is z=t, vz=1). If I impose the Dirichlet boundary, the solver works very well. But if it is periodic, the boundary becomes problematic again. This case is higher priority for us to get it right (than the previous case of Poisson). Any suggestions?

Thanks a lot! PS, if things go as planned, you will hear more from us soon.

def main():
    def pde(x, y): 
        Z, VZ = y[:, 0:1], y[:, 1:2]

        dZ = tf.gradients(Z, x)[0]
        dZ_t = dZ[:, 1:2]

        dVZ = tf.gradients(VZ, x)[0]
        dVZ_theta = dVZ[:, 0:1]
        dVZ_t = dVZ[:, 1:2]

        return [  
            dZ_t - VZ,
            dVZ_t + dVZ_theta 
        ]

    geom = dde.geometry.Interval(0, 2*np.pi) 
    timedomain = dde.geometry.TimeDomain(0, 1) 
    geomtime = dde.geometry.GeometryXTime(geom, timedomain)

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

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

    def fun_bc0(x):
        return x[:, 1:]

    def fun_bc1(x):
        return 1+0*x[:, 1:]

    # testing BCs for each component
    bc1 = dde.PeriodicBC(geomtime, 0, boundary_r)
    bc2 = dde.PeriodicBC(geomtime, 0, boundary_l)
    bc3 = dde.PeriodicBC(geomtime, 1, boundary_r)
    bc4 = dde.PeriodicBC(geomtime, 1, boundary_l)

    bc_a = dde.DirichletBC(
        geomtime, fun_bc0, lambda _, on_boundary: on_boundary, component=0
    )
    bc_b = dde.DirichletBC(
        geomtime, fun_bc1, lambda _, on_boundary: on_boundary, component=1
    )

    def boundary(_, on_initial):
        return on_initial

    ic1 = dde.IC(geomtime, lambda X: np.zeros(X[:,0:1].shape), boundary, component=0) # z0(theta)
    ic2 = dde.IC(geomtime, lambda X: np.ones(X[:,0:1].shape), boundary, component=1) # vz0(theta)

    data = dde.data.TimePDE(
    #geomtime, 2, pde, [bc_a, bc_b, ic1,ic2], num_domain = 1000, num_boundary = 100, num_initial = 100
    geomtime, 2, pde, [bc1,bc2,bc3,bc4, ic1,ic2], num_domain = 1000, num_boundary = 100, num_initial = 100
    )

    net = dde.maps.FNN([2] + [40] * 4 + [2], "tanh", "Glorot normal")
    model = dde.Model(data, net)

    model.compile("adam", lr=1e-3)
    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()

[lululxvi]

The grammar looks OK. Could you provide the PDE and BCs?

[tangqi]

Sorry, what do you mean? This was the entire code.

[lululxvi]

I mean the math equations.

[tangqi]

OK, this is some simplified version of some physics problem: z and v depends on theta and time dz_dt = v dv_dt = -dv_dtheta

To simplify the problem and identify the issue, there is no spatial dependence in the initial condition of (z0=0, v0=1). BC can be either dirichelt or periodic.

[tangqi]

FYI, the Dirichlet code runs fine, the periodic code has trouble around the boundary.

[lululxvi]

• First, there is a bug in the your periodic BC. It should be

  bc1 = dde.PeriodicBC(geomtime, 0, boundary_r, component=0)
  bc2 = dde.PeriodicBC(geomtime, 0, boundary_l, component=0)
  bc3 = dde.PeriodicBC(geomtime, 0, boundary_r, component=1)
  bc4 = dde.PeriodicBC(geomtime, 0, boundary_l, component=1)

  Note: For bc1, the first 0 means that this periodic BC is in the direction of the first coordinate, i.e., x, and component=0 means this BC is for y0.

• Second, it works if you train for a smaller time and use more points, e.g.,

  timedomain = dde.geometry.TimeDomain(0, 0.1)
  data = dde.data.TimePDE(..., num_initial =500)

  Hence, for solving longer time, it needs more training points, and perhaps some tuning to make the training loss small.

[tangqi]

Thanks a lot for the help! I will turn the parameters.

I want to make sure I understand you. Say I want to impose the third component of the solution vector along y direction in a 2D case, it would be dde.PeriodicBC(geomtime, 1, boundary_r, component=2)?

[lululxvi]

Yes, right.

[tangqi]

I did some tests. But it seems the training error decays very slowly in the periodic case while in the Dirichlet case decays very fast. So I am wondering if enforcing an extra constraint of u'(0)=u'(1) would be helpful. (I believe first order derivative should be a reasonable constraint to force for periodic, but probably not higher derivatives.)

You mentioned that it is tricky to implement it. But shouldn't it be just a combination of the current implementation of Neumann and Periodic boundary condition? If you can give some hint or point out some hidden caveat, we may be able to implement by ourself.

Thanks again, -qi

[lululxvi]

It makes sense that training with periodic BC is harder and slower than training with Dirichlet. Intuitively, when training with periodic BC, the information of the solution is only from IC; when training with periodic BC, the information of the solution is from both BC and IC. Even Neumann BC is harder to train than Dirichlet.

Yes, the implementation of first order periodic BC is a combination of Neumann and Periodic BC. It is not hard, but you need some understanding of the code.

[tangqi]

Thanks. After implemented the periodic bc to enforce u'(0)=u'(L), the above two cases work well. I was able to obtain the solutions almost as good as the exact solutions (no need to increase training points). Your package is quite clear, so it is not too hard to add one more boundary condition.

One more question: the training points consist of internal points, boundary points and collocation boundary points? I believe the collocation boundary points are used to enforce periodic boundary condition. What happens if we use other boundary conditions, Dirichlet for instance? Are they still being generated (I believe so) and what cost function are used on those points? Or will they become redundant in such a case? Thanks.

[lululxvi]

Great to know that it works. If you would like to share the code of the new BC, you are welcomed to open a pull request, or you can just paste the code here.

Here are the details in DeepXDE:

1. Generate the required number of training points: (1) internal points, (2) boundary points, (3) initial points;
2. All the points are used to force PDE losses;
3. DeepXDE loops over all BCs and ICs, and for each BC/IC DeepXDE checks which points generated in Step 1 should be used for this BC, and thus apply the BC for these points. Hence, a point could be used in more than one BC.

[tangqi]

continue in https://github.com/lululxvi/deepxde/pull/27

[athtareq]

Hello @tangqi, can you post the code and PDE problem (equations and boundary conditions) that worked for you here ? I'd like to understand more about PeriodicBC and it looks like your program is very extensive. Thanks !

[Fayaud123]

I think I am making an error on including the periodicity condition. This is the Burger equation in one dimension. I only want to include u(0, x)=u(1, x). Can someone see the mistake? thank you!

from future import absolute_import from future import division from future import print_function

import numpy as np

import deepxde as dde

import matplotlib.pyplot as plt

from deepxde.backend import tf nu = 0.01/np.pi a=1 b=5

'Define the exact solution' def func(x): return np.sin(anp.pix[:,1:2])np.cos(bx[:,0:1])

def pde(x, y): dy_x = dde.grad.jacobian(y, x, i=0, j=0) dy_t = dde.grad.jacobian(y, x, i=0, j=1) dy_xx = dde.grad.hessian(y, x, i=0, j=0) return( dy_t + y dy_x - nu dy_xx - ( a(np.pi)tf.cos(anp.pix[:,1:2])tf.cos(bx[:,0:1])

• btf.sin(anp.pix[:,1:2])tf.cos(bx[:,0:1])tf.sin(anp.pix[:,1:2])tf.sin(bx[:,0:1]) +b*2nutf.sin(anp.pix[:,1:2])tf.cos(b*x[:,0:1]) ) )

def boundary_t0(x, on_boundary): return on_boundary and np.isclose(x[1], 0)

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

bc1 = dde.DirichletBC(geomtime, lambda x: func(x), lambda _, on_boundary: on_boundary)

bc2= dde.PeriodicBC(geomtime, 1, boundary_t0, component=0) #here is where I am not sure bc3= dde.PeriodicBC(geomtime, 1, boundary_t1, component=0) # and here data = dde.data.TimePDE(geomtime, pde, [bc1, bc2, bc3], num_domain=2540, num_boundary=80, num_initial=160)

[lululxvi]

Is the BC should be u(0, t)=u(1, t), or u(0, x)=u(1, x)? t is always the last coordinate.

[Fayaud123]

Yes thank you for replying. I want to set u(×, 0)= u(x, 1) for all x in the spatial domain. I am getting an error. Thanks.

Best,

Fayaud Mezatio PhD Student Mathematics and Statistics Ottawa University STEM building, #501

On Fri., Jun. 18, 2021, 9:20 p.m. Lu Lu, @.***> wrote:

Is the BC should be u(0, t)=u(1, t), or u(0, x)=u(1, x)? t is always the last coordinate.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/lululxvi/deepxde/issues/26#issuecomment-864336157, or unsubscribe https://github.com/notifications/unsubscribe-auth/AUQXM5MMBHSFR6QZDVPZWTTTTPWERANCNFSM4LZTKUMQ .

-- DISCLAIMER: The contents of this email and any attachments are confidential. They are intended for the named recipient(s) only. If you have received this email by mistake, please notify the sender immediately and you are herewith notified that the contents are legally privileged and that you do not have permission to disclose the contents to anyone, make copies thereof, retain or distribute or act upon it by any means, electronically, digitally or in print. The views expressed in this communication may be of a personal nature and not be representative of  AIMS-Rwanda and/or any of its Initiatives.

[lululxvi]

def boundary(x, on_boundary):
    return on_boundary

PeriodicBC(geomtime, 0, boundary)

[Fayaud123]

Thank you very much. I will try it.

Best,

On Tue., Jun. 22, 2021, 7:01 p.m. Lu Lu, @.***> wrote:

def boundary(x, on_boundary): return on_boundary PeriodicBC(geomtime, 0, boundary)

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/lululxvi/deepxde/issues/26#issuecomment-866391620, or unsubscribe https://github.com/notifications/unsubscribe-auth/AUQXM5P73GQXNYTJF23PWQTTUEI2ZANCNFSM4LZTKUMQ .

-- DISCLAIMER: The contents of this email and any attachments are confidential. They are intended for the named recipient(s) only. If you have received this email by mistake, please notify the sender immediately and you are herewith notified that the contents are legally privileged and that you do not have permission to disclose the contents to anyone, make copies thereof, retain or distribute or act upon it by any means, electronically, digitally or in print. The views expressed in this communication may be of a personal nature and not be representative of  AIMS-Rwanda and/or any of its Initiatives.

[Fayaud]

Sorry Professor Lu lu!

Why not: PeriodicBC(geomtime,1 , boundary)? we are setting the time periodicity condition which is the second coordinate. Or am I missing something?

best,

[lululxvi]

@Fayaud you are right. It should be 1.

[Fayaud]

Thank you very much Professor @lululxvi

[AishwaryaKOLANCHIRAJAN]

hello@lululxvi if i want to solve 1d problem ( in terms of x and t) where i want to impose periodic boundary condition as"2sech(2)*t" then should i have my periodic bc as follows: bc1 = dde.PeriodicBC(geomtime, 1, boundary_r, component=0)

[lululxvi]

What do you mean by "periodic boundary condition as "2sech(2)*t""?

[AishwaryaKOLANCHIRAJAN]

Dear sir, The exact solution is given as follows: ψ(x, t) = A sech [A(x − x0 )] exp i(A*A /2)t , where A =1, x0 = 0 The initial condition is given as follows : ψ(x, t = 0) = N sech(x) where N=2

Now if i want to impose periodic boundary condition as ψ (−L, t) = ψ (L, t), t ∈ [t0 , T ] (periodic boundary condition)

then should i have

bc_u0 = dde.icbc.PeriodicBC(geomtime, 0, lambda , on_boundary: on_boundary, derivative_order=0, component=0) bc_v0 = dde.icbc.PeriodicBC(geomtime, 0, lambda , on_boundary: on_boundary, derivative_order=0, component=1)

On Thu, Dec 29, 2022 at 10:41 AM Lu Lu @.***> wrote:

What do you mean by "periodic boundary condition as "2sech(2)*t""?

— Reply to this email directly, view it on GitHub https://github.com/lululxvi/deepxde/issues/26#issuecomment-1367076685, or unsubscribe https://github.com/notifications/unsubscribe-auth/AK2JBK5GRVH6KRFTFJSMYZTWPUMOPANCNFSM4LZTKUMQ . You are receiving this because you commented.Message ID: @.***>

[lululxvi]

Yes.