meliao / helmholtz_solvers

Pytorch implementations of high-resolution solvers for forward wave scattering problems.

helmholtz_solvers

Pytorch implementations of high-resolution solvers for forward scattering problems. The intention is to make this code efficient for rapidly solving a large number of PDEs by taking advantage of hardware accelerators (GPUs).

PDE Problem

Given a compactly-supported two-dimensional scattering potential $q(x)$, the code solves this problem:

$$\Delta u(x) + \omega^2 (1 - q(x)) u(x) = 0 \ \ \ \ x \in \mathbb{R}^2$$

where the solution $u(x) = u_s(x) + u_i(x)$. The incident part of the wave field is a plane wave with known direction and frequency $\omega$. The scattered part of the wave field satisfies the Sommerfeld radiation condition:

$$ \frac{\partial u_s(x)}{\partial | x | } - i \omega u_s(x) = o\left( | x |^{-1/2} \right) \ \ \ | x | \to \infty $$

Solver 1: Lippmann-Schwinger Equation (Stable)

Discretizes and solves the following integral equation:

$$ \sigma(x) - \omega^2(x) q(x) \int G_\omega(x, x') \sigma(x') dx' = \omega^2 u_i(x) q(x)$$

where $G_\omega(x,x')$ is the Green's function for the homogeneous Helmholtz operator $\Delta + \omega^2$. The solution $u_s(x)$ to the above problem is given by $u_s(x) = \int G_\omega(x, x') \sigma(x') dx'$.

The linear equations are discretized and implemented as sparse linear operators, and the system is inverted using an iterative solver GMRES or (optionally) BICGSTAB.

The object HelmholtzSolverAccelerated in the file src/lippmann_schwinger_eqn/HelmholtzSolver.py serves as the main interface.

The object should be initialized using the setup_accelerated_solver() method. Some of the methods that may be useful:

• HemholtzSolverAccelerated.Helmholtz_solve_exterior(): Given incident source directions and a scattering object, evalutates the scattered wave field $u_s(x)$ on a ring far away from the scattering object. This is often used as a forward model in inverse wave scattering problems.
• HelmholtzSolverAccelerated.Helmholtz_solve_interior(): Given incident source directions and a scattering object, evaluates the solution $u(x) = u_s(x) + u_i(x)$ on the scattering domain. Returns the total wave field, the incident wave field, and the scattered wave field.

Solver 2: Hierarchical Poincare-Steklov (Under Development)

The HPS method uses a heirarchical spatial decomposition to define and directly solve a linear system defined on the boundary of the scattering domain.

Reading

These are the papers I've read about the HPS method:

• (BGM15) A spectrally accurate direct solution technique for frequency-domain scattering problems with variable media
• (M15) The Hierarchical Poincare-Steklov (HPS) solver for elliptic PDEs: A tutorial

Other references:

• (T00) Spectral Methods in MATLAB
• (S11) Numerical Analysis

TODO

• Build and test 1D Gauss-Legendre quadrature object.
• Re-write 2D Cheby quadrature object to include a list of indices mapping from 1d point locations to indices in the 2d point array.
• Write tests for the 2D Cheby quadrature object
• Test the 1D differentiation matrix object.
• Write code to interpolate from equispaced grid to Chebyshev grid.
• Test code to interpolate from equispaced to Chebyshev grid.
• Build the LeafNode object.
• Test the D_x and D_y operators by differentiating polynomials.
• Test the LeafNode object by building a notebook and confirming that it produces acceptable local solutions to the scattering problem.
• Build objects for merging LeafNodes
• Read about BIE methods, single- and double-layer potetntials.
• Implement S, the single-layer potential from BGM15 eqn 3.2.
• Implement D, the double-layer potential from BGM15 eqn 3.2.
• Implement T_ext, the exterior DtN map.
• Implement the solve operation given T_int and T_ext.
• Check that things look reasonable on small problem instances.
• Refactor LeafNode, Node, and Merge
  • Node should have objects BoundaryQuad, R, is_leaf_bool, S (optionally, if not a Leaf).
  • Node should have methods get_R_submatrices(),
  • Merge should operate on two Node instances. It should have some method that returns a new Node instance.
• Make a BoundaryQuad object that has the following capabilities:
  • Need a function for merging 2 such quadrature objects together.
  • Is able to give quadrature weights.
  • Handles the submatrix indexing that is currently handled by LeafNode.
• Write a function that builds a tree of nodes.
• Write a function that maps q to T_int.

TODO: Optimization

• Precompute an interpolation matrix from a regularly-spaced grid to a Chebyshev grid. The code in Cheby2D.interp_to_2d_points() relies on scipy and does not take advantage of any precomputation.
• Factor the stuff in LeafNode.__init__() to be all pre-computed and have the diagonals of certain objects be updated given a new scattering object or frequency.
• Optimize the linear system solve in LeafNode.solve()
• There are multiple steps in LeafNode that require taking the Kronecker product between a dense matrix and 4x4 identity matrix. Can we accelerate applying these operators?