`0` and `Nan` values for `IC` defined as a `NeumannBC`

[Saransh-cpp]

I have been trying to tweak the wave_1D example with different ICs and BCs. I went through all the FAQs and figured out that ICs with derivatives can be defined as NeumannBCs. I tried doing this in my code but for some reason, the loss value associated with that IC is always 0, and loss_weights makes the value nan (because of division by 0). Could someone please help me with this? Is there something that I am doing wrong?

The BCs and ICs

def boundary_initial(x, _):
    return np.isclose(x[1], 0)

bc = dde.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary)
ic_1 = dde.IC(
    geomtime,
    lambda x: np.sin(n * np.pi * x[:, 0:1] / L),
    lambda _, on_initial: on_initial,
)
ic_2 = dde.NeumannBC(geomtime, lambda x: 0, boundary_initial)

The complete code

import numpy as np
import pandas as pd
import deepxde as dde
from deepxde.backend import tf
import matplotlib.pyplot as plt

A = 2
c = 10
L = 1
C = 2
n = 1
m = 1

def pde(x, y):
    dy_tt = dde.grad.hessian(y, x, i=1, j=1)
    dy_xx = dde.grad.hessian(y, x, i=0, j=0)
    return dy_tt - c ** 2 * dy_xx

def sol(x):
    x, t = np.split(x, 2, axis=1)
    return C * np.sin(n * np.pi * x / L) * np.cos(m * np.pi * t * c / L)

def boundary_initial(x, _):
    return np.isclose(x[1], 0)

def get_initial_loss(model):
    model.compile("adam", lr=0.001, metrics=["l2 relative error"])
    losshistory, train_state = model.train(0)
    return losshistory.loss_train[0]

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

bc = dde.DirichletBC(geomtime, lambda x: 0, lambda _, on_boundary: on_boundary)
ic_1 = dde.IC(
    geomtime,
    lambda x: np.sin(n * np.pi * x[:, 0:1] / L),
    lambda _, on_initial: on_initial,
)
ic_2 = dde.NeumannBC(geomtime, lambda x: 0, boundary_initial)

data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc, ic_1, ic_2],
    num_domain=360,
    num_boundary=360,
    num_initial=360,
    num_test=10000,
    solution=sol,
)

layer_size = [2] + [100] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.STMsFFN(
    layer_size, activation, initializer, sigmas_x=[1], sigmas_t=[1, 10]
)
net.apply_feature_transform(lambda x: (x - 0.5) * 2 * np.sqrt(3))

model = dde.Model(data, net)
initial_losses = get_initial_loss(model)
loss_weights = 5 / initial_losses
losshistory, train_state = model.train(0)
model.compile(
    "adam",
    lr=0.001,
    metrics=["l2 relative error"],
    decay=("inverse time", 2000, 0.9),
    loss_weights=loss_weights,
)
pde_residual_resampler = dde.callbacks.PDEResidualResampler(period=1)
losshistory, train_state = model.train(epochs=10000, callbacks=[pde_residual_resampler])

dde.saveplot(
    losshistory,
    train_state,
    issave=True,
    isplot=True,
    test_fname="../dat_data/wave_1D.dat",
)

Training output

Initializing variables...
Training model...

Step      Train loss                                  Test loss                                   Test metric   
0         [7.26e+05, 4.77e-02, 2.58e-01, 0.00e+00]    [8.66e+05, 4.77e-02, 2.58e-01, 0.00e+00]    [9.77e-01]    

Best model at step 0:
  train loss: 7.26e+05
  test loss: 8.66e+05
  test metric: [9.77e-01]

'train' took 1.184103 s

D:\OpenSource\python\PDEsWithPINNs\physics_informed_spatio_temporal_multi_scale_fourier_feature_networks\wave_1D.py:80: RuntimeWarning: divide by zero encountered in true_divide
  loss_weights = 5 / initial_losses
Initializing variables...
Training model...

Step      Train loss                                  Test loss                                   Test metric   
0         [9.36e+05, 1.32e-01, 2.16e-02, 0.00e+00]    [6.49e+05, 1.32e-01, 2.16e-02, 0.00e+00]    [1.05e+00]    

Best model at step 0:
  train loss: 7.26e+05
  test loss: 8.66e+05
  test metric: [9.77e-01]

'train' took 0.156740 s

Compiling model...
'compile' took 1.140433 s

Initializing variables...
Training model...

Step      Train loss                                  Test loss                                   Test metric
0         [6.36e+00, 1.01e+01, 7.51e+00, nan]         [8.70e+00, 1.01e+01, 7.51e+00, nan]         [1.06e+00]
1000      [6.64e-03, 2.05e-04, 2.01e-05, nan]         [1.80e-02, 2.05e-04, 2.01e-05, nan]         [8.39e-01]
2000      [2.77e-03, 8.66e-05, 1.32e-05, nan]         [7.61e-03, 8.66e-05, 1.32e-05, nan]         [7.39e-01]
3000      [1.33e-03, 1.88e-04, 8.96e-05, nan]         [3.54e-03, 1.88e-04, 8.96e-05, nan]         [6.52e-01]
4000      [5.32e-04, 5.26e-05, 4.56e-06, nan]         [1.55e-03, 5.26e-05, 4.56e-06, nan]         [5.83e-01]
5000      [3.36e-04, 1.34e-03, 1.55e-03, nan]         [7.93e-04, 1.34e-03, 1.55e-03, nan]         [5.43e-01]
6000      [1.92e-04, 2.09e-04, 7.03e-05, nan]         [5.18e-04, 2.09e-04, 7.03e-05, nan]         [5.22e-01]
7000      [1.09e-04, 5.56e-05, 2.22e-05, nan]         [3.45e-04, 5.56e-05, 2.22e-05, nan]         [5.14e-01]
8000      [1.03e-04, 7.60e-05, 1.66e-04, nan]         [2.82e-04, 7.60e-05, 1.66e-04, nan]         [5.10e-01]
9000      [9.47e-05, 2.85e-04, 6.99e-05, nan]         [2.46e-04, 2.85e-04, 6.99e-05, nan]         [5.08e-01]
10000     [5.59e-05, 1.19e-05, 8.72e-07, nan]         [1.92e-04, 1.19e-05, 8.72e-07, nan]         [5.07e-01]

Best model at step 0:
  train loss: 7.26e+05
  test loss: 8.66e+05
  test metric: [9.77e-01]

'train' took 584.973021 s

[Saransh-cpp]

I'll try using OperatorBC here. Thank you, @lululxvi! Though I was able to get rid of nan values by adding a small factor in the denominator in loss_weights = 5 / initial_losses.

[Saransh-cpp]

Thank you, @lululxvi! This was solved by using OperatorBC.

[tjboise]

I'll try using OperatorBC here. Thank you, @lululxvi! Though I was able to get rid of nan values by adding a small factor in the denominator in loss_weights = 5 / initial_losses.

Hi Saransh, I think I meet the same problem with you. Would you like show me how to use the OperatorBC?

[Saransh-cpp]

@TianjieZhang1993 you can do something like this -

def boundary_initial(x, _):
    return np.isclose(x[1], 0)  # when t = 0

ic_2 = dde.OperatorBC(
    geomtime, lambda x, u, _: dde.grad.jacobian(u, x, i=0, j=1), boundary_initial,
)

OperatorBc defines the BC in the way -

func(inputs, outputs, X) = 0.

[tjboise]

@TianjieZhang1993 you can do something like this -

def boundary_initial(x, _):
    return np.isclose(x[1], 0)  # when t = 0

ic_2 = dde.OperatorBC(
    geomtime, lambda x, u, _: dde.grad.jacobian(u, x, i=0, j=1), boundary_initial,
)

OperatorBc defines the BC in the way -

func(inputs, outputs, X) = 0.

Thanks so much! Though I still figure out the problem in my code~

[engsbk]

I'm trying to tweak that example for a 2D space wave in 1D time. I used:

S_BC = dde.icbc.OperatorBC(
    geomtime,
    lambda x, u, _: Amp*tf.sin(2*np.pi*freq*x[:,2:3]),
    lambda x, _: np.isclose(x[0],0.5) and np.isclose(x[1], 2.5),
)

def get_initial_loss(model):
    model.compile("adam", lr=0.001)
    losshistory, train_state = model.train(0)
    return losshistory.loss_train[0]

def pde(x, u):
    c = 1 
    du_tt = dde.grad.hessian(u, x, i=2, j=2)
    du_xx = dde.grad.hessian(u, x, i=0, j=0)
    du_yy = dde.grad.hessian(u, x, i=1, j=1)

    return du_tt - (c ** 2) * (du_xx + du_yy) 

def func(x):
    return 0

t_end = 5
d_start = 1
d_end = 5
geom = dde.geometry.Rectangle([d_start,d_start], [d_end,d_end])
timedomain = dde.geometry.TimeDomain(0, t_end)
geomtime = dde.geometry.GeometryXTime(geom, timedomain)

bc = dde.icbc.DirichletBC(geomtime, func, lambda _, on_boundary: on_boundary)
ic = dde.icbc.IC(geomtime, func, lambda _, on_initial: on_initial)
# do not use dde.NeumannBC here, since `normal_derivative` does not work with temporal coordinate.
Amp = 2
freq = 1
S_BC = dde.icbc.OperatorBC(
    geomtime,
    lambda x, u, _: Amp*tf.sin(2*np.pi*freq*x[:,2:3]),
    lambda x, _: np.isclose(x[0],0.5) and np.isclose(x[1], 2.5),
)
data = dde.data.TimePDE(
    geomtime,
    pde,
    [bc, ic, S_BC],
    num_domain=500,
    num_boundary=500,
    num_initial=500,
)

layer_size = [3] + [100] * 3 + [1]
n = 10
activation = "tanh" #f"LAAF-{n} tanh"  # "LAAF-10 tanh"
initializer = "Glorot uniform"
net = dde.nn.STMsFFN(
    layer_size, activation, initializer, sigmas_x=[1, 1], sigmas_t=[1]
)

model = dde.Model(data, net)
initial_losses = get_initial_loss(model)
loss_weights = 5 / initial_losses
model.compile(
    "adam",
    lr=0.001,
    loss_weights=loss_weights,
    decay=("inverse time", 2000, 0.9),
)
pde_residual_resampler = dde.callbacks.PDEResidualResampler(period=1)
losshistory, train_state = model.train(
    epochs=10000, callbacks=[pde_residual_resampler], display_every=500
)

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

But I received this error:

Initializing variables..
.
Training model...

Step      Train loss                                  Test loss                                   Test metric
0         [8.35e-01, 6.64e-02, 1.05e-01, nan]         [8.35e-01, 6.64e-02, 1.05e-01, nan]         []  

Best model at step 0:
  train loss: inf
  test loss: inf
  test metric: 

'train' took 0.782480 s

Compiling model...
'compile' took 2.669732 s

Initializing variables...
Training model...

Step      Train loss                                  Test loss                                   Test metric
0         [4.02e+00, 4.19e+00, 3.88e+00, nan]         [4.02e+00, 4.19e+00, 3.88e+00, nan]         []  

Best model at step 0:
  train loss: inf
  test loss: inf
  test metric: 

'train' took 2.134405 s

Saving loss history to /content/loss.dat ...
Saving training data to /content/train.dat ...
Saving test data to /content/test.dat ...

---------------------------------------------------------------------------

ValueError                                Traceback (most recent call last)

[<ipython-input-4-f4f9f2108362>](https://localhost:8080/#) in <module>()
     47 )
     48 
---> 49 dde.saveplot(losshistory, train_state, issave=True, isplot=True)

2 frames

[/usr/local/lib/python3.7/dist-packages/deepxde/utils/external.py](https://localhost:8080/#) in saveplot(loss_history, train_state, issave, isplot, loss_fname, train_fname, test_fname, output_dir)
    173         test_fname = os.path.join(output_dir, test_fname)
    174         save_loss_history(loss_history, loss_fname)
--> 175         save_best_state(train_state, train_fname, test_fname)
    176 
    177     if isplot:

[/usr/local/lib/python3.7/dist-packages/deepxde/utils/external.py](https://localhost:8080/#) in save_best_state(train_state, fname_train, fname_test)
    342     print("Saving test data to {} ...".format(fname_test))
    343     if y_test is None:
--> 344         test = np.hstack((train_state.X_test, best_y))
    345         if best_ystd is None:
    346             np.savetxt(fname_test, test, header="x, y_pred")

<__array_function__ internals> in hstack(*args, **kwargs)

[/usr/local/lib/python3.7/dist-packages/numpy/core/shape_base.py](https://localhost:8080/#) in hstack(tup)
    343         return _nx.concatenate(arrs, 0)
    344     else:
--> 345         return _nx.concatenate(arrs, 1)
    346 
    347 

<__array_function__ internals> in concatenate(*args, **kwargs)

ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 1 has 1 dimension(s)

Any idea how to solve this?

[lululxvi]

Do you have residual points for S_BC?