1D Helmholtz equation with NBC and RBC

[juliandwain]

Hello,

I am currently trying to implement the Helmholtz equation in 1D (evaluating an acoustical problem) given as:

with a NBC at the left end and a RBC at the right end of the interval.
The value of the NBC equals and the value of the RBC equals
.
In general, the solution given the mentioned BCs is stated as .
Since only the real part is considered, the solution is given as .
The code I used is shown below, but it is not working correctly and I don't know what I did wrong.

# -*- coding: utf-8 -*-
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import deepxde as dde
import numpy as np
import tensorflow as tf


def main():
    length = 3.4

    v_s = [-1, 0]  # structural velocity at x=0 and x=length in [m/s]
    rho = 1.3  # fluid density in [kg/m^3]
    c = 340  # speed of sound in [m/s]
    admittance = [0, 1 / (rho * c)]  # boundary admittance at x=0 and x=length in [m^3/(N s)]
    frequency = 200  # frequency in [Hz]
    wavelength = c / frequency  # wavelength in [m]
    k = (2 * np.pi * frequency) / c  # wavenumber in [m^-1]

    def pde(x, y):
        """Helmholtz PDE 1D.

        Parameters
        ----------
        x :
            x-coordinates:
        y :
            Sound pressure solution of Helmholtz PDE.

        Returns
        -------

        """
        dy_x = tf.gradients(y, x)[0]
        dy_xx = tf.gradients(dy_x, x)[0]
        return dy_xx + (k ** 2) * y

    def boundary_left(x, on_boundary):
        """Neumann BC on left boundary.

        Parameters
        ----------
        x
        on_boundary

        Returns
        -------

        """
        return on_boundary and np.isclose(x[0], 0)

    def boundary_right(x, on_boundary):
        """Robin BC on right boundary.

        Parameters
        ----------
        x
        on_boundary

        Returns
        -------

        """
        return on_boundary and np.isclose(x[0], length)

    def func(x):
        """Reference solution real part.

        Parameters
        ----------
        x

        Returns
        -------

        """
        return -v_s[0] * rho * c * np.cos(k * x)

    geom = dde.geometry.Interval(0, length)
    bc_left = dde.NeumannBC(geom, lambda X: v_s[0], boundary_left)
    bc_right = dde.RobinBC(geom, lambda X, y: admittance[1] * y + v_s[1], boundary_right)
    data = dde.data.PDE(geom, pde, [bc_left, bc_right], num_domain=16, num_boundary=2, solution=func, num_test=100)

    layer_size = [1] + [50] * 3 + [1]
    activation = "tanh"
    initializer = "Glorot uniform"
    net = dde.maps.FNN(layer_size, activation, initializer)

    model = dde.Model(data, net)
    model.compile("adam", lr=0.001, metrics=["l2 relative error"])  # metrics=["l2 relative error"]
    losshistory, train_state = model.train(epochs=10000)

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


if __name__ == "__main__":
    main()

Thank you very much for your help!

[lululxvi]

Could you try:

bc_left = dde.NeumannBC(geom, lambda X: np.full((len(X), 1), v_s[0]), boundary_left)

[juliandwain]

Hi Lu,

thank you very much for your quick answer, I edited the file but it does not seem to work correctly.
The loss for the training data is at O(1e-5), the loss for the test data is at O(1e0), and the test metric stays constant at 1.
It is always predicting only 0 at each position.
Acutally, I don't really understand why this is happening like this, the configuration, the analytical solution as well as the BCs are (in my best knowledge) correctly.
This is the output I get when training:

Step      Train loss                        Test loss                         Test metric   
0         [1.10e+01, 1.42e+00, 4.58e-04]    [1.10e+01, 0.00e+00, 0.00e+00]    [1.00e+00]    
1000      [8.77e-02, 7.06e-01, 8.00e-05]    [8.66e-02, 0.00e+00, 0.00e+00]    [1.00e+00]    
2000      [2.34e-01, 1.66e-01, 4.52e-03]    [2.34e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
3000      [2.08e-01, 1.23e-01, 6.97e-04]    [2.03e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
4000      [1.89e-01, 9.48e-02, 4.93e-04]    [1.86e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
5000      [1.75e-01, 6.23e-02, 5.99e-03]    [1.74e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
6000      [1.69e-01, 5.28e-02, 6.46e-03]    [1.68e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
7000      [1.64e-01, 4.90e-02, 4.97e-03]    [1.62e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
8000      [1.54e-01, 4.41e-02, 3.20e-03]    [1.51e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
9000      [1.41e-01, 3.64e-02, 5.05e-04]    [1.41e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
10000     [1.29e-01, 2.94e-02, 4.79e-04]    [1.29e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
Best model at step 10000:
  train loss: 1.59e-01
  test loss: 1.29e-01
  test metric: [1.00e+00]

Thank you very much for your help, if you can provide me any tips I would be very thankful.

[lululxvi]

Actually your train loss is too large. There are 3 columns in train loss: the first column is PDE loss, the second column is left BC loss, and the third column is right BC loss. So the main issue is that the PDE loss doesn't go down.

The possible reason is that your PDE solution is not of scale O(1), which makes the training hard. There are two solutions:

1. re-scale your PDE, such that the solution is O(1);
2. If you assume the PDE solution is of scale O(10), then multiply the network output by 10:

net.outputs_modify(lambda x, y: 10 * y))

[juliandwain]

Thank you for your help!
Unfortunately, as I do not know on how to scale the PDE properly for now (I will come back to this), I tried to modify my code according to your 2nd suggestions, but this did not work either.
I put your line in right after creating the FNN

layer_size = [1] + [50] * 3 + [1]
activation = "tanh"
initializer = "Glorot uniform"
net = dde.maps.FNN(layer_size, activation, initializer)
net.outputs_modify(lambda x, y: 10 * y)

and I did not change any other parameter, but the loss for the PDE is still too high:

Step      Train loss                        Test loss                         Test metric   
0         [2.62e+02, 4.45e+00, 2.97e-02]    [2.77e+02, 0.00e+00, 0.00e+00]    [1.00e+00]    
1000      [2.08e-01, 2.38e-01, 8.59e-04]    [1.86e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
2000      [1.85e-01, 1.04e-01, 5.35e-03]    [1.86e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
3000      [1.82e-01, 7.54e-02, 8.64e-03]    [1.92e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
4000      [1.82e-01, 7.53e-02, 7.41e-03]    [1.91e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
5000      [8.35e-01, 4.25e-02, 4.77e-03]    [8.98e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
6000      [1.80e-01, 7.23e-02, 6.70e-03]    [1.90e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
7000      [1.80e-01, 6.88e-02, 6.79e-03]    [1.90e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
8000      [3.66e-01, 8.63e-02, 5.91e-03]    [3.78e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
9000      [1.71e-01, 6.24e-02, 4.21e-03]    [1.80e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
10000     [1.56e-01, 5.24e-02, 1.32e-03]    [1.64e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
11000     [1.80e-01, 3.73e-02, 2.93e-04]    [1.93e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
12000     [1.18e+00, 4.68e-02, 8.33e-03]    [1.25e+00, 0.00e+00, 0.00e+00]    [1.00e+00]    
13000     [1.64e-01, 1.72e-02, 7.67e-03]    [1.72e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
14000     [9.79e-02, 1.72e-02, 1.30e-02]    [1.05e-01, 0.00e+00, 0.00e+00]    [1.00e+00]    
15000     [1.45e-01, 1.33e-02, 1.40e-02]    [1.52e-01, 0.00e+00, 0.00e+00]    [1.00e+00]  

[lululxvi]

There is some problem in the reference solution you provided. -v_s[0] * rho * c * np.cos(k * x) does not satisfy the left BC.

[voloddia]

Actually your train loss is too large. There are 3 columns in train loss: the first column is PDE loss, the second column is left BC loss, and the third column is right BC loss. So the main issue is that the PDE loss doesn't go down.

The possible reason is that your PDE solution is not of scale O(1), which makes the training hard. There are two solutions:

1. re-scale your PDE, such that the solution is O(1);

2. If you assume the PDE solution is of scale O(10), then multiply the network output by 10:

net.outputs_modify(lambda x, y: 10 * y))

I'm just wondering

• Why is the training hard when the solution is not O(1)?
• apply_output_transform(lambda y: 10*y) is just a linear transformation, why can't this be learned by the last layer?
• If the solution has log-uniform distribution, how to make it O(1)?

[lululxvi]

For example, see https://machinelearningmastery.com/how-to-improve-neural-network-stability-and-modeling-performance-with-data-scaling/

If you do exp on log-uniform, then is uniform.

[voloddia]

For example, see https://machinelearningmastery.com/how-to-improve-neural-network-stability-and-modeling-performance-with-data-scaling/

If you do exp on log-uniform, then is uniform.

Thanks, Lu. I've read that article.
The point it made for scaling output is: "You must ensure that the scale of your output variable matches the scale of the activation function (transfer function) on the output layer of your network."
I find this not persuasive because:

• The output layer often has no activation function (in other words, identity activation function) so it has arbitrary scale.
• If I use swish instead of tanh or sigmoid, then even outputs of hidden layers are not upper-bounded.