diff --git a/.github/workflows/codeql-analysis.yml b/.github/workflows/codeql-analysis.yml index a35d47d02..8272d1d36 100644 --- a/.github/workflows/codeql-analysis.yml +++ b/.github/workflows/codeql-analysis.yml @@ -41,7 +41,7 @@ jobs: # ✏️ If the Autobuild fails above, remove it and uncomment the following three lines # and modify them (or add more) to build your code if your project # uses a compiled language - + - name: Build run: | cmake -S multi_physics/QED \ diff --git a/.github/workflows/dependencies/dependencies.sh b/.github/workflows/dependencies/dependencies.sh index c1ec45e53..6b52d7e83 100755 --- a/.github/workflows/dependencies/dependencies.sh +++ b/.github/workflows/dependencies/dependencies.sh @@ -3,7 +3,7 @@ # Copyright 2020 The PICSAR Community # # License: BSD-3-Clause-LBNL -# Authors: Axel Huebl +# Authors: Axel Huebl, Luca Fedeli set -eu -o pipefail @@ -15,4 +15,7 @@ sudo apt-get install -y --no-install-recommends \ libboost-dev \ libboost-math-dev \ libboost-test-dev \ - g++ gfortran + g++ gfortran \ + python3-dev + +pip install "pybind11[global]" diff --git a/.github/workflows/dependencies/dependencies_mac.sh b/.github/workflows/dependencies/dependencies_mac.sh index a7e9de985..96b4b984d 100755 --- a/.github/workflows/dependencies/dependencies_mac.sh +++ b/.github/workflows/dependencies/dependencies_mac.sh @@ -1,13 +1,14 @@ #!/usr/bin/env bash # -# Copyright 2020 The AMReX Community +# Copyright 2020 The PICSAR Community # # License: BSD-3-Clause-LBNL -# Authors: Axel Huebl +# Authors: Axel Huebl, Luca Fedeli set -eu -o pipefail brew update brew install boost brew install libomp +brew install pybind11 #brew install open-mpi diff --git a/multi_physics/QED/CMakeLists.txt b/multi_physics/QED/CMakeLists.txt index 11f996466..f47ab7cb1 100644 --- a/multi_physics/QED/CMakeLists.txt +++ b/multi_physics/QED/CMakeLists.txt @@ -20,6 +20,7 @@ option(PXRMP_QED_OMP "Enable OpenMP support" ON) option(PXRMP_QED_TABLEGEN "Enable table generation (needs Boost)" ON) option(PXRMP_QED_TEST "Build PICSAR QED tests " ${IS_TOPLEVEL}) option(PXRMP_QED_TOOLS "Build PICSAR QED tools " ${IS_TOPLEVEL}) +option(PXRMP_QED_PYTHON_BINDINGS "Build PICSAR QED python bindings " ${IS_TOPLEVEL}) option(PXRMP_BOOST_TEST_DYN_LINK "Link against the Boost Unit Test Framework shared library" ON) option(PXRMP_DPCPP_FIX "Use cl::sycl::floor cl::sycl::floorf on device" OFF) @@ -50,6 +51,10 @@ if(PXRMP_QED_TABLEGEN) endif() endif() +if(PXRMP_QED_PYTHON_BINDINGS) + add_subdirectory(python_bindings) +endif() + if(PXRMP_DPCPP_FIX) target_compile_definitions(PXRMP_QED INTERFACE PXRMP_DPCPP_FIX=1) endif() diff --git a/multi_physics/QED/python_bindings/CMakeLists.txt b/multi_physics/QED/python_bindings/CMakeLists.txt new file mode 100644 index 000000000..96cc74b76 --- /dev/null +++ b/multi_physics/QED/python_bindings/CMakeLists.txt @@ -0,0 +1,71 @@ +set(name pxr_qed) + +set(_PY_DEV_MODULE Development.Module) +if(CMAKE_VERSION VERSION_LESS 3.18.0) + # over-specification needed for CMake<3.18 + # https://pybind11.readthedocs.io/en/latest/compiling.html#findpython-mode + # https://cmake.org/cmake/help/v3.18/module/FindPython.html + set(_PY_DEV_MODULE Development) +endif() +find_package(Python COMPONENTS Interpreter ${_PY_DEV_MODULE} REQUIRED) +find_package(pybind11 2.6 CONFIG REQUIRED) + +pybind11_add_module(${name} "${name}.cpp") + +target_link_libraries(${name} PRIVATE PXRMP_QED) + +option( + PXRQEDPY_DOUBLE_PRECISION + "Compile pxr_qed python module in double precision" ON) + +if(PXRQEDPY_DOUBLE_PRECISION) + target_compile_definitions(${name} PRIVATE PXRQEDPY_DOUBLE_PRECISION=1) +endif() + +set(PXRQEDPY_UNITS_VALUES SI NORM_OMEGA NORM_LAMBDA HEAVISIDE_LORENTZ) +set(PXRQEDPY_UNITS SI CACHE STRING "Unit system for the pxr_qed python module (SI/NORM_OMEGA/NORM_LAMBDA/HEAVISIDE_LORENTZ)") +set_property(CACHE PXRQEDPY_UNITS PROPERTY STRINGS ${PXRQEDPY_UNITS_VALUES}) + +if(NOT PXRQEDPY_UNITS IN_LIST PXRQEDPY_UNITS_VALUES) + message(FATAL_ERROR "PXRQEDPY_UNITS (${PXRQEDPY_UNITS}) must be one of ${PXRQEDPY_UNITS_VALUES}") +endif() + +if(PXRQEDPY_UNITS STREQUAL "SI") + target_compile_definitions(${name} PRIVATE PXRQEDPY_SI=1) +elseif(PXRQEDPY_UNITS STREQUAL "NORM_OMEGA") + target_compile_definitions(${name} PRIVATE PXRQEDPY_NORM_OMEGA=1) +elseif(PXRQEDPY_UNITS STREQUAL "NORM_LAMBDA") + target_compile_definitions(${name} PRIVATE PXRQEDPY_NORM_LAMBDA=1) +elseif(PXRQEDPY_UNITS STREQUAL "HEAVISIDE_LORENTZ") + target_compile_definitions(${name} PRIVATE PXRQEDPY_HEAVISIDE_LORENTZ=1) +else() + message(FATAL_ERROR "PXRQEDPY_UNITS is incorrectly defined") +endif() + +target_link_libraries(${name} PRIVATE pybind11::module pybind11::lto pybind11::windows_extras) + +# OpenMP support +if(PXRMP_QED_OMP) + find_package(OpenMP REQUIRED) + target_link_libraries(${name} PRIVATE OpenMP::OpenMP_CXX) + target_compile_definitions(${name} PRIVATE PXQEDPY_HAS_OPENMP=1) +endif() + +# Move module in "python_bindings" subdirectory +set(PXRQEDPY_INSTALL_DIR ${CMAKE_BINARY_DIR}/python_bindings CACHE PATH "Installation directory for the python bindings") +set_target_properties(${name} PROPERTIES + RUNTIME_OUTPUT_DIRECTORY ${PXRQEDPY_INSTALL_DIR}) + +configure_file(demo_python_bindings.ipynb + ${CMAKE_BINARY_DIR}/python_bindings/demo_python_bindings.ipynb COPYONLY) + +# Require C++14 or newer +target_compile_features(${name} PUBLIC cxx_std_14) +set_target_properties(${name} PROPERTIES CXX_EXTENSIONS OFF) + +# Enable warnings +if(MSVC) + target_compile_options(${name} PRIVATE /W4) +else() + target_compile_options(${name} PRIVATE -Wall -Wextra -pedantic) +endif() diff --git a/multi_physics/QED/python_bindings/demo_python_bindings.ipynb b/multi_physics/QED/python_bindings/demo_python_bindings.ipynb new file mode 100644 index 000000000..dc839594e --- /dev/null +++ b/multi_physics/QED/python_bindings/demo_python_bindings.ipynb @@ -0,0 +1,473 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "greatest-rally", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "pxr_qed is compiled in double precision with SI units\n", + "pxr_qed has openMP support!\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import scipy as sci\n", + "import scipy.constants as scicon\n", + "import pxr_qed\n", + "\n", + "print(\"pxr_qed is compiled in {:s} precision with {:s} units\".\n", + " format(pxr_qed.PRECISION, pxr_qed.UNITS))\n", + "if(pxr_qed.HAS_OPENMP):\n", + " print(\"pxr_qed has openMP support!\")\n", + "else:\n", + " print(\"pxr_qed does not have openMP support!\")\n", + "\n", + "#Warning: this notebook is conceived for SI units!!" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "parliamentary-installation", + "metadata": {}, + "outputs": [], + "source": [ + "m_e = scicon.electron_mass\n", + "q_e = scicon.elementary_charge\n", + "hbar = scicon.hbar\n", + "c = scicon.c\n", + "um = scicon.micron\n", + "fs = scicon.femto\n", + "\n", + "E_s = m_e**2 * c**3 / (hbar * q_e)\n", + "B_s = E_s/c" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "attended-stock", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "chi_phot: [ 5.69987791 5.82570482 10.97872065 5.89781209 14.09537311 1.19244247\n", + " 6.74005305 6.19747655 11.90240205 16.65693906]\n", + "chi_ele_pos: [ 5.74174727 5.80739474 11.002853 5.97281582 14.13122611 1.34429715\n", + " 6.78460622 6.22443299 11.95747209 16.69832577]\n" + ] + } + ], + "source": [ + "#Chi functions\n", + "\n", + "chi_how_many = 10\n", + "chi_gamma = 10\n", + "\n", + "chi_Ex = (np.random.rand(chi_how_many)-0.5)*2*E_s\n", + "chi_Ey = (np.random.rand(chi_how_many)-0.5)*2*E_s\n", + "chi_Ez = (np.random.rand(chi_how_many)-0.5)*2*E_s\n", + "chi_Bx = (np.random.rand(chi_how_many)-0.5)*2*B_s\n", + "chi_By = (np.random.rand(chi_how_many)-0.5)*2*B_s\n", + "chi_Bz = (np.random.rand(chi_how_many)-0.5)*2*B_s\n", + "chi_px = (np.random.rand(chi_how_many)-0.5)*2*chi_gamma*m_e*c;\n", + "chi_py = (np.random.rand(chi_how_many)-0.5)*2*chi_gamma*m_e*c;\n", + "chi_pz = (np.random.rand(chi_how_many)-0.5)*2*chi_gamma*m_e*c;\n", + "\n", + "print(\"chi_phot: \", pxr_qed.chi_photon(\n", + " chi_px, chi_py, chi_pz,\n", + " chi_Ex, chi_Ey, chi_Ez,\n", + " chi_Bx, chi_By, chi_Bz, 1))\n", + "\n", + "print(\"chi_ele_pos: \", pxr_qed.chi_ele_pos(\n", + " chi_px, chi_py, chi_pz, \n", + " chi_Ex, chi_Ey, chi_Ez,\n", + " chi_Bx, chi_By, chi_Bz, 1))" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "indie-screen", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "bw.dndt_lookup_table_params:\n", + "\tchi_phot_min : 0.01\n", + "\tchi_phot_max : 100\n", + "\tchi_phot_how_many: 128\n", + "bw.dndt_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "bw.dndt_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "bw.dndt_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "bw.dndt_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "bw.pair_prod_lookup_table_params:\n", + "\tchi_phot_min : 0.01\n", + "\tchi_phot_max : 100\n", + "\tchi_phot_how_many: 64\n", + "\tfrac_how_many : 64\n", + "bw.pair_prod_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "bw.pair_prod_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "bw.pair_prod_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "bw.pair_prod_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "Chi parameters: [ 4.52085808 11.13066108 14.24967391 4.23136133] \n", + "\n", + "dN/dt lookup table interpolation: [0.0920503 0.10880774 0.10828566 0.08944404] \n", + "\n", + "Pair production lookup table interpolation: [ 0.76526949 10.60995375 2.80380622 3.21306961] \n", + "\n", + "Breit-Wheeler optical depth: [0.11317728 3.97356138 0.24272302 1.66933813] \n", + "\n", + "Breit-Wheeler dNdt: [2.66096115e+17 7.73092531e+17 8.33218698e+17 2.92127171e+17] \n", + "\n", + "Breit-Wheeler opt depth before evolution: [0.11317728 3.97356138 0.24272302 1.66933813] \n", + "\n", + "Breit-Wheeler opt depth after evolution: [-2.54778386 -3.75736393 -8.08946396 -1.25193358] \n", + "\n", + "Breit-Wheeler electron px: [-6.67695454e-23 5.44758270e-22 -3.23205737e-22 -1.27335703e-21]\n", + "Breit-Wheeler electron py: [ 3.54747146e-22 1.92691707e-21 -3.41196791e-22 -3.50969928e-22]\n", + "Breit-Wheeler electron pz: [-3.78765620e-22 4.13137730e-22 4.86347765e-22 2.93773362e-22]\n", + "Breit-Wheeler positron px: [-2.30829557e-22 6.28739332e-23 -1.01227706e-21 -5.28018577e-22]\n", + "Breit-Wheeler positron py: [ 1.22639935e-21 2.22397459e-22 -1.06862485e-21 -1.45535492e-22]\n", + "Breit-Wheeler positron pz: [-1.30943381e-21 4.76827898e-23 1.52323622e-21 1.21817989e-22]\n" + ] + } + ], + "source": [ + "#Breit-Wheeler pair production\n", + "\n", + "bw_dndt_params = pxr_qed.bw.dndt_lookup_table_params()\n", + "bw_dndt_params.chi_phot_min = 0.01\n", + "bw_dndt_params.chi_phot_max = 100\n", + "bw_dndt_params.chi_phot_how_many = 128\n", + "print(bw_dndt_params)\n", + "\n", + "bw_dndt_lookup_table = pxr_qed.bw.dndt_lookup_table(bw_dndt_params)\n", + "print(bw_dndt_lookup_table)\n", + "bw_dndt_lookup_table.generate()\n", + "print(bw_dndt_lookup_table)\n", + "\n", + "bw_dndt_lookup_table.save_as(\"bw_dndt.bin\")\n", + "bw_dndt_lookup_table_2 = pxr_qed.bw.dndt_lookup_table()\n", + "print(bw_dndt_lookup_table_2)\n", + "bw_dndt_lookup_table_2.load_from(\"bw_dndt.bin\");\n", + "print(bw_dndt_lookup_table_2)\n", + "\n", + "bw_pair_prod_params = pxr_qed.bw.pair_prod_lookup_table_params()\n", + "bw_pair_prod_params.chi_phot_min = 0.01\n", + "bw_pair_prod_params.chi_phot_max = 100\n", + "bw_pair_prod_params.chi_phot_how_many = 64\n", + "bw_pair_prod_params.frac_how_many = 64\n", + "print(bw_pair_prod_params)\n", + "\n", + "bw_pair_prod_lookup_table = pxr_qed.bw.pair_prod_lookup_table(bw_pair_prod_params)\n", + "print(bw_pair_prod_lookup_table)\n", + "bw_pair_prod_lookup_table.generate()\n", + "print(bw_pair_prod_lookup_table)\n", + "\n", + "bw_pair_prod_lookup_table.save_as(\"bw_pairprod.bin\")\n", + "bw_pair_prod_lookup_table_2 = pxr_qed.bw.pair_prod_lookup_table()\n", + "print(bw_pair_prod_lookup_table_2)\n", + "bw_pair_prod_lookup_table_2.load_from(\"bw_pairprod.bin\");\n", + "print(bw_pair_prod_lookup_table_2)\n", + "\n", + "bw_how_many = 4\n", + "bw_gamma = 10\n", + "bw_dt = 1e-2*fs\n", + "\n", + "bw_Ex = (np.random.rand(bw_how_many)-0.5)*2*E_s\n", + "bw_Ey = (np.random.rand(bw_how_many)-0.5)*2*E_s\n", + "bw_Ez = (np.random.rand(bw_how_many)-0.5)*2*E_s\n", + "bw_Bx = (np.random.rand(bw_how_many)-0.5)*2*B_s\n", + "bw_By = (np.random.rand(bw_how_many)-0.5)*2*B_s\n", + "bw_Bz = (np.random.rand(bw_how_many)-0.5)*2*B_s\n", + "bw_px = (np.random.rand(bw_how_many)-0.5)*2*bw_gamma*m_e*c;\n", + "bw_py = (np.random.rand(bw_how_many)-0.5)*2*bw_gamma*m_e*c;\n", + "bw_pz = (np.random.rand(bw_how_many)-0.5)*2*bw_gamma*m_e*c\n", + "bw_rand = np.random.rand(bw_how_many)\n", + "bw_ee = np.sqrt(bw_px**2 + bw_py**2 + bw_pz**2)*c;\n", + "\n", + "\n", + "bw_chi = pxr_qed.chi_photon(bw_px, bw_py, bw_pz, bw_Ex, bw_Ey, bw_Ez, bw_Bx, bw_By, bw_Bz)\n", + "print(\"Chi parameters: \", bw_chi, \"\\n\")\n", + "\n", + "bw_interp_dndt = bw_dndt_lookup_table.interp(bw_chi)\n", + "print(\"dN/dt lookup table interpolation: \", bw_interp_dndt, \"\\n\")\n", + "\n", + "bw_interp_pair_prod = bw_pair_prod_lookup_table.interp(bw_chi, bw_rand)\n", + "print(\"Pair production lookup table interpolation: \", bw_interp_pair_prod, \"\\n\")\n", + "\n", + "bw_opt = pxr_qed.bw.get_optical_depth(bw_rand)\n", + "print(\"Breit-Wheeler optical depth: \", bw_opt, \"\\n\")\n", + "\n", + "bw_dndt = pxr_qed.bw.get_dn_dt(bw_ee, bw_chi, bw_dndt_lookup_table)\n", + "print(\"Breit-Wheeler dNdt: \", bw_dndt, \"\\n\")\n", + "\n", + "print(\"Breit-Wheeler opt depth before evolution: \", bw_opt, \"\\n\")\n", + "pxr_qed.bw.evolve_optical_depth(bw_ee, bw_chi, bw_dt, bw_opt, bw_dndt_lookup_table)\n", + "print(\"Breit-Wheeler opt depth after evolution: \", bw_opt, \"\\n\")\n", + "\n", + "bw_ele_px, bw_ele_py, bw_ele_pz, bw_pos_px, bw_pos_py, bw_pos_pz = pxr_qed.bw.generate_breit_wheeler_pairs(\n", + " bw_chi, bw_px, bw_py, bw_pz, bw_rand, bw_pair_prod_lookup_table)\n", + "\n", + "print(\"Breit-Wheeler electron px: \", bw_ele_px)\n", + "print(\"Breit-Wheeler electron py: \", bw_ele_py)\n", + "print(\"Breit-Wheeler electron pz: \", bw_ele_pz)\n", + "print(\"Breit-Wheeler positron px: \", bw_pos_px)\n", + "print(\"Breit-Wheeler positron py: \", bw_pos_py)\n", + "print(\"Breit-Wheeler positron pz: \", bw_pos_pz)" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "settled-humor", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "qs.dndt_lookup_table_params:\n", + "\tchi_part_min : 0.01\n", + "\tchi_part_max : 100\n", + "\tchi_part_how_many: 128\n", + "qs.dndt_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "qs.dndt_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "qs.dndt_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "qs.dndt_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "qs.photon_emission_lookup_table_params:\n", + "\tchi_part_min : 0.01\n", + "\tchi_part_max : 100\n", + "\tfrac_min : 1e-12\n", + "\tchi_part_how_many: 64\n", + "\tfrac_how_many : 64\n", + "qs.dndt_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "qs.dndt_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "qs.pair_prod_lookup_table:\n", + "\tis initialized? : False\n", + "\n", + "qs.pair_prod_lookup_table:\n", + "\tis initialized? : True\n", + "\n", + "Chi parameters: [14.31260562 18.67401583 19.20375739 10.46322055] \n", + "\n", + "dN/dt lookup table interpolation: [11.87633168 14.35585401 14.64345668 9.47657392] \n", + "\n", + "Pair production lookup table interpolation: [10.62380866 0.0798416 1.49208283 0.07908836] \n", + "\n", + "Quantum synchrotron optical depth: [0.07761866 1.6869038 0.73837553 1.46927273] \n", + "\n", + "Quantum synchrotron dNdt: [3.67619254e+18 6.09509990e+18 3.81990033e+18 3.58546285e+18] \n", + "\n", + "Quantum synchrotron opt depth before evolution: [0.07761866 1.6869038 0.73837553 1.46927273] \n", + "\n", + "Quantum synchrotron opt depth after evolution: [-36.68430676 -59.26409522 -37.46062774 -34.38535576] \n", + "\n", + "Quantum synchrotron part px (before): [ 2.13806811e-21 1.83617015e-21 2.12451506e-21 -2.09011890e-21]\n", + "Quantum synchrotron part py (before): [-2.36370983e-21 8.10882961e-22 -2.39056996e-21 1.46504596e-21]\n", + "Quantum synchrotron part pz (before): [-9.32769361e-22 -1.34089696e-21 2.30881741e-21 9.17695841e-22]\n", + "Quantum synchrotron part px: [ 6.76195755e-22 1.82915759e-21 1.97047892e-21 -2.07583108e-21]\n", + "Quantum synchrotron part py: [-7.47558297e-22 8.07786103e-22 -2.21724374e-21 1.45503107e-21]\n", + "Quantum synchrotron part pz: [-2.95002147e-22 -1.33577592e-21 2.14141859e-21 9.11422579e-22]\n", + "Quantum synchrotron photon px: [ 1.46187235e-21 7.01255172e-24 1.54036140e-22 -1.42878106e-23]\n", + "Quantum synchrotron photon py: [-1.61615153e-21 3.09685827e-24 -1.73326222e-22 1.00148844e-23]\n", + "Quantum synchrotron photon pz: [-6.37767214e-22 -5.12104464e-24 1.67398824e-22 6.27326247e-24]\n" + ] + } + ], + "source": [ + "# Quantum Synchrotron radiation\n", + "\n", + "qs_dndt_params = pxr_qed.qs.dndt_lookup_table_params()\n", + "qs_dndt_params.chi_part_min = 0.01\n", + "qs_dndt_params.chi_part_max = 100\n", + "qs_dndt_params.chi_part_how_many = 128\n", + "print(qs_dndt_params)\n", + "\n", + "qs_dndt_lookup_table = pxr_qed.qs.dndt_lookup_table(qs_dndt_params)\n", + "print(qs_dndt_lookup_table)\n", + "qs_dndt_lookup_table.generate()\n", + "print(qs_dndt_lookup_table)\n", + "\n", + "qs_dndt_lookup_table.save_as(\"qs_dndt.bin\")\n", + "qs_dndt_lookup_table_2 = pxr_qed.qs.dndt_lookup_table()\n", + "print(qs_dndt_lookup_table_2)\n", + "qs_dndt_lookup_table_2.load_from(\"qs_dndt.bin\");\n", + "print(qs_dndt_lookup_table_2)\n", + "\n", + "qs_photem_params = pxr_qed.qs.photon_emission_lookup_table_params()\n", + "qs_photem_params.chi_part_min = 0.01\n", + "qs_photem_params.chi_part_max = 100\n", + "qs_photem_params.frac_min = 1e-12\n", + "qs_photem_params.chi_part_how_many = 64\n", + "qs_photem_params.frac_how_many = 64\n", + "print(qs_photem_params)\n", + "\n", + "qs_dndt_lookup_table.save_as(\"qs_dndt.bin\")\n", + "qs_dndt_lookup_table_2 = pxr_qed.qs.dndt_lookup_table()\n", + "print(qs_dndt_lookup_table_2)\n", + "qs_dndt_lookup_table_2.load_from(\"qs_dndt.bin\");\n", + "print(qs_dndt_lookup_table_2)\n", + "\n", + "qs_photon_emission_lookup_table = pxr_qed.qs.photon_emission_lookup_table(qs_photem_params)\n", + "print(qs_photon_emission_lookup_table)\n", + "qs_photon_emission_lookup_table.generate()\n", + "print(qs_photon_emission_lookup_table)\n", + "\n", + "qs_how_many = 4\n", + "qs_gamma = 10\n", + "qs_dt = 1e-2*fs\n", + "\n", + "qs_Ex = (np.random.rand(qs_how_many)-0.5)*2*E_s\n", + "qs_Ey = (np.random.rand(qs_how_many)-0.5)*2*E_s\n", + "qs_Ez = (np.random.rand(qs_how_many)-0.5)*2*E_s\n", + "qs_Bx = (np.random.rand(qs_how_many)-0.5)*2*B_s\n", + "qs_By = (np.random.rand(qs_how_many)-0.5)*2*B_s\n", + "qs_Bz = (np.random.rand(qs_how_many)-0.5)*2*B_s\n", + "qs_px = (np.random.rand(qs_how_many)-0.5)*2*qs_gamma*m_e*c;\n", + "qs_py = (np.random.rand(qs_how_many)-0.5)*2*qs_gamma*m_e*c;\n", + "qs_pz = (np.random.rand(qs_how_many)-0.5)*2*qs_gamma*m_e*c;\n", + "qs_rand = np.random.rand(qs_how_many)\n", + "qs_ee = np.sqrt(1 + (qs_px**2 + qs_py**2 + qs_pz**2)/((m_e*c)**2))*m_e*c**2;\n", + "\n", + "\n", + "qs_chi = pxr_qed.chi_ele_pos(qs_px, qs_py, qs_pz, qs_Ex, qs_Ey, qs_Ez, qs_Bx, qs_By, qs_Bz)\n", + "print(\"Chi parameters: \", qs_chi, \"\\n\")\n", + "\n", + "qs_interp_dndt = qs_dndt_lookup_table.interp(qs_chi)\n", + "print(\"dN/dt lookup table interpolation: \", qs_interp_dndt, \"\\n\")\n", + "\n", + "qs_interp_photon_emission = qs_photon_emission_lookup_table.interp(qs_chi, qs_rand)\n", + "print(\"Pair production lookup table interpolation: \", qs_interp_photon_emission, \"\\n\")\n", + "\n", + "qs_opt = pxr_qed.qs.get_optical_depth(qs_rand)\n", + "print(\"Quantum synchrotron optical depth: \", qs_opt, \"\\n\")\n", + "\n", + "qs_dndt = pxr_qed.qs.get_dn_dt(qs_ee, qs_chi, qs_dndt_lookup_table)\n", + "print(\"Quantum synchrotron dNdt: \", qs_dndt, \"\\n\")\n", + "\n", + "print(\"Quantum synchrotron opt depth before evolution: \", qs_opt, \"\\n\")\n", + "pxr_qed.qs.evolve_optical_depth(qs_ee, qs_chi, qs_dt, qs_opt, qs_dndt_lookup_table)\n", + "print(\"Quantum synchrotron opt depth after evolution: \", qs_opt, \"\\n\")\n", + "\n", + "print(\"Quantum synchrotron part px (before): \", qs_px)\n", + "print(\"Quantum synchrotron part py (before): \", qs_py)\n", + "print(\"Quantum synchrotron part pz (before): \", qs_pz)\n", + "\n", + "qs_phot_px, qs_phot_py, qs_phot_pz = pxr_qed.qs.generate_photon_update_momentum(\n", + " qs_chi, qs_px, qs_py, qs_pz, qs_rand, qs_photon_emission_lookup_table)\n", + "\n", + "print(\"Quantum synchrotron part px: \", qs_px)\n", + "print(\"Quantum synchrotron part py: \", qs_py)\n", + "print(\"Quantum synchrotron part pz: \", qs_pz)\n", + "print(\"Quantum synchrotron photon px: \", qs_phot_px)\n", + "print(\"Quantum synchrotron photon py: \", qs_phot_py)\n", + "print(\"Quantum synchrotron photon pz: \", qs_phot_pz)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "gross-excerpt", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "pair_production_rate: [4.89384424e+13 3.59995908e+16 1.23279622e-02 1.10208859e+17\n", + " 1.48899978e+16 4.91680000e-15 4.56176309e+15 2.66908513e+17\n", + " 9.06950963e+14 1.96993009e+16]\n", + "expected pair number: [6.36928574e+18 4.68530810e+21 1.60447105e+03 1.43435647e+22\n", + " 1.93791723e+21 6.39916242e-10 5.93708568e+20 3.47378565e+22\n", + " 1.18038694e+20 2.56384286e+21]\n" + ] + } + ], + "source": [ + "#Schwinger pair production\n", + "\n", + "sc_how_many = 10\n", + "\n", + "sc_Ex = (np.random.rand(sc_how_many)-0.5)*2*E_s\n", + "sc_Ey = (np.random.rand(sc_how_many)-0.5)*2*E_s\n", + "sc_Ez = (np.random.rand(sc_how_many)-0.5)*2*E_s\n", + "sc_Bx = (np.random.rand(sc_how_many)-0.5)*2*B_s\n", + "sc_By = (np.random.rand(sc_how_many)-0.5)*2*B_s\n", + "sc_Bz = (np.random.rand(sc_how_many)-0.5)*2*B_s\n", + "\n", + "sc_volume = um**3\n", + "sc_dt = fs\n", + "\n", + "pair_production_rate = pxr_qed.sc.pair_production_rate(sc_Ex, sc_Ey, sc_Ez, sc_Bx, sc_By, sc_Bz, 1)\n", + "exp_pair_number = pxr_qed.sc.expected_pair_number(sc_Ex, sc_Ey, sc_Ez, sc_Bx, sc_By, sc_Bz, sc_volume, sc_dt, 1)\n", + "\n", + "print(\"pair_production_rate:\", pair_production_rate)\n", + "\n", + "print(\"expected pair number:\", exp_pair_number)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/multi_physics/QED/python_bindings/pxr_qed.cpp b/multi_physics/QED/python_bindings/pxr_qed.cpp new file mode 100644 index 000000000..fdc07a481 --- /dev/null +++ b/multi_physics/QED/python_bindings/pxr_qed.cpp @@ -0,0 +1,1358 @@ +/** +* This file contains the python bindings for the PICSAR QED library +*/ + +#include +#include +#include + +#ifdef PXQEDPY_HAS_OPENMP + #include +#endif + +#include "picsar_qed/math/vec_functions.hpp" + +#include "picsar_qed/physics/chi_functions.hpp" + +#include "picsar_qed/physics/breit_wheeler/breit_wheeler_engine_core.hpp" +#include "picsar_qed/physics/breit_wheeler/breit_wheeler_engine_tables_generator.hpp" + +#include "picsar_qed/physics/quantum_sync/quantum_sync_engine_core.hpp" +#include "picsar_qed/physics/quantum_sync/quantum_sync_engine_tables_generator.hpp" + +#include "picsar_qed/physics/schwinger/schwinger_pair_engine_core.hpp" + +#include +#include +#include +#include +#include +#include + +// Some useful aliases +namespace pxr_phys = picsar::multi_physics::phys; +namespace pxr_math = picsar::multi_physics::math; +namespace pxr_bw = picsar::multi_physics::phys::breit_wheeler; +namespace pxr_qs = picsar::multi_physics::phys::quantum_sync; +namespace pxr_sc = picsar::multi_physics::phys::schwinger; +//___________________________________________________________________________________ + + +// Double or single precision +#ifdef PXRQEDPY_DOUBLE_PRECISION + using REAL = double; + const auto PXRQEDPY_PRECISION_STRING = std::string{"double"}; +#else + using REAL = float; + const auto PXRQEDPY_PRECISION_STRING = std::string{"single"}; +#endif +//___________________________________________________________________________________ + + +// Choice of physical units +#if defined(PXRQEDPY_SI) + const auto UU = pxr_phys::unit_system::SI; + const auto PXRQEDPY_USTRING = std::string{"SI"}; +#elif defined(PXRQEDPY_NORM_OMEGA) + const auto UU = pxr_phys::unit_system::norm_omega; + const auto PXRQEDPY_USTRING = std::string{"NORM_OMEGA"}; +#elif defined(PXRQEDPY_NORM_LAMBDA) + const auto UU = pxr_phys::unit_system::norm_lambda; + const auto PXRQEDPY_USTRING = std::string{"NORM_LAMBDA"}; +#elif defined(PXRQEDPY_UNITS) + const auto UU = pxr_phys::unit_system::heaviside_lorentz; + const auto PXRQEDPY_USTRING = std::string{"HEAVISIDE_LORENTZ"}; +#else + #error Incorrect units choice! +#endif +//___________________________________________________________________________________ + + +// PXRQEDPY_FOR can replace a for cycle and uses +// openMP if the library is compiled with openMP support. +#ifdef PXQEDPY_HAS_OPENMP + template + void PXRQEDPY_FOR(int N, const Func& func){ + #pragma omp parallel for + for (int i = 0; i < N; ++i) func(i); + } + const auto PXRQEDPY_OPENMP_FLAG = true; +#else + template + void PXRQEDPY_FOR(int N, const Func& func){ + for (int i = 0; i < N; ++i) func(i); + } + const auto PXRQEDPY_OPENMP_FLAG = true; +#endif +//___________________________________________________________________________________ + + +// Few helper functions + +/** +* Converts a real number into a string using stringstream +* (std::to_string has only 6-digits precision) +* +* @tparam RealType the floating point type +* @param[in] num the number +* @return num as a string +*/ +template +std::string float_to_string(Real num) +{ + std::stringstream ss; + ss << num; + return ss.str(); +} + +/** +* Converts a bool into a string +* +* @param[in] val the bool +* @return either True or False +*/ +std::string bool_to_string(bool val) +{ + return (val)?"True":"False"; +} + +/** +* Throws an error message +* +* @param[in] err_msg the error message +*/ +void throw_error(const std::string& err_msg) +{ + throw std::runtime_error(" [ Error! ] " + err_msg); +} +//___________________________________________________________________________________ + + +// Other useful aliases +namespace py = pybind11; +using pyArr = py::array_t; +using pyBufInfo = py::buffer_info; +using stdVec = std::vector; +using rawVec = std::vector; +//___________________________________________________________________________________ + + +// Helper functions to get raw data pointers from pyArr objects + +/** +* Checks that 'last' is one-dimensional with length len and returns +* a tuple containing a pointer to 'last' raw data +* (uses recursive template metaprogramming, see below) +* +* @param[in] len the length of the array +* @param[in] last a py::array_t object +* @return a tuple containing a pointer to the raw data of last +*/ +auto aux_check_and_get_pointers(const long int len, const pyArr& last) +{ + const auto last_buf = last.request(); + if (last_buf.ndim != 1 || last_buf.shape[0] != len) + throw_error("All arrays must be one-dimensional with equal size"); + + const auto cptr = static_cast(last_buf.ptr); + return std::make_tuple(cptr); +} + +/** +* Checks that 'arg' and 'args' are one-dimensional with length len +* and puts pointers to their raw data into a tuple +* (uses recursive template metaprogramming, see above and below) +* +* @tparam ...Args a variable number of const py::array_t & +* @param[in] len the length of the arrays +* @param[in] arg the first py::array_t in the argument list +* @param[in] args the other arguments +* @return a tuple containing pointers to 'arg' and 'args' raw data +*/ +template +auto aux_check_and_get_pointers(const long int len, const pyArr& arg, const Args& ...args) +{ + const auto arg_buf = arg.request(); + if (arg_buf.ndim != 1 || arg_buf.shape[0] != len) + throw_error("All arrays must be one-dimensional with equal size"); + + const auto cptr = static_cast(arg_buf.ptr); + return std::tuple_cat(std::make_tuple(cptr), + aux_check_and_get_pointers(len, args...)); +} + +/** +* Checks that 'first' and 'args' are one-dimensional with the same length +* and puts pointers to their raw data into a tuple +* (uses recursive template metaprogramming, see above) +* +* @tparam ...Args a variable number of const py::array_t & +* @param[in] len the length of the arrays +* @param[in] first the first py::array_t in the argument list +* @param[in] args the other arguments +* @return a tuple containing the length of the arrays and pointers to 'first' and 'args' raw data +*/ +template +auto check_and_get_pointers(const pyArr& first, const Args& ...args) +{ + const auto first_buf = first.request(); + if (first_buf.ndim != 1) + throw_error("All arrays must be one-dimensional with equal size"); + + const auto len = first_buf.shape[0]; + + const auto cptr = static_cast(first_buf.ptr); + + return std::tuple_cat(std::make_tuple(len, cptr), + aux_check_and_get_pointers(len, args...)); +} + +/** +* Overload of 'check_and_get_pointers' for one array: checks that +* the array is one-dimensional and returns a tuple +* +* @param[in] arr a py::array_t +* @return a tuple containing the length of the array and a pointer to 'arr' raw data +*/ +auto check_and_get_pointers(const pyArr& arr) +{ + const auto arr_buf = arr.request(); + if (arr_buf.ndim != 1) + throw_error("Array must be one-dimensional"); + + const auto len = arr_buf.shape[0]; + + const auto cptr = static_cast(arr_buf.ptr); + + return std::make_tuple(len, cptr); +} + +/** +* Checks that arr has a given length and returns a pointer to +* its raw data. +* +* @param[in] arr a py::array_t +* @param[in] len the array length +* @return a pointer to 'arr' raw data +*/ +auto check_and_get_pointer_nonconst(pyArr& arr, const long int len) +{ + const auto arr_buf = arr.request(); + if (arr_buf.ndim != 1 || arr_buf.shape[0] != len) + throw_error("Array must be one-dimensional with size " + std::to_string(len)); + + const auto cptr = static_cast(arr_buf.ptr); + + return cptr; +} +//___________________________________________________________________________________ + + +// ******************************* Chi parameters wrappers ************************** + +/** +* Wrapper for chi_photon function +* +* @param[in] px x components of photon momenta +* @param[in] py y components of photon momenta +* @param[in] pz z components of photon momenta +* @param[in] ex x components of electric field +* @param[in] ey y components of electric field +* @param[in] ez z components of electric field +* @param[in] bx x components of magnetic field +* @param[in] by y components of magnetic field +* @param[in] bz z components of magnetic field +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return the chi parameter of the photons +*/ +pyArr +chi_photon_wrapper( + const pyArr& px, const pyArr& py, const pyArr& pz, + const pyArr& ex, const pyArr& ey, const pyArr& ez, + const pyArr& bx, const pyArr& by, const pyArr& bz, + const REAL ref_quantity) +{ + const REAL + *p_px = nullptr, *p_py = nullptr, *p_pz = nullptr, + *p_ex = nullptr, *p_ey = nullptr, *p_ez = nullptr, + *p_bx = nullptr, *p_by = nullptr, *p_bz = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_px, p_py, p_pz, + p_ex, p_ey, p_ez, + p_bx, p_by, p_bz) = + check_and_get_pointers( + px, py, pz, + ex, ey, ez, + bx, by, bz); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_phys::chi_photon( + p_px[i], p_py[i], p_pz[i], + p_ex[i], p_ey[i], p_ez[i], + p_bx[i], p_by[i], p_bz[i], + ref_quantity); + }); + + return res; +} + +/** +* Wrapper for chi_ele_pos function +* +* @param[in] px x components of particle momenta +* @param[in] py y components of particle momenta +* @param[in] pz z components of particle momenta +* @param[in] ex x components of electric field +* @param[in] ey y components of electric field +* @param[in] ez z components of electric field +* @param[in] bx x components of magnetic field +* @param[in] by y components of magnetic field +* @param[in] bz z components of magnetic field +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return the chi parameter of the particles +*/ +pyArr +chi_ele_pos_wrapper( + const pyArr& px, const pyArr& py, const pyArr& pz, + const pyArr& ex, const pyArr& ey, const pyArr& ez, + const pyArr& bx, const pyArr& by, const pyArr& bz, + const REAL ref_quantity) +{ + const REAL + *p_px = nullptr, *p_py = nullptr, *p_pz = nullptr, + *p_ex = nullptr, *p_ey = nullptr, *p_ez = nullptr, + *p_bx = nullptr, *p_by = nullptr, *p_bz = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_px, p_py, p_pz, + p_ex, p_ey, p_ez, + p_bx, p_by, p_bz) = + check_and_get_pointers( + px, py, pz, + ex, ey, ez, + bx, by, bz); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_phys::chi_ele_pos( + p_px[i], p_py[i], p_pz[i], + p_ex[i], p_ey[i], p_ez[i], + p_bx[i], p_by[i], p_bz[i], + ref_quantity); + }); + + return res; +} +// ______________________________________________________________________________________________ + + +// ******************************* Breit-Wheeler pair production wrappers *********************** + +// Aliases for table objects and parameter structs +using bw_dndt_lookup_table_params = + pxr_bw::dndt_lookup_table_params; + +using bw_pair_prod_lookup_table_params = + pxr_bw::pair_prod_lookup_table_params; + +using bw_dndt_lookup_table = + pxr_bw::dndt_lookup_table; + +using bw_pair_prod_lookup_table = + pxr_bw::pair_prod_lookup_table; + +// Aliases for table generation policies +const auto bw_regular = + pxr_bw::generation_policy::regular; + +const auto bw_force_double = + pxr_bw::generation_policy::force_internal_double; + +/** +* Wrapper for Breit-Wheeler get_optical_depth function +* +* @param[in] unf_zero_one_minus_epsi an array of random numbers uniformly distributed in [0,1) +* @return the optical depths drawn from an exponential distribution +*/ +pyArr +bw_get_optical_depth_wrapper( + const pyArr& unf_zero_one_minus_epsi) +{ + const REAL + *p_unf_zero_one_minus_epsi = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, p_unf_zero_one_minus_epsi)= + check_and_get_pointers(unf_zero_one_minus_epsi); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_bw::get_optical_depth( + p_unf_zero_one_minus_epsi[i]); + }); + + return res; +} + +/** +* Wrapper for Breit-Wheeler get_dn_dt function +* +* @param[in] energy_phot the energy of the photons +* @param[in] chi_phot the chi parameters of the photons +* @param[in] ref_table the dn_dt lookup table +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return the Breit-Wheeler pair production rates +*/ +pyArr +bw_get_dn_dt_wrapper( + const pyArr& energy_phot, const pyArr& chi_phot, + const bw_dndt_lookup_table& ref_table, + const REAL ref_quantity) +{ + const REAL + *p_energy_phot = nullptr, *p_chi_phot = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_energy_phot, p_chi_phot) = + check_and_get_pointers( + energy_phot, chi_phot); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_bw::get_dN_dt( + p_energy_phot[i], p_chi_phot[i], + ref_table, ref_quantity); + }); + + return res; +} + +/** +* Wrapper for Breit-Wheeler evolve_optical_depth function +* +* @param[in] energy_phot the energy of the photons +* @param[in] chi_phot the chi parameters of the photons +* @param[in] dt the timestep +* @param[in, out] optical_depth the optical depth of the photons +* @param[in] ref_table the dn_dt lookup table +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +*/ +void +bw_evolve_optical_depth_wrapper( + const pyArr& energy_phot, const pyArr& chi_phot, + const REAL dt, pyArr& optical_depth, + const bw_dndt_lookup_table& ref_table, + const REAL ref_quantity) +{ + const REAL + *p_energy_phot = nullptr, *p_chi_phot = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_energy_phot, p_chi_phot) = + check_and_get_pointers( + energy_phot, chi_phot); + + auto p_optical_depth = + check_and_get_pointer_nonconst(optical_depth, how_many); + + PXRQEDPY_FOR(how_many, [&](int i){ + pxr_bw::evolve_optical_depth( + p_energy_phot[i], p_chi_phot[i], + dt, p_optical_depth[i], + ref_table, ref_quantity); + }); +} + +/** +* Wrapper for Breit-Wheeler generate_breit_wheeler_pairs function +* +* @param[in] chi_phot the chi parameters of the photons +* @param[in] phot_px x components of photon momenta +* @param[in] phot_py y components of photon momenta +* @param[in] phot_pz z components of photon momenta +* @param[in] unf_zero_one_minus_epsi an array of random numbers uniformly distributed in [0,1) +* @param[in] ref_table the pair production lookup table +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return a tuple containing the components of the momenta of the generated particles (ele_px, ele_py,..,pos_pz) +*/ +auto +bw_generate_breit_wheeler_pairs_wrapper( + const pyArr& chi_phot, + const pyArr& phot_px, const pyArr& phot_py, const pyArr& phot_pz, + const pyArr& unf_zero_one_minus_epsi, + const bw_pair_prod_lookup_table& ref_table, + const REAL ref_quantity) +{ + const REAL + *p_chi_phot = nullptr, + *p_phot_px = nullptr, *p_phot_py = nullptr, *p_phot_pz = nullptr, + *p_unf_zero_one_minus_epsi; + + size_t how_many = 0; + + std::tie( + how_many, + p_chi_phot, p_phot_px, p_phot_py, p_phot_pz, + p_unf_zero_one_minus_epsi) = + check_and_get_pointers( + chi_phot, phot_px, phot_py, phot_pz, unf_zero_one_minus_epsi); + + auto ele_px = pyArr(how_many); + auto ele_py = pyArr(how_many); + auto ele_pz = pyArr(how_many); + auto pos_px = pyArr(how_many); + auto pos_py = pyArr(how_many); + auto pos_pz = pyArr(how_many); + auto p_ele_px = static_cast(ele_px.request().ptr); + auto p_ele_py = static_cast(ele_py.request().ptr); + auto p_ele_pz = static_cast(ele_pz.request().ptr); + auto p_pos_px = static_cast(pos_px.request().ptr); + auto p_pos_py = static_cast(pos_py.request().ptr); + auto p_pos_pz = static_cast(pos_pz.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + auto ele_mom = pxr_math::vec3{}; + auto pos_mom = pxr_math::vec3{}; + pxr_bw::generate_breit_wheeler_pairs( + p_chi_phot[i], + pxr_math::vec3{p_phot_px[i], p_phot_py[i], p_phot_pz[i]}, + p_unf_zero_one_minus_epsi[i], + ref_table, + ele_mom, pos_mom, + ref_quantity); + + p_ele_px[i] = ele_mom[0]; + p_ele_py[i] = ele_mom[1]; + p_ele_pz[i] = ele_mom[2]; + + p_pos_px[i] = pos_mom[0]; + p_pos_py[i] = pos_mom[1]; + p_pos_pz[i] = pos_mom[2]; + }); + + return std::make_tuple( + std::move(ele_px),std::move(ele_py), std::move(ele_pz), + std::move(pos_px),std::move(pos_py), std::move(pos_pz)); +} +// ______________________________________________________________________________________________ + + +// ******************************* Quantum synchrotron emission wrappers ************************ + +// Aliases for table objects and parameter structs +using qs_dndt_lookup_table_params = + pxr_qs::dndt_lookup_table_params; + +using qs_photon_emission_lookup_table_params = + pxr_qs::photon_emission_lookup_table_params; + +using qs_dndt_lookup_table = + pxr_qs::dndt_lookup_table; + +using qs_photon_emission_lookup_table = + pxr_qs::photon_emission_lookup_table; + +// Aliases for table generation policies +const auto qs_regular = + pxr_qs::generation_policy::regular; + +const auto qs_force_double = + pxr_qs::generation_policy::force_internal_double; + +/** +* Wrapper for Quantum Synchrotron get_optical_depth function +* +* @param[in] unf_zero_one_minus_epsi an array of random numbers uniformly distributed in [0,1) +* @return the optical depths drawn from an exponential distribution +*/ +pyArr +qs_get_optical_depth_wrapper( + const pyArr& unf_zero_one_minus_epsi) +{ + const REAL + *p_unf_zero_one_minus_epsi = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, p_unf_zero_one_minus_epsi)= + check_and_get_pointers(unf_zero_one_minus_epsi); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_qs::get_optical_depth( + p_unf_zero_one_minus_epsi[i]); + }); + + return res; +} + +/** +* Wrapper for Quantum Synchrotron get_dn_dt function +* +* @param[in] energy_part the energy of the particles +* @param[in] chi_part the chi parameters of the particles +* @param[in] ref_table the dn_dt lookup table +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return the Quantum Synchrotron photon emission +*/ +pyArr +qs_get_dn_dt_wrapper( + const pyArr& energy_part, const pyArr& chi_part, + const qs_dndt_lookup_table& ref_table, + const REAL ref_quantity) +{ + const REAL + *p_energy_part = nullptr, *p_chi_part = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_energy_part, p_chi_part) = + check_and_get_pointers( + energy_part, chi_part); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_qs::get_dN_dt( + p_energy_part[i], p_chi_part[i], + ref_table, ref_quantity); + }); + + return res; +} + +/** +* Wrapper for Quantum Synchrotron evolve_optical_depth function +* +* @param[in] energy_part the energy of the particles +* @param[in] chi_phot the chi parameters of the particles +* @param[in] dt the timestep +* @param[in, out] optical_depth the optical depth of the particles +* @param[in] ref_table the dn_dt lookup table +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +*/ +void +qs_evolve_optical_depth_wrapper( + const pyArr& energy_part, const pyArr& chi_part, + const REAL dt, pyArr& optical_depth, + const qs_dndt_lookup_table& ref_table, + const REAL ref_quantity) +{ + const REAL + *p_energy_part = nullptr, *p_chi_part = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_energy_part, p_chi_part) = + check_and_get_pointers( + energy_part, chi_part); + + auto p_optical_depth = + check_and_get_pointer_nonconst(optical_depth, how_many); + + PXRQEDPY_FOR(how_many, [&](int i){ + pxr_qs::evolve_optical_depth( + p_energy_part[i], p_chi_part[i], + dt, p_optical_depth[i], + ref_table, ref_quantity); + }); +} + +/** +* Wrapper for Quantum Synchrotron generate_photon_update_momentum function +* +* @param[in] chi_part the chi parameters of the particles +* @param[in,out] part_px x components of particle momenta +* @param[in,out] part_py y components of particle momenta +* @param[in,out] part_pz z components of particle momenta +* @param[in] unf_zero_one_minus_epsi an array of random numbers uniformly distributed in [0,1) +* @param[in] ref_table the photon emission lookup table +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return a tuple containing the components of the momenta of the generated photons +*/ +auto +qs_generate_photon_update_momentum_wrapper( + const pyArr& chi_part, + pyArr& part_px, pyArr& part_py, pyArr& part_pz, + const pyArr& unf_zero_one_minus_epsi, + const qs_photon_emission_lookup_table& ref_table, + const REAL ref_quantity) +{ + const REAL + *p_chi_part = nullptr, *p_unf_zero_one_minus_epsi = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_chi_part, p_unf_zero_one_minus_epsi) = + check_and_get_pointers(chi_part, unf_zero_one_minus_epsi); + + auto p_part_px = + check_and_get_pointer_nonconst(part_px, how_many); + auto p_part_py = + check_and_get_pointer_nonconst(part_py, how_many); + auto p_part_pz = + check_and_get_pointer_nonconst(part_pz, how_many); + + auto phot_px = pyArr(how_many); + auto phot_py = pyArr(how_many); + auto phot_pz = pyArr(how_many); + auto p_phot_px = static_cast(phot_px.request().ptr); + auto p_phot_py = static_cast(phot_py.request().ptr); + auto p_phot_pz = static_cast(phot_pz.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + auto part_mom = pxr_math::vec3{p_part_px[i], p_part_py[i], p_part_pz[i]}; + auto phot_mom = pxr_math::vec3{}; + pxr_qs::generate_photon_update_momentum( + p_chi_part[i], + part_mom, + p_unf_zero_one_minus_epsi[i], + ref_table, + phot_mom, + ref_quantity); + + p_part_px[i] = part_mom[0]; + p_part_py[i] = part_mom[1]; + p_part_pz[i] = part_mom[2]; + p_phot_px[i] = phot_mom[0]; + p_phot_py[i] = phot_mom[1]; + p_phot_pz[i] = phot_mom[2]; + }); + + return std::make_tuple( + std::move(phot_px),std::move(phot_py), std::move(phot_pz)); +} + +// ______________________________________________________________________________________________ + + +// ******************************* Schwinger pair production ************************************* + +/** +* Wrapper for Schwinger pair_production_rate function +* +* @param[in] ex x components of electric field +* @param[in] ey y components of electric field +* @param[in] ez z components of electric field +* @param[in] bx x components of magnetic field +* @param[in] by y components of magnetic field +* @param[in] bz z components of magnetic field +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return the Schwinger pair production rate per unit volume +*/ +pyArr +sc_pair_production_rate_wrapper( + const pyArr& ex, const pyArr& ey, const pyArr& ez, + const pyArr& bx, const pyArr& by, const pyArr& bz, + const REAL ref_quantity) +{ + const REAL + *p_ex = nullptr, *p_ey = nullptr, *p_ez = nullptr, + *p_bx = nullptr, *p_by = nullptr, *p_bz = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_ex, p_ey, p_ez, + p_bx, p_by, p_bz) = + check_and_get_pointers( + ex, ey, ez, + bx, by, bz); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_sc::pair_production_rate( + p_ex[i], p_ey[i], p_ez[i], + p_bx[i], p_by[i], p_bz[i], + ref_quantity); + }); + + return res; +} + +/** +* Wrapper for Schwinger expected_pair_number function +* +* @param[in] ex x components of electric field +* @param[in] ey y components of electric field +* @param[in] ez z components of electric field +* @param[in] bx x components of magnetic field +* @param[in] by y components of magnetic field +* @param[in] bz z components of magnetic field +* @param[in] volume the volume of the region where pair production takes place +* @param[in] dt the timestep +* @param[in] ref_quantity reference quantity (for NORM_LAMBDA or NORM_OMEGA units) +* @return the expected number of Schwinger pairs +*/ +pyArr +sc_expected_pair_number_wrapper( + const pyArr& ex, const pyArr& ey, const pyArr& ez, + const pyArr& bx, const pyArr& by, const pyArr& bz, + const REAL volume, const REAL dt, + const REAL ref_quantity) +{ + const REAL + *p_ex = nullptr, *p_ey = nullptr, *p_ez = nullptr, + *p_bx = nullptr, *p_by = nullptr, *p_bz = nullptr; + + size_t how_many = 0; + + std::tie( + how_many, + p_ex, p_ey, p_ez, + p_bx, p_by, p_bz) = + check_and_get_pointers( + ex, ey, ez, + bx, by, bz); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = + pxr_sc::expected_pair_number( + p_ex[i], p_ey[i], p_ez[i], + p_bx[i], p_by[i], p_bz[i], + volume, dt, + ref_quantity); + }); + + return res; +} + +// ______________________________________________________________________________________________ + + +// ******************************* Python module ************************************************* + +PYBIND11_MODULE(pxr_qed, m) { + // A doc string and few useful variables + m.doc() = "pybind11 pxr_qed plugin"; + + m.attr("PRECISION") = py::str(PXRQEDPY_PRECISION_STRING); + m.attr("HAS_OPENMP") = py::bool_(PXRQEDPY_OPENMP_FLAG); + m.attr("UNITS") = py::str(PXRQEDPY_USTRING); + //________________________________________ + + + // Functions to calculate chi parameters + m.def( + "chi_photon", + &chi_photon_wrapper, + py::arg("px").noconvert(true), py::arg("py").noconvert(true), py::arg("pz").noconvert(true), + py::arg("ex").noconvert(true), py::arg("ey").noconvert(true), py::arg("ez").noconvert(true), + py::arg("bx").noconvert(true), py::arg("by").noconvert(true), py::arg("bz").noconvert(true), + py::arg("ref_quantity") = py::float_(1.0) + ); + + m.def( + "chi_ele_pos", + &chi_ele_pos_wrapper, + py::arg("px").noconvert(true), py::arg("py").noconvert(true), py::arg("pz").noconvert(true), + py::arg("ex").noconvert(true), py::arg("ey").noconvert(true), py::arg("ez").noconvert(true), + py::arg("bx").noconvert(true), py::arg("by").noconvert(true), py::arg("bz").noconvert(true), + py::arg("ref_quantity") = py::float_(1.0) + ); + //________________________________________ + + // Breit-Wheeler submodule + auto bw = m.def_submodule( "bw" ); + + bw.def( + "get_optical_depth", + &bw_get_optical_depth_wrapper, + "Computes the optical depth of a new photon", + py::arg("unf_zero_one_minus_epsi").noconvert(true)); + + py::class_(bw, + "dndt_lookup_table_params", + "Parameters to generate a dN/dt lookup table") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &bw_dndt_lookup_table_params::operator==) + .def_readwrite("chi_phot_min", &bw_dndt_lookup_table_params::chi_phot_min) + .def_readwrite("chi_phot_max", &bw_dndt_lookup_table_params::chi_phot_max) + .def_readwrite("chi_phot_how_many", &bw_dndt_lookup_table_params::chi_phot_how_many) + .def("__repr__", + [](const bw_dndt_lookup_table_params &a) { + return + std::string("bw.dndt_lookup_table_params:\n")+ + std::string("\tchi_phot_min : ") + float_to_string(a.chi_phot_min)+"\n"+ + std::string("\tchi_phot_max : ") + float_to_string(a.chi_phot_max)+"\n"+ + std::string("\tchi_phot_how_many: ") + std::to_string(a.chi_phot_how_many); + }); + + py::class_(bw, + "pair_prod_lookup_table_params") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &bw_pair_prod_lookup_table_params::operator==) + .def_readwrite("chi_phot_min", &bw_pair_prod_lookup_table_params::chi_phot_min) + .def_readwrite("chi_phot_max", &bw_pair_prod_lookup_table_params::chi_phot_max) + .def_readwrite("chi_phot_how_many", &bw_pair_prod_lookup_table_params::chi_phot_how_many) + .def_readwrite("frac_how_many", &bw_pair_prod_lookup_table_params::frac_how_many) + .def("__repr__", + [](const bw_pair_prod_lookup_table_params &a) { + return + std::string("bw.pair_prod_lookup_table_params:\n")+ + std::string("\tchi_phot_min : ") + float_to_string(a.chi_phot_min)+"\n"+ + std::string("\tchi_phot_max : ") + float_to_string(a.chi_phot_max)+"\n"+ + std::string("\tchi_phot_how_many: ") + std::to_string(a.chi_phot_how_many)+"\n"+ + std::string("\tfrac_how_many : ") + std::to_string(a.frac_how_many); + + }); + + py::class_(bw, + "dndt_lookup_table", + "dN/dt lookup table") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &bw_dndt_lookup_table::operator==) + .def("generate", + [&](bw_dndt_lookup_table &self, + bool do_regular, bool verbose){ + if(do_regular) + self.generate(verbose); + else + self.generate(verbose); + }, + py::arg("do_regular") = py::bool_(true), + py::arg("verbose") = py::bool_(true)) + .def("save_as", + [&](const bw_dndt_lookup_table &self, const std::string file_name){ + if(!self.is_init()) + throw_error("Table must be initialized!"); + const auto raw = self.serialize(); + auto of = std::fstream(file_name, + std::ios::out | std::ios::binary); + if( !of ) + throw_error("Opening file failed!"); + of.write(raw.data(), raw.size()); + of.close(); + }, + py::arg("file_name")) + .def("load_from", + [&](bw_dndt_lookup_table &self, const std::string file_name){ + auto input = std::ifstream(file_name, + std::ios::ate | std::ios::binary); + if( !input ) + throw_error("Opening file failed!"); + const auto pos = input.tellg(); + auto raw = rawVec(pos); + + input.seekg(0, std::ios::beg); + input.read(raw.data(), pos); + + self = bw_dndt_lookup_table{raw}; + input.close(); + }, + py::arg("file_name")) + .def("interp", + [&](bw_dndt_lookup_table &self, const pyArr& chi_phot){ + const REAL* p_chi_phot = nullptr; + size_t how_many = 0; + std::tie(how_many, p_chi_phot)= + check_and_get_pointers(chi_phot); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = self.interp(p_chi_phot[i]); + }); + return res; + }, + py::arg("chi_phot")) + .def("__repr__", + [](const bw_dndt_lookup_table &a) { + return + std::string("bw.dndt_lookup_table:\n")+ + std::string("\tis initialized? : ") + bool_to_string(a.is_init())+"\n"; + }); + + py::class_(bw, + "pair_prod_lookup_table", + "Pair production lookup table") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &bw_pair_prod_lookup_table::operator==) + .def("generate", + [&](bw_pair_prod_lookup_table &self, + bool do_regular, bool verbose){ + if(do_regular) + self.generate(verbose); + else + self.generate(verbose); + }, + py::arg("do_regular") = py::bool_(true), + py::arg("verbose") = py::bool_(true)) + .def("save_as", + [&](const bw_pair_prod_lookup_table &self, const std::string file_name){ + if(!self.is_init()) + throw_error("Table must be initialized!"); + const auto raw = self.serialize(); + auto of = std::fstream(file_name, + std::ios::out | std::ios::binary); + if( !of ) + throw_error("Opening file failed!"); + of.write(raw.data(), raw.size()); + of.close(); + }, + py::arg("file_name")) + .def("load_from", + [&](bw_pair_prod_lookup_table &self, const std::string file_name){ + auto input = std::ifstream(file_name, + std::ios::ate | std::ios::binary); + if( !input ) + throw_error("Opening file failed!"); + const auto pos = input.tellg(); + auto raw = rawVec(pos); + + input.seekg(0, std::ios::beg); + input.read(raw.data(), pos); + + self = bw_pair_prod_lookup_table{raw}; + input.close(); + }, + py::arg("file_name")) + .def("interp", + [&](bw_pair_prod_lookup_table &self, + const pyArr& chi_phot, const pyArr& unf_zero_one_minus_epsi){ + const REAL + *p_chi_phot = nullptr, *p_unf_zero_one_minus_epsi = nullptr; + size_t how_many = 0; + std::tie(how_many, p_chi_phot, p_unf_zero_one_minus_epsi)= + check_and_get_pointers(chi_phot, unf_zero_one_minus_epsi); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = self.interp(p_chi_phot[i], p_unf_zero_one_minus_epsi[i]); + }); + return res; + }, + py::arg("chi_phot"), + py::arg("unf_zero_one_minus_epsi")) + .def("__repr__", + [](const bw_pair_prod_lookup_table &a) { + return + std::string("bw.pair_prod_lookup_table:\n")+ + std::string("\tis initialized? : ") + bool_to_string(a.is_init())+"\n"; + }); + + bw.def( + "get_dn_dt", + &bw_get_dn_dt_wrapper, + py::arg("energy_phot").noconvert(true), + py::arg("chi_phot").noconvert(true), + py::arg("ref_table"), py::arg("ref_quantity") = py::float_(1.0) + ); + + bw.def( + "evolve_optical_depth", + &bw_evolve_optical_depth_wrapper, + py::arg("energy_phot").noconvert(true), + py::arg("chi_phot").noconvert(true), + py::arg("dt"), py::arg("optical_depth").noconvert(true), + py::arg("ref_table"), py::arg("ref_quantity") = py::float_(1.0) + ); + + bw.def( + "generate_breit_wheeler_pairs", + &bw_generate_breit_wheeler_pairs_wrapper, + py::arg("chi_phot").noconvert(true), + py::arg("phot_px").noconvert(true), + py::arg("phot_py").noconvert(true), + py::arg("phot_pz").noconvert(true), + py::arg("unf_zero_one_minus_epsi").noconvert(true), + py::arg("ref_table"), py::arg("ref_quantity") = py::float_(1.0) + ); + + //________________________________________ + + + //Quantum Synchrotron submodule + auto qs = m.def_submodule( "qs" ); + + qs.def( + "get_optical_depth", + &bw_get_optical_depth_wrapper, + "Computes the optical depth of a new electron or positron", + py::arg("unf_zero_one_minus_epsi").noconvert(true)); + + py::class_(qs, + "dndt_lookup_table_params", + "Parameters to generate a dN/dt lookup table") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &qs_dndt_lookup_table_params::operator==) + .def_readwrite("chi_part_min", &qs_dndt_lookup_table_params::chi_part_min) + .def_readwrite("chi_part_max", &qs_dndt_lookup_table_params::chi_part_max) + .def_readwrite("chi_part_how_many", &qs_dndt_lookup_table_params::chi_part_how_many) + .def("__repr__", + [](const qs_dndt_lookup_table_params &a) { + return + std::string("qs.dndt_lookup_table_params:\n")+ + std::string("\tchi_part_min : ") + float_to_string(a.chi_part_min)+"\n"+ + std::string("\tchi_part_max : ") + float_to_string(a.chi_part_max)+"\n"+ + std::string("\tchi_part_how_many: ") + std::to_string(a.chi_part_how_many); + }); + + py::class_(qs, + "photon_emission_lookup_table_params") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &qs_photon_emission_lookup_table_params::operator==) + .def_readwrite("chi_part_min", &qs_photon_emission_lookup_table_params::chi_part_min) + .def_readwrite("chi_part_max", &qs_photon_emission_lookup_table_params::chi_part_max) + .def_readwrite("frac_min", &qs_photon_emission_lookup_table_params::frac_min) + .def_readwrite("chi_part_how_many", &qs_photon_emission_lookup_table_params::chi_part_how_many) + .def_readwrite("frac_how_many", &qs_photon_emission_lookup_table_params::frac_how_many) + .def("__repr__", + [](const qs_photon_emission_lookup_table_params &a) { + return + std::string("qs.photon_emission_lookup_table_params:\n")+ + std::string("\tchi_part_min : ") + float_to_string(a.chi_part_min)+"\n"+ + std::string("\tchi_part_max : ") + float_to_string(a.chi_part_max)+"\n"+ + std::string("\tfrac_min : ") + float_to_string(a.frac_min)+"\n"+ + std::string("\tchi_part_how_many: ") + std::to_string(a.chi_part_how_many)+"\n"+ + std::string("\tfrac_how_many : ") + std::to_string(a.frac_how_many); + }); + + py::class_(qs, + "dndt_lookup_table", + "dN/dt lookup table") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &qs_dndt_lookup_table::operator==) + .def("generate", + [&](qs_dndt_lookup_table &self, + bool do_regular, bool verbose){ + if(do_regular) + self.generate(verbose); + else + self.generate(verbose); + }, + py::arg("do_regular") = py::bool_(true), + py::arg("verbose") = py::bool_(true)) + .def("save_as", + [&](const qs_dndt_lookup_table &self, const std::string file_name){ + if(!self.is_init()) + throw_error("Table must be initialized!"); + const auto raw = self.serialize(); + auto of = std::fstream(file_name, + std::ios::out | std::ios::binary); + if( !of ) + throw_error("Opening file failed!"); + of.write(raw.data(), raw.size()); + of.close(); + }, + py::arg("file_name")) + .def("load_from", + [&](qs_dndt_lookup_table &self, const std::string file_name){ + auto input = std::ifstream(file_name, + std::ios::ate | std::ios::binary); + if( !input ) + throw_error("Opening file failed!"); + const auto pos = input.tellg(); + auto raw = rawVec(pos); + + input.seekg(0, std::ios::beg); + input.read(raw.data(), pos); + + self = qs_dndt_lookup_table{raw}; + input.close(); + }, + py::arg("file_name")) + .def("interp", + [&](qs_dndt_lookup_table &self, const pyArr& chi_part){ + const REAL* p_chi_part = nullptr; + size_t how_many = 0; + std::tie(how_many, p_chi_part)= + check_and_get_pointers(chi_part); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = self.interp(p_chi_part[i]); + }); + return res; + }, + py::arg("chi_part")) + .def("__repr__", + [](const qs_dndt_lookup_table &a) { + return + std::string("qs.dndt_lookup_table:\n")+ + std::string("\tis initialized? : ") + bool_to_string(a.is_init())+"\n"; + }); + + py::class_(qs, + "photon_emission_lookup_table", + "Photon emission lookup table") + .def(py::init<>()) + .def(py::init()) + .def("__eq__", &qs_photon_emission_lookup_table::operator==) + .def("generate", + [&](qs_photon_emission_lookup_table &self, + bool do_regular, bool verbose){ + if(do_regular) + self.generate(verbose); + else + self.generate(verbose); + }, + py::arg("do_regular") = py::bool_(true), + py::arg("verbose") = py::bool_(true)) + .def("save_as", + [&](const qs_photon_emission_lookup_table &self, const std::string file_name){ + if(!self.is_init()) + throw_error("Table must be initialized!"); + const auto raw = self.serialize(); + auto of = std::fstream(file_name, + std::ios::out | std::ios::binary); + if( !of ) + throw_error("Opening file failed!"); + of.write(raw.data(), raw.size()); + of.close(); + }, + py::arg("file_name")) + .def("load_from", + [&](qs_photon_emission_lookup_table &self, const std::string file_name){ + auto input = std::ifstream(file_name, + std::ios::ate | std::ios::binary); + if( !input ) + throw_error("Opening file failed!"); + const auto pos = input.tellg(); + auto raw = rawVec(pos); + + input.seekg(0, std::ios::beg); + input.read(raw.data(), pos); + + self = qs_photon_emission_lookup_table{raw}; + input.close(); + }, + py::arg("file_name")) + .def("interp", + [&](qs_photon_emission_lookup_table &self, + const pyArr& chi_part, const pyArr& unf_zero_one_minus_epsi){ + const REAL + *p_chi_part = nullptr, *p_unf_zero_one_minus_epsi = nullptr; + size_t how_many = 0; + std::tie(how_many, p_chi_part, p_unf_zero_one_minus_epsi)= + check_and_get_pointers(chi_part, unf_zero_one_minus_epsi); + + auto res = pyArr(how_many); + auto p_res = static_cast(res.request().ptr); + + PXRQEDPY_FOR(how_many, [&](int i){ + p_res[i] = self.interp(p_chi_part[i], p_unf_zero_one_minus_epsi[i]); + }); + return res; + }, + py::arg("chi_part"), + py::arg("unf_zero_one_minus_epsi")) + .def("__repr__", + [](const qs_photon_emission_lookup_table &a) { + return + std::string("qs.pair_prod_lookup_table:\n")+ + std::string("\tis initialized? : ") + bool_to_string(a.is_init())+"\n"; + }); + + qs.def( + "get_dn_dt", + &qs_get_dn_dt_wrapper, + py::arg("energy_part").noconvert(true), + py::arg("chi_part").noconvert(true), + py::arg("ref_table"), py::arg("ref_quantity") = py::float_(1.0) + ); + + qs.def( + "evolve_optical_depth", + &qs_evolve_optical_depth_wrapper, + py::arg("energy_part").noconvert(true), + py::arg("chi_part").noconvert(true), + py::arg("dt"), py::arg("optical_depth").noconvert(true), + py::arg("ref_table"), py::arg("ref_quantity") = py::float_(1.0) + ); + + qs.def( + "generate_photon_update_momentum", + &qs_generate_photon_update_momentum_wrapper, + py::arg("chi_part").noconvert(true), + py::arg("part_px").noconvert(true), + py::arg("part_py").noconvert(true), + py::arg("part_pz").noconvert(true), + py::arg("unf_zero_one_minus_epsi").noconvert(true), + py::arg("ref_table"), py::arg("ref_quantity") = py::float_(1.0) + ); + //________________________________________ + + //Schwinger submodule + auto sc = m.def_submodule( "sc" ); + + sc.def("pair_production_rate", + &sc_pair_production_rate_wrapper, + "Computes the Schwinger pair production rate using the Nikishov formula", + py::arg("ex").noconvert(true), py::arg("ey").noconvert(true), py::arg("ez").noconvert(true), + py::arg("bx").noconvert(true), py::arg("by").noconvert(true), py::arg("bz").noconvert(true), + py::arg("ref_quantity") = py::float_(1.0) + ); + + sc.def("expected_pair_number", + &sc_expected_pair_number_wrapper, + "Computes the expected number of Schwinger pairs using the Nikishov formula", + py::arg("ex").noconvert(true), py::arg("ey").noconvert(true), py::arg("ez").noconvert(true), + py::arg("bx").noconvert(true), py::arg("by").noconvert(true), py::arg("bz").noconvert(true), + py::arg("volume"), py::arg("dt"), + py::arg("ref_quantity") = py::float_(1.0) + ); + //________________________________________ + +} + +// ______________________________________________________________________________________________