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).
Given a compactly-supported two-dimensional scattering potential
where the solution
Discretizes and solves the following integral equation:
where
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.
The HPS method uses a heirarchical spatial decomposition to define and directly solve a linear system defined on the boundary of the scattering domain.
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
- 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
LeafNode
s - 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
, andMerge
-
Node
should have objectsBoundaryQuad
,R
,is_leaf_bool
,S
(optionally, if not a Leaf). -
Node
should have methodsget_R_submatrices()
, -
Merge
should operate on twoNode
instances. It should have some method that returns a newNode
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.
- 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?