2D acoustic wave equation with FDTD data for training?

[itskobold]

Hello,

I'm looking at implementing a PINN to solve the PDE dy_tt - C ** 2 * (dy_xx + dy_yy) where C is wave propagation speed. The initial condition is a simple Dirac pulse at the center of the rectangular domain and boundary conditions are dirichlet for now. These are both "hardly" enforced using net.apply_output_transform().

As far as I'm aware, I can't provide a solution function to the network as the domain must be discretised for solution. Instead, I have arrays of FDTD simulation data in the form (x, y, t). Is there a way to reference this data as the "solution" for the network to train from? I've looked at some examples for other PDE implementations but am still a little lost.

Thanks in advance!

[lululxvi]

You don't have to use the solution function.

[lululxvi]

Yes

[itskobold]

Great - thanks for your help. I'll mark this as answered

[itskobold]

Hi again Lu,

I've got the PINN solution to the wave equation coded up in DeepXDE and it seems to be converging.
However, I'm struggling to obtain the result I'm expecting (a 'blob' at the center of the square 2D domain which ripples out and reflects from the boundaries). The boundary conditions are dirichlet and enforced hardly, the ICs are softly enforced via loss functions.
I'm starting to think that the problem requires far too many sampling points to be tractable. There are a lot of discontinuities in the solution caused by waves colliding with each other after reflecting from the boundaries and I've read elsewhere that convergence of PINNs often suffer in situations like this.
I've attached the code I'm working with. Are there any problems that you can see or advice you could offer? After 45,000 training epochs I can obtain a very rough (but encouraging) model of waves rippling out from the center. The combined loss products after training were around 1e-07. I've uploaded a video of this attempt here: https://www.youtube.com/shorts/5ELwG6d4fsA

Thanks again for your assistance.

dde.config.set_default_float("float64")

x_len = 1
y_len = 1
lr = 0.001
C = 1
epochs_adam = 25000
epochs_l_bfgs_b = 20000

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

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

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

def initial(x, y):
a = 1
b_x = 0
b_y = 0
c = 0.02

top_half = (x - b_x) ** 2 + (y - b_y) ** 2
bottom_half = 2 * c ** 2
return a * np.e ** -(top_half / bottom_half)

def func(x):
x, y, _ = np.split(x, 3, axis=1)
return initial(x, y)

ic_1 = dde.icbc.IC(geomtime, func, lambda _, on_initial: on_initial)
ic_2 = dde.icbc.OperatorBC(
geomtime,
lambda x, y, _: dde.grad.jacobian(y, x, i=0, j=1),
lambda x, _: np.isclose(x[1], 0),
)

data = dde.data.TimePDE(
geomtime,
pde,
[ic_1, ic_2],
num_domain=10000,
num_boundary=1000,
num_initial=1000,
num_test=1000,
)

layer_size = [3] + [200] + [100] * 5 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.nn.STMsFFN(
layer_size, activation, initializer, sigmas_x=[1], sigmas_t=[1, 10], dropout_rate=0
)

max_len = max(x_len, y_len)
net.apply_feature_transform(lambda x: (x - max_len) * 2 * np.sqrt(3))

# Dirichlet boundary hard enforcement
net.apply_output_transform(
lambda x, y: (1 - x[:, 0:1] ** 2) *
(1 - x[:, 1:2] ** 2) *
y
)

# Create and train model
model = dde.Model(data, net)
initial_losses = get_initial_loss(model)
loss_weights = 5 / initial_losses

# First train with adam optimiser
model.compile(
"adam",
lr=lr,
loss_weights=loss_weights,
decay=("inverse time", 1000, 0.7),
)
pde_residual_resampler = dde.callbacks.PDEResidualResampler(period=1)
losshistory, train_state = model.train(
epochs=epochs_adam, callbacks=[pde_residual_resampler], display_every=100
)

# Then train with L-BFGS-B
dde.optimizers.set_LBFGS_options(gtol=1e-12, maxiter=epochs_l_bfgs_b)
model.compile("L-BFGS-B", loss_weights=losshistory.loss_weights)
losshistory, train_state = model.train(
callbacks=[pde_residual_resampler], display_every=100
)

[lululxvi]

Discontinuity is always difficult to handle, and I don't have a silver bullet for this.

[JinglaiZheng]

Have you solved this problem? I'm struggling to solve hyperbolic equations, but the model loss converges to a very high value,thanks. @itskobold

[whtwu]

Hello,

I'm looking at implementing a PINN to solve the PDE dy_tt - C ** 2 * (dy_xx + dy_yy) where C is wave propagation speed. The initial condition is a simple Dirac pulse at the center of the rectangular domain and boundary conditions are dirichlet for now. These are both "hardly" enforced using net.apply_output_transform().

As far as I'm aware, I can't provide a solution function to the network as the domain must be discretised for solution. Instead, I have arrays of FDTD simulation data in the form (x, y, t). Is there a way to reference this data as the "solution" for the network to train from? I've looked at some examples for other PDE implementations but am still a little lost.

Thanks in advance!

Hello, I notice that you use Dirac pulse, and i want to use the impulse excitation as the initial condition on my beam problem, can you show me your IC code? thanks a lot.

[TethysSun]

@lululxvi Hi Lu, just a quick question. How to include an exogenous factor for pde? for example, the equation looks like: A + B * dy_tt - C ** 2 * (dy_xx + dy_yy), where A, B and C are tensors with size nx * ny?

A,B,C cannot be computed as functions of x, y and t. If this is the case, how to introduce ABC in the training of PINN using the deepxde?

Thank you!

[lululxvi]

https://deepxde.readthedocs.io/en/latest/demos/pinn_inverse/lorenz.inverse.forced.html