Helmholtz equation with Neumann boundary conditions over a 2D square domain
Problem setup
For a wavenumber \(k_0 = 2\pi n\) with \(n = 2\), we will solve a Helmholtz equation:
with the Neumann boundary conditions
with \(n\) the normal exterior vector and a source term \(f(x,y) = k_0^2 \cos(k_0 x)\cos(k_0 y)\).
Remark that the exact solution reads:
This example is the Neumann boundary condition conterpart to this Dolfinx tutorial. One can refer to Ihlenburg's book "Finite Element Analysis of Acoustic Scattering" p138-139 for more details.
Implementation
This description goes through the implementation of a solver for the above described Helmholtz equation step-by-step.
First, the DeepXDE and Numpy modules are imported:
import deepxde as dde
import numpy as np
We begin by definying the general parameters for the problem. We use a collocation points density of 10 (resp. 30) points per wavelength for the training (resp. testing) data along each direction. This code allows to use both soft and hard boundary conditions.
n = 2
precision_train = 10
precision_test = 30
weights = 100
The PINN will be trained over 5000 epochs. We define the learning rate, the number of dense layers and nodes, and the activation function. Moreover, we import the cosine function.
epochs = 5000
parameters = [1e-3, 3, 150, "sin"]
# Define sine function
if dde.backend.backend_name == "pytorch":
cos = dde.backend.pytorch.cos
else:
from deepxde.backend import tf
cos = tf.cos
learning_rate, num_dense_layers, num_dense_nodes, activation = parameters
Next, we express the PDE residual of the Helmholtz equation:
def pde(x, y):
dy_xx = dde.grad.hessian(y, x, i=0, j=0)
dy_yy = dde.grad.hessian(y, x, i=1, j=1)
f = k0 ** 2 * cos(k0 * x[:, 0:1]) * cos(k0 * x[:, 1:2])
return -dy_xx - dy_yy - k0 ** 2 * y - f
The first argument to pde is the network input, i.e., the \(x\)-coordinate and \(y\)-coordinate. The second argument is the network output, i.e., the solution \(u(x)\), but here we use y as the name of the variable.
Next, we introduce the exact solution and the Neumann boundary condition.
def func(x):
return np.cos(k0 * x[:, 0:1]) * np.cos(k0 * x[:, 1:2])
def boundary(_, on_boundary):
return on_boundary
Now, we define the geometry and evaluate the number of training and test random collocation points. The values allow to obtain collocation points density of 10 (resp. 30) points per wavelength along each direction. We define the boundary and the Neumann boundary conditions.
geom = dde.geometry.Rectangle([0, 0], [1, 1])
k0 = 2 * np.pi * n
wave_len = 1 / n
hx_train = wave_len / precision_train
nx_train = int(1 / hx_train)
hx_test = wave_len / precision_test
nx_test = int(1 / hx_test)
bc = dde.icbc.NeumannBC(geom, lambda x: 0, boundary)
Next, we generate the training and testing points.
data = dde.data.PDE(
geom,
pde,
bc,
num_domain=nx_train ** 2,
num_boundary=4 * nx_train,
solution=func,
num_test=nx_test ** 2,
)
Next, we choose the network. Here, we use a fully connected neural network of depth 4 (i.e., 3 hidden layers) and width 150. Besides, we choose sin as activation function and Glorot uniform as initializer :
net = dde.nn.FNN(
[2] + [num_dense_nodes] * num_dense_layers + [1], activation, "Glorot uniform"
)
Now, we have the PDE problem and the network. We build a Model and define the optimizer and learning rate.
model = dde.Model(data, net)
if hard_constraint == True:
model.compile("adam", lr=learning_rate, metrics=["l2 relative error"])
else:
loss_weights = [1, weights]
model.compile(
"adam",
lr=learning_rate,
metrics=["l2 relative error"],
loss_weights=loss_weights,
)
We first train the model for 5000 iterations with Adam optimizer:
losshistory, train_state = model.train(epochs=epochs)
Complete code
"""Backend supported: tensorflow.compat.v1, tensorflow, pytorch"""
import deepxde as dde
import numpy as np
# General parameters
n = 2
precision_train = 10
precision_test = 10
hard_constraint = True
weights = 100 # if hard_constraint == False
epochs = 5000
parameters = [1e-3, 3, 150, "sin"]
# Define sine function
if dde.backend.backend_name == "pytorch":
cos = dde.backend.pytorch.cos
else:
from deepxde.backend import tf
cos = tf.cos
learning_rate, num_dense_layers, num_dense_nodes, activation = parameters
def pde(x, y):
dy_xx = dde.grad.hessian(y, x, i=0, j=0)
dy_yy = dde.grad.hessian(y, x, i=1, j=1)
f = k0**2 * cos(k0 * x[:, 0:1]) * cos(k0 * x[:, 1:2])
return -dy_xx - dy_yy - k0**2 * y - f
# WORK / TODO /WORKING ON TRAVIS WHY
def func(x):
return np.cos(k0 * x[:, 0:1]) * np.cos(k0 * x[:, 1:2])
def boundary(_, on_boundary):
return on_boundary
geom = dde.geometry.Rectangle([0, 0], [1, 1])
k0 = 2 * np.pi * n
wave_len = 1 / n
hx_train = wave_len / precision_train
nx_train = int(1 / hx_train)
hx_test = wave_len / precision_test
nx_test = int(1 / hx_test)
bc = dde.icbc.NeumannBC(geom, lambda x: 0, boundary)
data = dde.data.PDE(
geom,
pde,
bc,
num_domain=nx_train**2,
num_boundary=4 * nx_train,
solution=func,
num_test=nx_test**2,
)
net = dde.nn.FNN(
[2] + [num_dense_nodes] * num_dense_layers + [1], activation, "Glorot uniform"
)
model = dde.Model(data, net)
loss_weights = [1, weights]
model.compile(
"adam",
lr=learning_rate,
metrics=["l2 relative error"],
loss_weights=loss_weights,
)
losshistory, train_state = model.train(epochs=epochs)
dde.saveplot(losshistory, train_state, issave=True, isplot=True)