Commit c0aa297c authored by Leon Merten Lohse's avatar Leon Merten Lohse
Browse files

optimize

parent 3a5edabd
......@@ -2,15 +2,23 @@ cmake_minimum_required (VERSION 3.11)
project (fresnel VERSION 0.1)
set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD_REQUIRED True)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(default_build_type "Release")
set(PYBIND11_PYTHON_VERSION 3.6)
add_subdirectory (extern/pybind11)
find_package (Eigen3 3.3 REQUIRED NO_MODULE)
find_package(OpenMP)
pybind11_add_module (_solver source/python.cpp)
target_link_libraries (_solver PRIVATE Eigen3::Eigen)
target_compile_options(_solver PRIVATE -Wpedantic -Wall -O3 -ffast-math)
pybind11_add_module (_propagate source/python.cpp)
#target_compile_definitions (_propagate PRIVATE VERSION_INFO=${FRESNEL_VERSION_INFO})
target_link_libraries (_propagate PUBLIC Eigen3::Eigen)
if(OpenMP_CXX_FOUND)
target_link_libraries(_solver PRIVATE OpenMP::OpenMP_CXX)
endif()
install(TARGETS _propagate DESTINATION fresnel)
install(TARGETS _solver DESTINATION fresnel)
install(DIRECTORY fresnel/ DESTINATION fresnel
FILES_MATCHING PATTERN "*.py")
This diff is collapsed.
This diff is collapsed.
import numpy as np
import _propagator
import _solver
_k = 2 * np.pi # all lengths are in units of the wavelength
class Propagator2d():
class Solver2d():
"""
Solver for 2d PDEs of the form
u_z = A * u_xx + F * u
"""
ndim = 1
dtype = np.complex128
dtype = np.complex128
def __init__(self, n0, u0, dz, dx):
# TODO: validate input
def __init__(self, A0, F0, u0, dz, dx):
self._nx = u0.shape[-1]
self._ones = np.ones((self._nx,), dtype = self.dtype)
......@@ -27,24 +23,22 @@ class Propagator2d():
self.dx = dx
self.dz = dz
self._compute_coefficients(n0)
self._compute_coefficients(A0, F0)
def _compute_coefficients(self, n):
A = -1j / (2. * _k)
F = -1j * _k * (n*n - 1) / 2
def _compute_coefficients(self, A, F):
self.rx = A * self.dz / ( 2. * self.dx * self.dx) * self._ones
self.rx = A * self.dz / ( 2. * self.dx**2 ) * self._ones
self.rc = F * self.dz / 2. * self._ones
def step(self, n, boundary):
def step(self, A, F, boundary):
up = self.u
rxp = self.rx
rcp = self.rc
self._compute_coefficients(n)
self._compute_coefficients(A, F)
u = self._boundary_slice
rx = self.rx
......@@ -53,14 +47,39 @@ class Propagator2d():
# set boundary values
u[0] = boundary[0]
u[-1] = boundary[1]
self.u = _fresnel.step_AF(up, rxp, rcp, u, rx, rc)
self.u = _solver.step_AF(up, rxp, rcp, u, rx, rc)
return self.u
class Propagator2d(Solver2d):
def __init__(self, n0, u0, dz, dx):
class Propagator3d():
A, F = self._compute_helmholtz_coefficients(n0)
super().__init__(A, F, u0, dz, dx)
def _compute_helmholtz_coefficients(self,n):
A = -1j / (2. * _k)
F = -1j * _k * (n*n - 1) / 2
return A, F
def step(self, n, boundary):
A, F = self._compute_helmholtz_coefficients(n)
return super().step(A, F, boundary)
class Solver3d():
"""
Solver for equations of the form
u_z = A * u_xx + C * u_yy + F * u
......@@ -71,15 +90,11 @@ class Propagator3d():
'rf : F * dz / 2
"""
ndim = 2
dtype = np.complex128
def __init__(self, n0, u0, dz, dy, dx):
def __init__(self, A0, F0, u0, dz, dy, dx):
# TODO: validate input
self._nx = u0.shape[-1]
self._ny = u0.shape[-2]
self._ones = np.ones((self._ny, self._nx,), dtype = self.dtype)
......@@ -91,26 +106,24 @@ class Propagator3d():
self.dy = dy
self.dz = dz
self._compute_coefficients(n0)
self._compute_coefficients(A0, F0)
def _compute_coefficients(self, n):
A = -1j / (2. * _k)
F = -1j * _k * (n*n - 1) / 2
def _compute_coefficients(self, A, F):
self.rx = A * self.dz / ( 2. * self.dx * self.dx) * self._ones
self.ry = A * self.dz / ( 2. * self.dy * self.dy) * self._ones
self.rc = F * self.dz / 4. * self._ones
def step(self, n, boundary):
def step(self, A, F, boundary):
up = self.u
rxp = self.rx
ryp = self.ry
rcp = self.rc
self._compute_coefficients(n)
self._compute_coefficients(A, F)
u = self._boundary_slice
rx = self.rx
......@@ -123,6 +136,30 @@ class Propagator3d():
u[:, 0] = boundary[2]
u[:,-1] = boundary[3]
self.u = _fresnel.step_ACF(up, rxp, ryp, rcp, u, rx, ry, rc)
self.u = _solver.step_ACF(up, rxp, ryp, rcp, u, rx, ry, rc)
return self.u
class Propagator3d(Solver3d):
def __init__(self, n0, u0, dz, dy, dx):
A, F = self._compute_helmholtz_coefficients(n0)
super().__init__(A, F, u0, dz, dy, dx)
def _compute_helmholtz_coefficients(self,n):
A = -1j / (2. * _k)
F = -1j * _k * (n*n - 1) / 2
return A, F
def step(self, n, boundary):
A, F = self._compute_helmholtz_coefficients(n)
return self.u
\ No newline at end of file
return super().step(A, F, boundary)
\ No newline at end of file
......@@ -143,6 +143,7 @@ namespace finite_differences{
array_2D uhalf_t = array_2D::Zero(u.cols(), u.rows()); // transposed shape
#pragma omp parallel for
for (size_t iy = 1; iy <= ny; ++iy) {
algebra::TriMatrix<scalar> m = algebra::TriMatrix<scalar>::Zero(3, nx);
......@@ -170,6 +171,7 @@ namespace finite_differences{
array_2D u_new { u };
#pragma omp parallel for
for (size_t ix = 1; ix <= nx; ++ix) {
algebra::TriMatrix<scalar> m = algebra::TriMatrix<scalar>::Zero(3, ny);
......
......@@ -8,7 +8,7 @@ namespace py = pybind11;
using namespace finite_differences;
PYBIND11_MODULE(_fresnel, m){
PYBIND11_MODULE(_solver, m){
m.def("step_AF",
&finite_differences::step_AF
......@@ -17,56 +17,5 @@ PYBIND11_MODULE(_fresnel, m){
m.def("step_ACF",
&finite_differences::step_ACF
);
/*
py::class_<finite_difference_AF>(m, "finite_difference_AF")
.def(py::init<size_t>())
.def("step", &finite_difference_AF::step)
.def("update", &finite_difference_AF::update)
.def("get_field", &finite_difference_AF::get_field, py::return_value_policy::reference_internal)
.def("set_field", &finite_difference_AF::set_field)
.def("set_rx", &finite_difference_AF::set_rx)
.def("set_rc", &finite_difference_AF::set_rc)
;
py::class_<finite_difference_ACF>(m, "finite_difference_ACF")
.def_readwrite("thread_count",&finite_difference_ACF::thread_count)
.def(py::init<size_t,size_t>())
.def("step_1", &finite_difference_ACF::step_1)
.def("step_2", &finite_difference_ACF::step_2)
.def("update", &finite_difference_ACF::update)
.def("get_field", &finite_difference_ACF::get_field, py::return_value_policy::reference_internal)
.def("set_field", &finite_difference_ACF::set_field)
.def("set_ra", &finite_difference_ACF::set_ra)
.def("set_rc", &finite_difference_ACF::set_rc)
.def("set_rf", &finite_difference_ACF::set_rf)
;
py::class_<finite_difference_A0F>(m, "finite_difference_A0F")
.def_readwrite("thread_count",&finite_difference_A0F::thread_count)
.def(py::init<size_t, size_t>())
.def("step", &finite_difference_A0F::step)
.def("update", &finite_difference_A0F::update)
.def("get_field", &finite_difference_A0F::get_field, py::return_value_policy::reference_internal)
.def("set_field", &finite_difference_A0F::set_field)
.def("set_ra", &finite_difference_A0F::set_ra)
.def("set_rf", &finite_difference_A0F::set_rf)
;
py::class_<finite_difference_ABC>(m, "finite_difference_ABC")
.def_readwrite("thread_count",&finite_difference_ABC::thread_count)
.def(py::init<size_t,size_t>())
.def("step", &finite_difference_ABC::step)
.def("update", &finite_difference_ABC::update)
.def("get_field", &finite_difference_ABC::get_field, py::return_value_policy::reference_internal)
.def("set_field", &finite_difference_ABC::set_field)
.def("set_ra", &finite_difference_ABC::set_ra)
.def("set_rb", &finite_difference_ABC::set_rb)
.def("set_rc", &finite_difference_ABC::set_rc)
.def("set_rz", &finite_difference_ABC::set_rz)
;
*/
}
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