1D and 2D wave Equation with a Source Term

[ShaikhaTheGreen]

Hello @lululxvi and @jeremyyu8 !!
Thanks for the great enhancement on PINN!

I am modifying burgers.py to run for the 1D wave equation.

du_tt - c**2 * du_xx - S

where,

S(0,t) = sin (pi *f *t) #Source function at point 0

I have a couple of questions to better understand the model.

So I changed the gPINN pde to be:

def pde(x, y):
    c = 1
    ###############################Source##################################
    Amp = 1
    frequency = 1
    sigma = 0.5  #width of source signal
    source_x_coord = 0
    Gaussian_impulse =  Amp * tf.exp(-(1/(sigma**2))*(x[:,0:1]-source_x_coord)**2)
    S = Gaussian_impulse *tf.sin( 1* np.pi  * frequency * x[:, 1:2] )
    ###############################PDE####################################
    #Using new DeepXDE Jacobians and Hessians
    du_xx = dde.grad.hessian(u, x, i=0, j=0)
    du_tt = dde.grad.hessian(u, x, i=1, j=1)

    du_ttx = dde.grad.jacobian(du_tt, x, j=0)
    du_xxx = dde.grad.jacobian(du_xx, x, j=0)
    dS_x = -4 * x[0:1] * tf.exp(-2*x[0:1]**2) * tf.sin(6.28318530717959*x[1:2])

    du_ttt = dde.grad.jacobian(du_tt, x, j=1)
    du_xxt = dde.grad.jacobian(du_xx, x, j=1)
    dS_t = 6.28318530717959 * tf.exp(-2*x[0:1]**2) * tf.cos(6.28318530717959*x[1:2])

    return [
        du_tt - (c**2) * du_xx - S,
        du_ttx - (c**2) * du_xxx - dS_x,
        du_ttt - (c**2) * du_xxt - dS_t,
    ]

1. Please let me know if this change seems correct to you?
2. Can you please elaborate more on how to use output_transform(x, y) ? And how can I change it to suit my 1D wave equation?
3. Also, to expand to a 2D wave equation with a time-dependent source term what are the changes that I need to do? This will be a problem with 2D (x,y) +1D (t).

Thank you very much!
I'm looking forward to your reply.

[jeremyyu8]

Hi @engsbk,

The new loss terms are correct but there seem to be errors in the expressions for dS_x and dS_t, as you use x[0:1] instead of x[:,0:1], etc. The output transform is used to enforce the boundary conditions using hard constraints, so in this case you can have an output transform that returns x * NN(x,t) + sin(pift), so that whenever x = 0 the expression just becomes sin(pift) as you'd like. To expand to 2D (x,y) + 1D (t), you just need one more additional loss term, as in you'll have df/dx, df/dy, df/dt as 3 additional loss terms to the standard PINN loss term.

Hope this helps!

[ShaikhaTheGreen]

Thank you so much @jeremyyu8 for the clarifying explanation and advice!

The results I got seem to be "amplified" of some sort...

I got
lr = 1e-6
epochs = 15000

Any suggestions on how to improve this?

:::Update:::
I tried adding observation points around x = 0, like this:

##########################Source observed############################
point_location = 0
x_source = np.zeros(500) + point_location
t_source = np.linspace(0, 5, num=500)
observe_x = np.vstack((x_source, t_source)).T

Then I included that in the model without using the output transform, like:

folds = 5
data = dde.data.TimePDE(
    geomtime, pde, 
    [bc, ic], 
    num_domain = int(folds*2500), 
    num_boundary = int(folds*600), 
    num_initial = int(folds*3200),
    anchors = observe_x
)
gPINNRARmodel = dde.Model(data, net)
gPINNRARmodel.compile("adam", lr=1.0e-6, loss_weights=[1, 0.0001, 0.0001, 1, 1])
losshistory, train_state = gPINNRARmodel.train(epochs=15000)
gPINNRARmodel.compile("L-BFGS-B", loss_weights=[1, 0.0001, 0.0001, 1, 1])
losshistory, train_state = gPINNRARmodel.train()

I had the same number of epochs and training points in both the previous run and this one. I still think the L2 error can be reduced somehow.
Please advise.

[jeremyyu8]

To further decrease the L2 error, perhaps you can try using the Residual-based Adaptive Refinement (RAR) method, which is shown in burgers.py, which adaptively adds training points during training to locations with high PDE residual.

[ShaikhaTheGreen]

To further decrease the L2 error, perhaps you can try using the Residual-based Adaptive Refinement (RAR) method, which is shown in burgers.py, which adaptively adds training points during training to locations with high PDE residual.

I was already using RAR when I produced the above errors. Do you have any suggestions to improve RAR parameters (patience, min_delta, or number of random points)?

[jeremyyu8]

I'm surprised the error didn't decrease more even though you used RAR. Would you mind sharing what the training loss looks like?

[ShaikhaTheGreen]

Sure.

For the 1D case, I did not use the output transform function because it disturbed the results somehow:

This is PINNs:

and this is gPINNS:

I also tried it for 2D wave function:
This is the PINN training:

and this is the gPINN's:

I can also send you the python notebooks in an email for you to run.

Thank you for your reply!
Awaiting your feedback.

[ShaikhaTheGreen]

For the 1D case, I did not use the output transform function because it disturbed the results somehow:

I mean disturbed like this:

[ShaikhaTheGreen]

Hi @engsbk,

The new loss terms are correct but there seem to be errors in the expressions for dS_x and dS_t, as you use x[0:1] instead of x[:,0:1], etc. The output transform is used to enforce the boundary conditions using hard constraints, so in this case you can have an output transform that returns x * NN(x,t) + sin(pi_f_t), so that whenever x = 0 the expression just becomes sin(pi_f_t) as you'd like. To expand to 2D (x,y) + 1D (t), you just need one more additional loss term, as in you'll have df/dx, df/dy, df/dt as 3 additional loss terms to the standard PINN loss term.

Hope this helps!

I went back to the Burger's example, to try improve my case. I have a question about the tranform function in Burger's case:

def output_transform(x, y):
    x_in = x[:, 0:1]
    t_in = x[:, 1:2]

    return (1 - x_in) * (1 + x_in) * (1 - tf.exp(-t_in)) * y - tf.sin(np.pi * x_in)

I figured that (1 - x_in) * (1 + x_in) controls the BC locations,
but what does the term (1 - tf.exp(-t_in)) control?

[lululxvi]

It is IC.

[ShaikhaTheGreen]

It is IC.

Thanks for the reply @lululxvi .
I have another question regarding the loss weights assigned when training gPINN. How are the loss weights being set?
In the Burgers example the gradients weights are set to 0.0001 but in other examples it's sometimes 0.01 or 0.1.
What could be the best value for the 1D or 2D wave equation?

My 1D wave gPINNpde looks like this:

def gPINNpde(x, u):
    c = 1
    
    #Using new DeepXDE Jacobians and Hessians
    du_xx = dde.grad.hessian(u, x, i=0, j=0)
    du_tt = dde.grad.hessian(u, x, i=1, j=1)

    du_ttx = dde.grad.jacobian(du_tt, x, j=0)
    du_xxx = dde.grad.jacobian(du_xx, x, j=0)
    
    du_ttt = dde.grad.jacobian(du_tt, x, j=1)
    du_xxt = dde.grad.jacobian(du_xx, x, j=1)
    
    return [
        du_tt - (c**2) * du_xx,
        du_ttx - (c**2) * du_xxx,
        du_ttt - (c**2) * du_xxt,
    ]

[lululxvi]

The best weight is problem dependent. You have to tune the weight.

[ShaikhaTheGreen]

I have another question about the transform function and hard constraints. If I have my transform function as follows:

def output_transform(x, y):
    x_in = x[:, 0:1]
    t_in = x[:, 1:2]

    return (x_in) * (5 - x_in) * (1 - tf.exp(-t_in)) * y - tf.sin(np.pi * t_in) + (5 - x_in) * tf.sin(2*np.pi*f*t_in)

To satisfy:

u(0,t) = sin(2 pi f t), and u(5,t) = 0,
u(x,0) = 0

If the values of x = 0 or 5, then the equation satisfies my BC, but if I have middle values such as 2,3,...etc., this will cause contamination to the prediction results of the NN( multiplies the prediction of the NN with a number 5-(2 or 3) = 3 or 2). I also face the same problem for middle values of t.

How can I avoid that?

[lululxvi]

I cannot understand your question.

[ShaikhaTheGreen]

I'll try to explain with this example that I'm trying to run. In the standard PINN configuration, I tried to run this example:

################################################################
########################HARD IC & BC############################
################################################################
def output_transform(x, u):
    x_in = x[:, 0:1]
    t_in = x[:, 1:2]
    freq = 1

    return  (1-x_in) * (x_in) * (t_in) * u + (1-x_in ) * tf.sin(4*np.pi * freq * t_in )

to apply these hard constraints:

u(0,t) = sin(4 pi f t), and u(1,t) = 0,
u(x,0) = 0

In this domain:

geom = dde.geometry.Interval(0,1)
timedomain = dde.geometry.TimeDomain(0,1)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

For this wave PDE:

def pde(x, u):
    c = 1
    #Using new DeepXDE Jacobians and Hessians
    du_xx = dde.grad.hessian(u, x, i=0, j=0)
    du_tt = dde.grad.hessian(u, x, i=1, j=1)
    
    return   du_tt - (c ** 2)*(du_xx)

My PINN model:

d = 2 #dimensions
level = 3  #dataset density level
nx_train = 10 * 2 ** (level - 1)
data = dde.data.TimePDE(
        geomtime, pde, 
        [], 
        num_domain = nx_train**d, 
        num_boundary = 2 * d * nx_train, 
        num_initial = nx_train**d,
        train_distribution="Sobol"
    ) 
n = 5
activation = f"LAAF-{n} sin"
net = dde.maps.FNN([2] + [75] * 7 + [1], activation, "Glorot normal")
net.apply_output_transform(output_transform)
the_model = dde.Model(data, net)
model.append(the_model)
model[model_i].compile("adam", lr=4.9e-5)
model[model_i].train(epochs=10000)
dde.optimizers.set_LBFGS_options(ftol=1e-18, gtol=1e-18, maxiter=100000)
model[model_i].compile("L-BFGS-B")
losshistory, train_state = model[model_i].train()
dde.saveplot(losshistory, train_state, issave=True, isplot=True)

After running this example 10 times I got these results ( I used 10 trials to make sure it did not have to do with randomization):

I'm not sure why the prediction looks a bit out of phase in the beginning. What would you recommend to fix this?