Commit 9b4a3bfe authored by lohse6's avatar lohse6
Browse files

Merge branch 'pml' into 'master'

Pml

See merge request !2
parents 0a1d9c95 07fbd59c
Pipeline #235841 canceled with stage
This diff is collapsed.
This diff is collapsed.
......@@ -3,15 +3,27 @@ import numpy as np
from . import _solver
def _compute_pml(width, pml_width, dtype=np.complex128):
sigma_rel = np.zeros(width, dtype=dtype)
if pml_width > 0:
sigma_ramp = np.linspace(1, 0, pml_width, endpoint=False)
sigma_rel[:pml_width] = sigma_ramp
sigma_rel[-pml_width:] = sigma_ramp[::-1]
return sigma_rel ** 2
class Solver2d:
"""
Solver for 2d PDEs of the form
A * u_z = -u_xx + F(x,z) * u
with perfectly matched layers
"""
dtype = np.complex128
def __init__(self, A, F0, u0, dz, dx):
def __init__(self, A, F0, u0, dz, dx, pml_width=0, sigma_max=0.1):
self._nx = u0.shape[-1]
self._ones = np.ones((self._nx,), dtype=self.dtype)
......@@ -24,6 +36,9 @@ class Solver2d:
self.a = self._compute_a(A)
self.f = self._compute_f(F0)
sigma_rel = _compute_pml(self._nx, pml_width, self.dtype)
self.eta = 1.0 / (1.0 + 1j * sigma_max * sigma_rel)
def _compute_a(self, A):
return A * self.dz / self.dx ** 2 * (1 + 0j)
......@@ -37,7 +52,7 @@ class Solver2d:
self.f = self._compute_f(F)
self.u = _solver.step1d_A0F(self.a, fp, self.f, up, boundary[0], boundary[1])
self.u = _solver.step1d_A0F_pml(self.a, fp, self.f, up, boundary[0], boundary[1], self.eta)
return self.u
......
......@@ -55,13 +55,14 @@ class MultislicePropagator:
class FDPropagator2d(fd.Solver2d):
def __init__(self, n0, u0, dz, dx, wl=1.0):
def __init__(self, n0, u0, dz, dx, pml_width=0, sigma_max=0.1, wl=1.0):
self._k = _k0 / wl
A = 1j / (2 * self._k)
F0 = _compute_potential(n0, self._k)
sigma_max_rel = sigma_max * wl
super().__init__(A, F0, u0, dz, dx)
super().__init__(A, F0, u0, dz, dx, pml_width, sigma_max_rel)
def step(self, n, boundary):
......
......@@ -117,6 +117,54 @@ namespace finite_differences {
return u_new;
}
const array_1D step1d_A0F_pml(const scalar a,
const Eigen::Ref<const array_1D> fp,
const Eigen::Ref<const array_1D> f,
const Eigen::Ref<const array_1D> up,
const scalar bl,
const scalar br,
const Eigen::Ref<const array_1D> eta) {
// _p : i
// _ : i+1
size_t n = up.size()-2;
// assert eta.size == n + 2
// setup tridiagonal n x n matrix M = tridiag(a,b,c) and rhs r
algebra::TriMatrix<scalar> m(3, n);
array_1D r(n, 1);
// Note that matrix indices go from 1 to up.size()-1, so that indices are off by 1
for (size_t i = 1; i <= n; ++i) {
m(0, i-1) = - a * eta(i) * (0.25 * eta(i-1) + eta(i) - 0.25 * eta(i+1)); // lower
m(1, i-1) = 2.0 + 2.0 * a * eta(i) * eta(i) - f(i); // diagonal
m(2, i-1) = - a * eta(i) * (-0.25 * eta(i-1) + eta(i) + 0.25 * eta(i+1)); // upper;
r(i-1) = a * eta(i) * (0.25 * eta(i-1) + eta(i) - 0.25 * eta(i+1)) * up(i-1)
+ (2.0 - 2.0 * a * eta(i) * eta(i) + fp(i)) * up(i)
+ a * eta(i) * (-0.25 * eta(i-1) + eta(i) + 0.25 * eta(i+1)) * up(i+1);
}
m(0, 0) = 0;
m(2, n-1) = 0;
// boundary
r(0) += a * bl * eta(1) * (0.25 * eta(0) + eta(1) - 0.25 * eta(2));
r(n-1) += a * br * eta(n) * (-0.25 * eta(n-1) + eta(n) + 0.25 * eta(n+1));
// initialize writeable copy
array_1D u_new(up.size(), 1);
// insert boundary values
u_new(0) = bl;
u_new(n+1) = br;
u_new.segment(1, n) = algebra::tridiagonal(std::move(m), std::move(r));
return u_new;
}
const array_1D step1d_AAF_sym(const complex a,
const Eigen::Ref<const array_1D> fp,
......
......@@ -17,6 +17,17 @@ PYBIND11_MODULE(_solver, m) {
py::arg("br").noconvert(),
py::return_value_policy::move);
m.def("step1d_A0F_pml",
&fd::step1d_A0F_pml,
py::arg("a").noconvert(),
py::arg("fp").noconvert(),
py::arg("f").noconvert(),
py::arg("up").noconvert(),
py::arg("bl").noconvert(),
py::arg("br").noconvert(),
py::arg("eta").noconvert(),
py::return_value_policy::move);
m.def("step2d_A0F",
&fd::step2d_A0F,
py::arg("ax").noconvert(),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment