Helmholtz Equation and PointSetBC Problem

[JakobEliasWagner]

Hi, @lululxvi I am pretty new to this framework and got a problem many people seem to be getting. I love DeepXDE and am trying to use it to deal with a acoustic Helmholtz problem. I created training data with FEniCS. I set up the training points only with Anchors and a PointSetBC. I prepared the data to be O(1). The problem is that my network does not seem to converge and not even display the trivial solution of the Helmholtz problem. Iused FEniCS in complex-mode therefore there is a complex-valued solution, but only used the real part to train the network. I used the following architecture:

• [1] + [50] * 5 +[1]
• trained with "adam"
• for 1e+4 epochs I got the following results
• 10000 [2.53e-07, 3.40e-01] [2.53e-07, 0.00e+00] []

Because of that i set up a Poisson Equation problem to understand the problem in a better way and eliminate the issue with a complex part of the solution. It did work a little bit better but did not achive any kind of accuracy either. Thus i must conclude there is another issue concerning my code. I tinkered with the methods apply_feature_transform and apply_output_transform as well. Both did not yield a solution i could interprete beeing better or wors in comparison to my first. What is wrong with my approach? What could I be tweaking? To not clutter this issue with code i uploaded my code to this repo.

Many thanks in advance, Jakob

[lululxvi]

There are many Poisson equations in the examples: https://github.com/lululxvi/deepxde/tree/master/examples. You may have a check.

[JakobEliasWagner]

Thank you for your answer @lululxvi. I've looked at the examples again and my network still is not converging as i would expect it to be doing. Here are some plots to illustrate my point. This is the best i got after tinkering with some hyperparameters. I can't really make sense of what is happening, as i tried to compare my code to the examples and tried implementing variants shown in those. I tried to stick in my formulation - not implementation - to the example Poisson_Dirichlet_1d.py. Thank you

from __future__ import division
from __future__ import print_function

import pathlib
import numpy as np
import matplotlib.pyplot as plt
import h5py

import deepxde as dde
from deepxde.backend import tf
import sklearn.preprocessing

def load_test_data(path_to_file: pathlib.Path):
    """
    :param path_to_file:
    :return: x - (collocation_points, 2) y - (collocation_points, 1)
    """
    with h5py.File(path_to_file, 'r') as file:
        x = file['Mesh']['mesh']['geometry'][()]
        key = [key for key in file['Function'].keys()]
        y = file['Function'][key[0]]['0'][()]
    return x[:, :1], x[:, 1:], y

def pde(x, y):
    """
    :param x: x-value
    :param y: sound pressure
    :return:
    """
    dy_xx = dde.grad.hessian(y, x, i=0, j=0)
    return -(-dy_xx - np.pi ** 2 * tf.sin(np.pi * x))

def sol(phi):
    return np.sin(np.pi * phi)

if __name__ == "__main__":
    FENICS_DATA = False

    if FENICS_DATA:
        curr_file_path = pathlib.Path(__file__)
        h5_file_path = curr_file_path.parent.joinpath('../../test_data/sol_poisson.h5')
        x_t, y_t, u_t = load_test_data(h5_file_path)
    else:
        x_t = np.linspace(-1, 1, 100)[:, None]
        u_t = sol(x_t)

    # prepare training data
    scaler = sklearn.preprocessing.MinMaxScaler()
    u_t = scaler.fit_transform(u_t)

    geom = dde.geometry.Interval(-1, 1)
    observe_y0 = dde.PointSetBC(x_t, u_t)

    bcs = [observe_y0]
    data = dde.data.PDE(geom, pde, bcs, num_domain=400, num_boundary=2, anchors=x_t)

    net = dde.maps.FNN([1] + [50] * 3 + [1], "tanh", "Glorot normal")

    model = dde.Model(data, net)
    model.compile("adam", lr=1e-3, loss_weights=[1, 100])

    h, t = model.train(epochs=int(1e+4))

    dde.saveplot(h, t, issave=False, isplot=True)

    x = geom.uniform_points(1000, True)
    y = model.predict(x, operator=pde)

    y = scaler.inverse_transform(y)
    u_t = scaler.inverse_transform(u_t)

    plt.figure()
    plt.plot(x, y, 'b')
    plt.plot(x_t, u_t, 'r*-')
    plt.show()

[lululxvi]

Remove the scaler.

[JakobEliasWagner]

Hey @lululxvi, with the scaler removed it looks like this: Thank you for your answer, Jakob

[lululxvi]

I am not sure why your code doesn't work. This is just a simple Poission eq. It is strange.

[JakobEliasWagner]

Okay, thank you anyways for looking at my code :)

[louhz]

I guess that maybe you can use the Poincaré–Steklov operator to propose an artificial boundary for this wave function. By looking into your graph, it looks like your predicted result is not in the correct boundary. There are two papers that may be helpful for you: On the implementation of the Dirichlet-to-Neumann radiation condition for iterative solution of the Helmholtz equation by Assad A. Oberai *, Manish Malhotra, Peter M. Pinsky and NON-REFLECTING BOUNDARY CONDITIONS FOR ELASTIC WAVES by Dan GIVOLI, Joseph B. KELLER