diff --git a/.github/workflows/cleanup-cache.yml b/.github/workflows/cleanup-cache.yml index 3abe232b879..c71a48cdad8 100644 --- a/.github/workflows/cleanup-cache.yml +++ b/.github/workflows/cleanup-cache.yml @@ -2,7 +2,7 @@ name: CleanUpCache on: workflow_run: - workflows: [🧹 clang-tidy, 🔍 CodeQL, 🐧 CUDA, 🐧 HIP, 🐧 Intel, 🍏 macOS, 🐧 OpenMP] + workflows: [🧴 clang sanitizers, 🧹 clang-tidy, 🔍 CodeQL, 🐧 CUDA, 🐧 HIP, 🐧 Intel, 🍏 macOS, 🐧 OpenMP] types: - completed diff --git a/.github/workflows/cuda.yml b/.github/workflows/cuda.yml index 11765013bb7..a1e7f5affda 100644 --- a/.github/workflows/cuda.yml +++ b/.github/workflows/cuda.yml @@ -131,7 +131,7 @@ jobs: which nvcc || echo "nvcc not in PATH!" git clone https://github.com/AMReX-Codes/amrex.git ../amrex - cd ../amrex && git checkout --detach 216ce6f37de4b65be57fc1006b3457b4fc318e03 && cd - + cd ../amrex && git checkout --detach 028638564f7be0694b9898f8d4088cdbf9a6f9f5 && cd - make COMP=gcc QED=FALSE USE_MPI=TRUE USE_GPU=TRUE USE_OMP=FALSE USE_FFT=TRUE USE_CCACHE=TRUE -j 4 ccache -s diff --git a/.github/workflows/insitu.yml b/.github/workflows/insitu.yml index be93dfb9beb..35f16842935 100644 --- a/.github/workflows/insitu.yml +++ b/.github/workflows/insitu.yml @@ -101,8 +101,8 @@ jobs: cmake --build build -j 10 - name: 2D Test run: | - cp Examples/Tests/ionization/inputs_test_2d_ionization_lab . - cp Examples/Tests/ionization/catalyst_pipeline.py . + cp Examples/Tests/field_ionization/inputs_test_2d_ionization_lab . + cp Examples/Tests/field_ionization/catalyst_pipeline.py . mpiexec -n 2 ./build/bin/warpx.2d \ inputs_test_2d_ionization_lab \ catalyst.script_paths = catalyst_pipeline.py\ diff --git a/Docs/source/developers/developers.rst b/Docs/source/developers/developers.rst index aa2e6196377..222f665b563 100644 --- a/Docs/source/developers/developers.rst +++ b/Docs/source/developers/developers.rst @@ -15,7 +15,6 @@ Implementation Details initialization diagnostics moving_window - qed portability warning_logger python diff --git a/Docs/source/developers/testing.rst b/Docs/source/developers/testing.rst index 8b85976c6f0..5bbc7d0fef4 100644 --- a/Docs/source/developers/testing.rst +++ b/Docs/source/developers/testing.rst @@ -81,6 +81,24 @@ For easier debugging, it can be convenient to run the tests on your local machin ctest --test-dir build -E laser_acceleration +* Sometimes two or more tests share a large number of input parameters and differ by a small set of options. + Such tests typically also share a base string in their names. + For example, you can find three different tests named ``test_3d_langmuir_multi``, ``test_3d_langmuir_multi_nodal`` and ``test_3d_langmuir_multi_picmi``. + In such a case, if you wish to run the test ``test_3d_langmuir_multi`` only, this can be done again with the ``-R`` regular `expression filter `__ via + + .. code-block:: sh + + ctest --test-dir build -R "test_3d_langmuir_multi\..*" + + Note that filtering with ``-R "test_3d_langmuir_multi"`` would include the additional tests that have the same substring in their name and would not be sufficient to isolate a single test. + Note also that the escaping ``\.`` in the regular expression is necessary in order to take into account the fact that each test is automatically appended with the strings ``.run``, ``.analysis`` and possibly ``.cleanup``. + +* Run only tests not labeled with the ``slow`` label: + + .. code-block:: sh + + ctest --test-dir build -LE slow + Once the execution of CTest is completed, you can find all files associated with each test in its corresponding directory under ``build/bin/``. For example, if you run the single test ``test_3d_laser_acceleration``, you can find all files associated with this test in the directory ``build/bin/test_3d_laser_acceleration/``. @@ -155,7 +173,7 @@ A new test can be added by adding a corresponding entry in ``CMakeLists.txt`` as If you need a new Python package dependency for testing, please add it in `Regression/requirements.txt `__. -Sometimes, two tests share a large number of input parameters. The shared input parameters can be collected in a "base" input file that can be passed as a runtime parameter in the actual test input files through the parameter ``FILE``. +Sometimes two or more tests share a large number of input parameters. The shared input parameters can be collected in a "base" input file that can be passed as a runtime parameter in the actual test input files through the parameter ``FILE``. Naming conventions for automated tests -------------------------------------- diff --git a/Docs/source/index.rst b/Docs/source/index.rst index 6216ad52290..9668620976a 100644 --- a/Docs/source/index.rst +++ b/Docs/source/index.rst @@ -111,8 +111,7 @@ Theory theory/amr theory/boundary_conditions theory/boosted_frame - theory/input_output - theory/collisions + theory/multiphysics_extensions theory/kinetic_fluid_hybrid_model theory/cold_fluid_model diff --git a/Docs/source/latex_theory/allbibs.bib b/Docs/source/latex_theory/allbibs.bib index b3475c5a81b..62810ca5d6a 100644 --- a/Docs/source/latex_theory/allbibs.bib +++ b/Docs/source/latex_theory/allbibs.bib @@ -2187,3 +2187,45 @@ @book{godfrey1985iprop title = {{The IPROP Three-Dimensional Beam Propagation Code}}, year = {1985} } + +@article{Ammosov1986, +title = {Tunnel ionization of complex atoms and of atomic ions in an alternating electromagnetic field}, +volume = {64}, +issn = {0044-4510}, +doi = {10.1117/12.938695}, +number = {December 1986}, +journal = {Sov. Phys. JETP}, +author = {Ammosov, M. V. and Delone, N. B. and Krainov, V. P.}, +year = {1986}, +pmid = {22232002}, +note = {ISBN: 0892526998}, +pages = {1191--1194}, +} + + +@article{zhang_empirical_2014, +title = {Empirical formula for over-barrier strong-field ionization}, +volume = {90}, +issn = {1050-2947, 1094-1622}, +doi = {10.1103/PhysRevA.90.043410}, +language = {en}, +number = {4}, +journal = {Physical Review A}, +author = {Zhang, Qingbin and Lan, Pengfei and Lu, Peixiang}, +month = oct, +year = {2014}, +pages = {043410}, +} + + +@book{Mulser2010, +title = {High {Power} {Laser}-{Matter} {Interaction}}, +volume = {238}, +isbn = {978-3-540-50669-0}, +publisher = {Springer Berlin Heidelberg}, +author = {Mulser, Peter and Bauer, Dieter}, +year = {2010}, +pmid = {25246403}, +doi = {10.1007/978-3-540-46065-7}, +note = {Series Title: Springer Tracts in Modern Physics}, +} diff --git a/Docs/source/theory/boosted_frame.rst b/Docs/source/theory/boosted_frame.rst index ea1f662bd30..04c30bedc98 100644 --- a/Docs/source/theory/boosted_frame.rst +++ b/Docs/source/theory/boosted_frame.rst @@ -12,7 +12,21 @@ The simulations of plasma accelerators from first principles are extremely compu A first principle simulation of a short driver beam (laser or charged particles) propagating through a plasma that is orders of magnitude longer necessitates a very large number of time steps. Recasting the simulation in a frame of reference that is moving close to the speed of light in the direction of the driver beam leads to simulating a driver beam that appears longer propagating through a plasma that appears shorter than in the laboratory. Thus, this relativistic transformation of space and time reduces the disparity of scales, and thereby the number of time steps to complete the simulation, by orders of magnitude. -Even using a moving window, however, a full PIC simulation of a plasma accelerator can be extraordinarily demanding computationally, as many time steps are needed to resolve the crossing of the short driver beam with the plasma column. As it turns out, choosing an optimal frame of reference that travels close to the speed of light in the direction of the laser or particle beam (as opposed to the usual choice of the laboratory frame) enables speedups by orders of magnitude :cite:p:`bf-Vayprl07,bf-Vaypop2011`. This is a result of the properties of Lorentz contraction and dilation of space and time. In the frame of the laboratory, a very short driver (laser or particle) beam propagates through a much longer plasma column, necessitating millions to tens of millions of time steps for parameters in the range of the BELLA or FACET-II experiments. As sketched in :numref:`fig_Boosted_frame`, in a frame moving with the driver beam in the plasma at velocity :math:`v=\beta c` (where :math:`c` is the speed of light in vacuum), the beam length is now elongated by :math:`\approx(1+\beta)\gamma` while the plasma contracts by :math:`\gamma` (where :math:`\gamma=1/\sqrt{1-\beta^2}` is the relativistic factor associated with the frame velocity). The number of time steps that is needed to simulate a “longer” beam through a “shorter” plasma is now reduced by up to :math:`\approx(1+\beta) \gamma^2` (a detailed derivation of the speedup is given below). +Even using a moving window, however, a full PIC simulation of a plasma accelerator can be extraordinarily demanding computationally, as many time steps are needed to resolve the crossing of the short driver beam with the plasma column. +As it turns out, choosing an optimal frame of reference that travels close to the speed of light in the direction of the laser or particle beam (as opposed to the usual choice of the laboratory frame) enables speedups by orders of magnitude :cite:p:`bf-Vayprl07,bf-Vaypop2011`. +This is a result of the properties of Lorentz contraction and dilation of space and time. +In the frame of the laboratory, a very short driver (laser or particle) beam propagates through a much longer plasma column, necessitating millions to tens of millions of time steps for parameters in the range of the BELLA or FACET-II experiments. +As sketched in :numref:`fig_Boosted_frame`, in a frame moving with the driver beam in the plasma at velocity :math:`v=\beta c` (where :math:`c` is the speed of light in vacuum), the beam length is now elongated by :math:`\approx(1+\beta)\gamma` while the plasma contracts by :math:`\gamma` (where :math:`\gamma=1/\sqrt{1-\beta^2}` is the relativistic factor associated with the frame velocity) +The number of time steps that is needed to simulate a “longer” beam through a “shorter” plasma is now reduced by up to :math:`\approx(1+\beta) \gamma^2` (a detailed derivation of the speedup is given below). + +.. note:: + + For additional reading on inputs and outputs in boosted frame simulations, consider the following pages: + + .. toctree:: + :maxdepth: 1 + + boosted_frame/input_output The modeling of a plasma acceleration stage in a boosted frame involves the fully electromagnetic modeling of a plasma propagating at near the speed of light, for which Numerical Cerenkov diff --git a/Docs/source/theory/Input_output.png b/Docs/source/theory/boosted_frame/Input_output.png similarity index 100% rename from Docs/source/theory/Input_output.png rename to Docs/source/theory/boosted_frame/Input_output.png diff --git a/Docs/source/theory/input_output.rst b/Docs/source/theory/boosted_frame/input_output.rst similarity index 99% rename from Docs/source/theory/input_output.rst rename to Docs/source/theory/boosted_frame/input_output.rst index 21a5f5c8d2c..f47915116df 100644 --- a/Docs/source/theory/input_output.rst +++ b/Docs/source/theory/boosted_frame/input_output.rst @@ -1,4 +1,4 @@ -.. _theory-io: +.. _boosted_frame-io: Inputs and Outputs ================== diff --git a/Docs/source/theory/collisions.rst b/Docs/source/theory/multiphysics/collisions.rst similarity index 98% rename from Docs/source/theory/collisions.rst rename to Docs/source/theory/multiphysics/collisions.rst index ef6b83a699b..08485345a13 100644 --- a/Docs/source/theory/collisions.rst +++ b/Docs/source/theory/multiphysics/collisions.rst @@ -1,4 +1,4 @@ -.. _theory-collisions: +.. _multiphysics-collisions: Collisions ========== @@ -8,7 +8,7 @@ including collisions between kinetic particles (Coulomb collisions, DSMC, nuclear fusion) as well as collisions between kinetic particles and a fixed (i.e. non-evolving) background species (MCC, background stopping). -.. _theory-collisions-mcc: +.. _multiphysics-collisions-mcc: Background Monte Carlo Collisions (MCC) --------------------------------------- @@ -52,7 +52,7 @@ for all scattering processes are evaluated at the energy as calculated above. Once a particle is selected for a specific collision process, that process determines how the particle is scattered as outlined below. -.. _theory-collisions-dsmc: +.. _multiphysics-collisions-dsmc: Direct Simulation Monte Carlo (DSMC) ------------------------------------ diff --git a/Docs/source/theory/multiphysics/ionization.rst b/Docs/source/theory/multiphysics/ionization.rst new file mode 100644 index 00000000000..11abea386c8 --- /dev/null +++ b/Docs/source/theory/multiphysics/ionization.rst @@ -0,0 +1,73 @@ +.. _multiphysics-ionization: + +Ionization +========== + +Field Ionization +---------------- + +Under the influence of a sufficiently strong external electric field atoms become ionized. +Particularly the dynamics of interactions between ultra-high intensity laser pulses and matter, e.g., Laser-Plasma Acceleration (LPA) with ionization injection, or Laser-Plasma Interactions with solid density targets (LPI) can depend on field ionization dynamics as well. + +WarpX models field ionization based on a description of the Ammosov-Delone-Krainov model:cite:p:`mpion-Ammosov1986` following :cite:t:`mpion-ChenPRSTAB13`. + +Implementation Details and Assumptions +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +.. note:: + + The current implementation makes the following assumptions + + * Energy for ionization processes is not removed from the electromagnetic fields + * Only one single-level ionization process can occur per macroparticle and time step + * Ionization happens at the beginning of the PIC loop before the field solve + * Angular momentum quantum number :math:`l = 0` and magnetic quantum number :math:`m = 0` + +The model implements the following equations (assumptions to :math:`l` and :math:`m` have already been applied). + +The electric field amplitude is calculated in the particle's frame of reference. + +.. math:: + + \begin{aligned} + \vec{E}_\mathrm{dc} &= \sqrt{ - \frac{1}{\mathrm{c}^2} \left( \vec{u} \cdot \vec{E} \right)^2 + + \left( \gamma \vec{E} + \vec{u} \times \vec{B} \right)^2 } + \\ + \gamma &= \sqrt{1 + \frac{\vec{u}^2}{\mathrm{c}^2}} + \end{aligned} + +Here, :math:`\vec{u} = (u_x, u_y, u_z)` is the momentum normalized to the particle mass, :math:`u_i = (\beta \gamma)_i \mathrm{c}`. +:math:`E_\mathrm{dc} = |\vec{E}_\mathrm{dc}|` is the DC-field in the frame of the particle. + +.. math:: + + \begin{aligned} + P &= 1 - \mathrm{e}^{-W\mathrm{d}\tau/\gamma} + \\ + W &= \omega_\mathrm{a} \mathcal{C}^2_{n^* l^*} \frac{U_\mathrm{ion}}{2 U_H} + \left[ 2 \frac{E_\mathrm{a}}{E_\mathrm{dc}} \left( \frac{U_\mathrm{ion}}{U_\mathrm{H}} \right)^{3/2} \right]^{2n^*-1} + \times \exp\left[ - \frac{2}{3} \frac{E_\mathrm{a}}{E_\mathrm{dc}} \left( \frac{U_\mathrm{ion}}{U_\mathrm{H}} \right)^{3/2} \right] + \\ + \mathcal{C}^2_{n^* l^*} &= \frac{2^{2n^*}}{n^* \Gamma(n^* + l^* + 1) \Gamma(n^* - l^*)} + \end{aligned} + +where :math:`\mathrm{d}\tau` is the simulation timestep, which is divided by the particle :math:`\gamma` to account for time dilation. The quantities are: :math:`\omega_\mathrm{a}`, the atomic unit frequency, :math:`U_\mathrm{ion}`, the ionization potential, :math:`U_\mathrm{H}`, Hydrogen ground state ionization potential, :math:`E_\mathrm{a}`, the atomic unit electric field, :math:`n^* = Z \sqrt{U_\mathrm{H}/U_\mathrm{ion}}`, the effective principal quantum number (*Attention!* :math:`Z` is the ionization state *after ionization*.) , :math:`l^* = n_0^* - 1`, the effective orbital quantum number. + +Empirical Extension to Over-the-Barrier Regime for Hydrogen +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +For hydrogen, WarpX offers the modified empirical ADK extension to the Over-the-Barrier (OTB) published in :cite:t:`mpion-zhang_empirical_2014` Eq. (8). + +.. math:: + + W_\mathrm{M} = \exp\left[ -\left( a_1 \frac{E^2}{E_\mathrm{b}} + a_2 \frac{E}{E_\mathrm{b}} + a_3 \right) \right] W_\mathrm{ADK} + +The parameters :math:`a_1` through :math:`a_3` are independent of :math:`E` and can be found in the same reference. :math:`E_\mathrm{b}` is the classical Barrier Suppresion Ionization (BSI) field strength :math:`E_\mathrm{b} = U_\mathrm{ion}^2 / (4 Z)` given here in atomic units (AU). For a detailed description of conversion between unit systems consider the book by :cite:t:`mpion-Mulser2010`. + +Testing +^^^^^^^ + +* `Testing the field ionization module <../../../../Examples/Tests/field_ionization/README.rst>`_. + +.. bibliography:: + :keyprefix: mpion- diff --git a/Docs/source/developers/qed.rst b/Docs/source/theory/multiphysics/qed.rst similarity index 86% rename from Docs/source/developers/qed.rst rename to Docs/source/theory/multiphysics/qed.rst index f509d6ea386..3e961bb2898 100644 --- a/Docs/source/developers/qed.rst +++ b/Docs/source/theory/multiphysics/qed.rst @@ -1,7 +1,7 @@ -.. _developers-qed: +.. _multiphysics-qed: -QED -==================== +Quantum Electrodynamics (QED) +============================= Quantum synchrotron ------------------- @@ -28,7 +28,8 @@ electron-positron pairs can be created in vacuum in the function ``MultiParticleContainer::doQEDSchwinger`` in turn calls the function ``filterCreateTransformFromFAB``: -.. doxygenfunction:: filterCreateTransformFromFAB(DstTile&, DstTile&, const amrex::Box, const FABs&, const Index, const Index, FilterFunc&&, CreateFunc1&&, CreateFunc2&&, TransFunc&&) +Filter Create Transform Function +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ``filterCreateTransformFromFAB`` proceeds in three steps. In the filter phase, we loop on every cell and calculate the number of physical pairs created within diff --git a/Docs/source/theory/multiphysics_extensions.rst b/Docs/source/theory/multiphysics_extensions.rst new file mode 100644 index 00000000000..6c9ab040ef2 --- /dev/null +++ b/Docs/source/theory/multiphysics_extensions.rst @@ -0,0 +1,13 @@ +.. _theory-multiphysics: + +Multi-Physics Extensions +======================== + +WarpX includes various extensions to the traditional PIC loop which enable it to model additional physics. + +.. toctree:: + :maxdepth: 1 + + multiphysics/collisions + multiphysics/ionization + multiphysics/qed diff --git a/Docs/source/usage/examples.rst b/Docs/source/usage/examples.rst index 244fbda6f75..f1bd2ec4266 100644 --- a/Docs/source/usage/examples.rst +++ b/Docs/source/usage/examples.rst @@ -137,6 +137,11 @@ An example of initializing the fields by accessing their data through Python, ad Many Further Examples, Demos and Tests -------------------------------------- +.. toctree:: + :maxdepth: 1 + + examples/field_ionization/README.rst + WarpX runs over 200 integration tests on a variety of modeling cases, which validate and demonstrate its functionality. Please see the `Examples/Tests/ `__ directory for many more examples. diff --git a/Docs/source/usage/examples/field_ionization b/Docs/source/usage/examples/field_ionization new file mode 120000 index 00000000000..b1c3e38dab2 --- /dev/null +++ b/Docs/source/usage/examples/field_ionization @@ -0,0 +1 @@ +../../../../Examples/Tests/field_ionization/ \ No newline at end of file diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 1d9c0c14bbf..86ab7594c5f 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -821,7 +821,7 @@ Particle initialization * ``particles.rigid_injected_species`` (`strings`, separated by spaces) List of species injected using the rigid injection method. The rigid injection method is useful when injecting a relativistic particle beam in boosted-frame - simulations; see the :ref:`input-output section ` for more details. + simulations; see the :ref:`input-output section ` for more details. For species injected using this method, particles are translated along the `+z` axis with constant velocity as long as their ``z`` coordinate verifies ``zzinject_plane``, @@ -1953,7 +1953,7 @@ Collision models ---------------- WarpX provides several particle collision models, using varying degrees of approximation. -Details about the collision models can be found in the :ref:`theory section `. +Details about the collision models can be found in the :ref:`theory section `. * ``collisions.collision_names`` (`strings`, separated by spaces) The name of each collision type. @@ -1976,10 +1976,10 @@ Details about the collision models can be found in the :ref:`theory section .species_type = 'deuterium'``) - ``dsmc`` for pair-wise, non-Coulomb collisions between kinetic species. This is a "direct simulation Monte Carlo" treatment of collisions between - kinetic species. See :ref:`DSMC section `. + kinetic species. See :ref:`DSMC section `. - ``background_mcc`` for collisions between particles and a neutral background. This is a relativistic Monte Carlo treatment for particles colliding - with a neutral background gas. See :ref:`MCC section `. + with a neutral background gas. See :ref:`MCC section `. - ``background_stopping`` for slowing of ions due to collisions with electrons or ions. This implements the approximate formulae as derived in Introduction to Plasma Physics, from Goldston and Rutherford, section 14.2. diff --git a/Examples/CMakeLists.txt b/Examples/CMakeLists.txt index 14943c93548..8bb02088a16 100644 --- a/Examples/CMakeLists.txt +++ b/Examples/CMakeLists.txt @@ -218,6 +218,22 @@ function(add_warpx_test endif() endfunction() +# Add a CTest label to a WarpX test set. +# +# Labeling it here will add the label to the run test, its analysis and its cleanup. +# +# name: unique name of this test +# label: ctest LABELS property value to be added +# +function(label_warpx_test name label) + set(_test_names "${name}.run;${name}.analysis;${name}.cleanup") + foreach(_test_name IN LISTS _test_names) + if(TEST ${_test_name}) + set_property(TEST ${_test_name} APPEND PROPERTY LABELS "${label}") + endif() + endforeach() +endfunction() + # Add tests (alphabetical order) ############################################## # diff --git a/Examples/Physics_applications/beam_beam_collision/CMakeLists.txt b/Examples/Physics_applications/beam_beam_collision/CMakeLists.txt index 0b34eeff865..09e96f04d7f 100644 --- a/Examples/Physics_applications/beam_beam_collision/CMakeLists.txt +++ b/Examples/Physics_applications/beam_beam_collision/CMakeLists.txt @@ -10,3 +10,4 @@ add_warpx_test( diags/diag1/ # output OFF # dependency ) +label_warpx_test(test_3d_beam_beam_collision slow) diff --git a/Examples/Physics_applications/beam_beam_collision/README.rst b/Examples/Physics_applications/beam_beam_collision/README.rst index 4f89365c8f0..a7a06521218 100644 --- a/Examples/Physics_applications/beam_beam_collision/README.rst +++ b/Examples/Physics_applications/beam_beam_collision/README.rst @@ -39,7 +39,7 @@ We compare different results: * (blue) large-scale WarpX simulation (high resolution and ad hoc generated tables ; * (black) literature results from :cite:t:`ex-Yakimenko2019`. -The small-scale simulation has been performed with a resolution of ``nx = 64, ny = 64, nz = 128`` grid cells, while the large-scale one has a much higher resolution of ``nx = 512, ny = 512, nz = 1024``. Moreover, the large-scale simulation uses dedicated QED lookup tables instead of the builtin tables. To generate the tables within WarpX, the code must be compiled with the flag ``-DWarpX_QED_TABLE_GEN=ON``. For the large-scale simulation we have used the following options: +The small-scale simulation has been performed with a resolution of ``nx = 64, ny = 64, nz = 64`` grid cells, while the large-scale one has a much higher resolution of ``nx = 512, ny = 512, nz = 1024``. Moreover, the large-scale simulation uses dedicated QED lookup tables instead of the builtin tables. To generate the tables within WarpX, the code must be compiled with the flag ``-DWarpX_QED_TABLE_GEN=ON``. For the large-scale simulation we have used the following options: .. code-block:: ini @@ -63,7 +63,7 @@ The small-scale simulation has been performed with a resolution of ``nx = 64, ny qed_bw.tab_pair_frac_how_many=512 qed_bw.save_table_in=my_bw_table.txt -.. figure:: https://user-images.githubusercontent.com/17280419/291749626-aa61fff2-e6d2-45a3-80ee-84b2851ea0bf.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTEiLCJleHAiOjE3MDMwMzQzNTEsIm5iZiI6MTcwMzAzNDA1MSwicGF0aCI6Ii8xNzI4MDQxOS8yOTE3NDk2MjYtYWE2MWZmZjItZTZkMi00NWEzLTgwZWUtODRiMjg1MWVhMGJmLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFJV05KWUFYNENTVkVINTNBJTJGMjAyMzEyMjAlMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjMxMjIwVDAxMDA1MVomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPWFiYzY2MGQyYzIyZGIzYzUxOWI3MzNjZTk5ZDM1YzgyNmY4ZDYxOGRlZjAyZTIwNTAyMTc3NTgwN2Q0YjEwNGMmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.I96LQpjqmFXirPDVnBlFQIkCuenR6IuOSY0OIIQvtCo +.. figure:: https://gist.github.com/user-attachments/assets/2dd43782-d039-4faa-9d27-e3cf8fb17352 :alt: Beam-beam collision benchmark against :cite:t:`ex-Yakimenko2019`. :width: 100% diff --git a/Examples/Physics_applications/beam_beam_collision/inputs_test_3d_beam_beam_collision b/Examples/Physics_applications/beam_beam_collision/inputs_test_3d_beam_beam_collision index 488e997f895..e856a078003 100644 --- a/Examples/Physics_applications/beam_beam_collision/inputs_test_3d_beam_beam_collision +++ b/Examples/Physics_applications/beam_beam_collision/inputs_test_3d_beam_beam_collision @@ -27,8 +27,7 @@ my_constants.Lz = 180.0*clight/omegab # for a full scale simulation use: nx, ny, nz = 512, 512, 1024 my_constants.nx = 64 my_constants.ny = 64 -my_constants.nz = 128 - +my_constants.nz = 64 # TIME my_constants.T = 0.7*Lz/clight diff --git a/Examples/Physics_applications/laser_acceleration/CMakeLists.txt b/Examples/Physics_applications/laser_acceleration/CMakeLists.txt index 1a09a669a6d..46e97a53d54 100644 --- a/Examples/Physics_applications/laser_acceleration/CMakeLists.txt +++ b/Examples/Physics_applications/laser_acceleration/CMakeLists.txt @@ -11,16 +11,6 @@ add_warpx_test( OFF # dependency ) -add_warpx_test( - test_1d_laser_acceleration_fluid # name - 1 # dims - 2 # nprocs - inputs_test_1d_laser_acceleration_fluid # inputs - analysis_1d_fluid.py # analysis - diags/diag1040000 # output - OFF # dependency -) - add_warpx_test( test_1d_laser_acceleration_fluid_boosted # name 1 # dims @@ -30,6 +20,7 @@ add_warpx_test( diags/diag1000001 # output OFF # dependency ) +label_warpx_test(test_1d_laser_acceleration_fluid_boosted slow) add_warpx_test( test_1d_laser_acceleration_picmi # name @@ -140,3 +131,4 @@ add_warpx_test( diags/diag1000010 # output OFF # dependency ) +label_warpx_test(test_rz_laser_acceleration_picmi slow) diff --git a/Examples/Physics_applications/laser_acceleration/analysis_1d_fluid.py b/Examples/Physics_applications/laser_acceleration/analysis_1d_fluid.py deleted file mode 100755 index 593036bc3f6..00000000000 --- a/Examples/Physics_applications/laser_acceleration/analysis_1d_fluid.py +++ /dev/null @@ -1,193 +0,0 @@ -#!/usr/bin/env python3 - -# Copyright 2019-2023 Grant Johnson, Remi Lehe -# -# -# This file is part of WarpX. -# -# License: BSD-3-Clause-LBNL -# -# This is a script that analyses the simulation results from -# the script `inputs_1d`. This simulates a 1D WFA with Pondermotive Envelope: -# REF: (Equations 20-23) https://journals.aps.org/rmp/pdf/10.1103/RevModPhys.81.1229 -import os -import sys - -import matplotlib - -matplotlib.use("Agg") -import matplotlib.pyplot as plt -import yt - -yt.funcs.mylog.setLevel(50) - -import numpy as np -from scipy.constants import c, e, epsilon_0, m_e - -sys.path.insert(1, "../../../../warpx/Regression/Checksum/") -import checksumAPI - -# this will be the name of the plot file -fn = sys.argv[1] - -# Parameters (these parameters must match the parameters in `inputs.multi.rt`) -n0 = 20.0e23 -# Plasma frequency -wp = np.sqrt((n0 * e**2) / (m_e * epsilon_0)) -kp = wp / c -tau = 15.0e-15 -a0 = 2.491668 -e = -e # Electrons -lambda_laser = 0.8e-6 - -zmin = -20e-6 -zmax = 100.0e-6 -Nz = 10240 - -# Compute the theory - -# Call the ode solver -from scipy.integrate import odeint - - -# ODE Function -def odefcn(phi, xi, kp, a0, c, tau, xi_0, lambda_laser): - phi1, phi2 = phi - a_sq = ( - a0**2 - * np.exp(-2 * (xi - xi_0) ** 2 / (c**2 * tau**2)) - * np.sin(2 * np.pi * (xi - xi_0) / lambda_laser) ** 2 - ) - dphi1_dxi = phi2 - dphi2_dxi = kp**2 * ((1 + a_sq) / (2 * (1 + phi1) ** 2) - 0.5) - return [dphi1_dxi, dphi2_dxi] - - -# Call odeint to solve the ODE -xi_span = [-20e-6, 100e-6] -xi_0 = 0e-6 -phi0 = [0.0, 0.0] -dxi = (zmax - zmin) / Nz -xi = zmin + dxi * (0.5 + np.arange(Nz)) -phi = odeint(odefcn, phi0, xi, args=(kp, a0, c, tau, xi_0, lambda_laser)) - -# Change array direction to match the simulations -xi = -xi[::-1] -phi = phi[::-1] -xi_0 = -0e-6 -phi2 = phi[:, 0] -Ez = -phi[:, 1] - -# Compute the derived quantities -a_sq = ( - a0**2 - * np.exp(-2 * (xi - xi_0) ** 2 / (c**2 * tau**2)) - * np.sin(2 * np.pi * (xi - xi_0) / lambda_laser) ** 2 -) -gamma_perp_sq = 1 + a_sq -n = n0 * (gamma_perp_sq + (1 + phi2) ** 2) / (2 * (1 + phi2) ** 2) -uz = (gamma_perp_sq - (1 + phi2) ** 2) / (2 * (1 + phi2)) -gamma = (gamma_perp_sq + (1 + phi2) ** 2) / (2 * (1 + phi2)) - -# Theory Components [convert to si] -uz *= c -J_th = np.multiply(np.divide(uz, gamma), n) -J_th *= e -rho_th = e * n -E_th = Ez -E_th *= (m_e * c * c) / e -V_th = np.divide(uz, gamma) -V_th /= c -# Remove the ions -rho_th = rho_th - e * n0 - -# Dictate which region to compare solutions over -# (Currently this is the full domain) -min_i = 0 -max_i = 10240 - -# Read the file -ds = yt.load(fn) -t0 = ds.current_time.to_value() -data = ds.covering_grid( - level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions -) -# Check the validity of the fields -error_rel = 0 -for field in ["Ez"]: - E_sim = data[("mesh", field)].to_ndarray()[:, 0, 0] - # E_th = get_theoretical_field(field, t0) - max_error = ( - abs(E_sim[min_i:max_i] - E_th[min_i:max_i]).max() / abs(E_th[min_i:max_i]).max() - ) - print("%s: Max error: %.2e" % (field, max_error)) - error_rel = max(error_rel, max_error) - -# Check the validity of the currents -for field in ["Jz"]: - J_sim = data[("mesh", field)].to_ndarray()[:, 0, 0] - # J_th = get_theoretical_J_field(field, t0) - max_error = ( - abs(J_sim[min_i:max_i] - J_th[min_i:max_i]).max() / abs(J_th[min_i:max_i]).max() - ) - print("%s: Max error: %.2e" % (field, max_error)) - error_rel = max(error_rel, max_error) - -# Check the validity of the charge -for field in ["rho"]: - rho_sim = data[("boxlib", field)].to_ndarray()[:, 0, 0] - # rho_th = get_theoretical_rho_field(field, t0) - max_error = ( - abs(rho_sim[min_i:max_i] - rho_th[min_i:max_i]).max() - / abs(rho_th[min_i:max_i]).max() - ) - print("%s: Max error: %.2e" % (field, max_error)) - error_rel = max(error_rel, max_error) - -V_sim = np.divide(J_sim, rho_sim) -V_sim /= c - -# Create a figure with 2 rows and 2 columns -fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12, 8)) - -# Titles and labels -titles = ["Ez", "rho", "Jz", "Vz/c"] -xlabel = r"Xi" -ylabel = ["Ez", "rho", "Jz", "Vz/c"] - -# Plotting loop -for i in range(3): - ax = axes[i // 2, i % 2] # Get the current subplot - - # Plot theoretical data - ax.plot(xi, [E_th, rho_th, J_th, V_th][i], label="Theoretical") - - # Plot simulated data - ax.plot(xi, [E_sim, rho_sim, J_sim, V_sim][i], label="Simulated") - - # Set titles and labels - ax.set_title(f"{titles[i]} vs Xi") - ax.set_xlabel(xlabel) - ax.set_ylabel(ylabel[i]) - - # Add legend - ax.legend() - -# Adjust subplot layout -plt.tight_layout() - -# Save the figure -plt.savefig("wfa_fluid_nonlinear_1d_analysis.png") - -plt.show() - - -tolerance_rel = 0.20 - -print("error_rel : " + str(error_rel)) -print("tolerance_rel: " + str(tolerance_rel)) - -assert error_rel < tolerance_rel - -test_name = os.path.split(os.getcwd())[1] -checksumAPI.evaluate_checksum(test_name, fn) diff --git a/Examples/Physics_applications/laser_acceleration/inputs_test_1d_laser_acceleration_fluid b/Examples/Physics_applications/laser_acceleration/inputs_test_1d_laser_acceleration_fluid deleted file mode 100644 index 73fa6b7283f..00000000000 --- a/Examples/Physics_applications/laser_acceleration/inputs_test_1d_laser_acceleration_fluid +++ /dev/null @@ -1,72 +0,0 @@ -################################# -####### GENERAL PARAMETERS ###### -################################# -max_step = 40000 -amr.n_cell = 10240 -amr.max_grid_size = 512 # maximum size of each AMReX box, used to decompose the domain -amr.blocking_factor = 512 # minimum size of each AMReX box, used to decompose the domain -geometry.dims = 1 -geometry.prob_lo = -120.e-6 # physical domain -geometry.prob_hi = 0.e-6 -amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported) - -################################# -####### Boundary condition ###### -################################# -boundary.field_lo = pec -boundary.field_hi = pec - -################################# -############ NUMERICS ########### -################################# -warpx.verbose = 1 -warpx.do_dive_cleaning = 0 -warpx.use_filter = 0 -warpx.cfl = 0.45 #Fluid CFL < 0.5 -warpx.do_moving_window = 1 -warpx.moving_window_dir = z -warpx.moving_window_v = 1.0 # units of speed of light -warpx.do_dynamic_scheduling = 0 -warpx.serialize_initial_conditions = 1 - -################################# -############ PLASMA ############# -################################# -fluids.species_names = electrons ions - -electrons.species_type = electron -electrons.profile = parse_density_function -electrons.density_function(x,y,z) = "1.0e10 + 20.e23*((z*5.e4 + -0.5)*(z>10.e-6)*(z<30.e-6)) + 20.e23*((z>30.e-6))" -electrons.momentum_distribution_type = "at_rest" - -ions.charge = q_e -ions.mass = m_p -ions.profile = parse_density_function -ions.density_function(x,y,z) = "1.0e10 + 20.e23*((z*5.e4 + -0.5)*(z>10.e-6)*(z<30.e-6)) + 20.e23*((z>30.e-6))" -ions.momentum_distribution_type = "at_rest" - -# Order of particle shape factors -algo.particle_shape = 3 - -################################# -############ LASER ############## -################################# -lasers.names = laser1 -laser1.profile = Gaussian -laser1.position = 0. 0. -11.e-6 # This point is on the laser plane -laser1.direction = 0. 0. 1. # The plane normal direction -laser1.polarization = 0. 1. 0. # The main polarization vector -laser1.e_max = 10.e12 # Maximum amplitude of the laser field (in V/m) -laser1.profile_waist = 5.e-6 # The waist of the laser (in m) -laser1.profile_duration = 15.e-15 # The duration of the laser (in s) -laser1.profile_t_peak = 30.e-15 # Time at which the laser reaches its peak (in s) -laser1.profile_focal_distance = 100.e-6 # Focal distance from the antenna (in m) -laser1.wavelength = 0.8e-6 # The wavelength of the laser (in m) - -# Diagnostics -diagnostics.diags_names = diag1 - -# LAB -diag1.intervals = 20000 -diag1.diag_type = Full -diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho diff --git a/Examples/Tests/CMakeLists.txt b/Examples/Tests/CMakeLists.txt index 108a28a6539..6fea9368e78 100644 --- a/Examples/Tests/CMakeLists.txt +++ b/Examples/Tests/CMakeLists.txt @@ -25,7 +25,7 @@ add_subdirectory(gaussian_beam) add_subdirectory(implicit) add_subdirectory(initial_distribution) add_subdirectory(initial_plasma_profile) -add_subdirectory(ionization) +add_subdirectory(field_ionization) add_subdirectory(ion_stopping) add_subdirectory(langmuir) add_subdirectory(langmuir_fluids) diff --git a/Examples/Tests/ionization/CMakeLists.txt b/Examples/Tests/field_ionization/CMakeLists.txt similarity index 100% rename from Examples/Tests/ionization/CMakeLists.txt rename to Examples/Tests/field_ionization/CMakeLists.txt diff --git a/Examples/Tests/field_ionization/README.rst b/Examples/Tests/field_ionization/README.rst new file mode 100644 index 00000000000..c2f11ed1a40 --- /dev/null +++ b/Examples/Tests/field_ionization/README.rst @@ -0,0 +1,59 @@ +.. _examples-tests-field_ionization: + +Field Ionization +================ + +Run Test +-------- + +For `MPI-parallel `__ runs, prefix these lines with ``mpiexec -n 4 ...`` or ``srun -n 4 ...``, depending on the system. + +.. tab-set:: + + .. tab-item:: lab frame + + This example can be run **either** as: + + * **Python** script: ``python3 inputs_test_2d_ionization_picmi.py`` or + * WarpX **executable** using an input file: ``warpx.2d inputs_test_2d_ionization_lab max_step=1600`` + + .. tab-set:: + + .. tab-item:: Python: Script + + .. literalinclude:: inputs_test_2d_ionization_picmi.py + :language: python3 + :caption: You can copy this file from ``Examples/Tests/field_ionization/inputs_test_2d_ionization_picmi.py``. + + .. tab-item:: Executable: Input File + + .. literalinclude:: inputs_test_2d_ionization_lab + :language: ini + :caption: You can copy this file from ``Examples/Tests/field_ionization/inputs_test_2d_ionization_lab``. + + .. tab-item:: boosted frame + + This example can be run as: + + * WarpX **executable** using an input file: ``warpx.2d inputs_test_2d_ionization_boost max_step=420`` + + .. literalinclude:: inputs_test_2d_ionization_boost + :language: ini + :caption: You can copy this file from ``Examples/Tests/field_ionization/inputs_test_2d_ionization_boost``. + +Analyze +------- + +.. dropdown:: Script ``analysis.py`` + + .. literalinclude:: analysis.py + :language: python3 + :caption: You can copy this file from ``Examples/Tests/field_ionization/analysis.py``. + +Visualize +--------- + +.. figure:: https://gist.githubusercontent.com/johvandewetering/48d092c003915f1d1689b507caa2865b/raw/29f5d12ed77831047ca12f456a07dbf3b99770d5/image_ionization.png + :alt: Electric field of the laser pulse with (top) ions with ionization levels and (bottom) ionized electrons. + + Electric field of the laser pulse with (top) ions with ionization levels and (bottom) ionized electrons. diff --git a/Examples/Tests/ionization/analysis.py b/Examples/Tests/field_ionization/analysis.py similarity index 100% rename from Examples/Tests/ionization/analysis.py rename to Examples/Tests/field_ionization/analysis.py diff --git a/Examples/Tests/ionization/catalyst_pipeline.py b/Examples/Tests/field_ionization/catalyst_pipeline.py similarity index 100% rename from Examples/Tests/ionization/catalyst_pipeline.py rename to Examples/Tests/field_ionization/catalyst_pipeline.py diff --git a/Examples/Tests/ionization/inputs_test_2d_ionization_boost b/Examples/Tests/field_ionization/inputs_test_2d_ionization_boost similarity index 100% rename from Examples/Tests/ionization/inputs_test_2d_ionization_boost rename to Examples/Tests/field_ionization/inputs_test_2d_ionization_boost diff --git a/Examples/Tests/ionization/inputs_test_2d_ionization_lab b/Examples/Tests/field_ionization/inputs_test_2d_ionization_lab similarity index 100% rename from Examples/Tests/ionization/inputs_test_2d_ionization_lab rename to Examples/Tests/field_ionization/inputs_test_2d_ionization_lab diff --git a/Examples/Tests/ionization/inputs_test_2d_ionization_picmi.py b/Examples/Tests/field_ionization/inputs_test_2d_ionization_picmi.py similarity index 100% rename from Examples/Tests/ionization/inputs_test_2d_ionization_picmi.py rename to Examples/Tests/field_ionization/inputs_test_2d_ionization_picmi.py diff --git a/Examples/Tests/langmuir/CMakeLists.txt b/Examples/Tests/langmuir/CMakeLists.txt index 3f44d364276..bd0cea79c7a 100644 --- a/Examples/Tests/langmuir/CMakeLists.txt +++ b/Examples/Tests/langmuir/CMakeLists.txt @@ -341,6 +341,7 @@ if(WarpX_FFT) diags/diag1000040 # output OFF # dependency ) + label_warpx_test(test_3d_langmuir_multi_psatd_vay_deposition_nodal slow) endif() add_warpx_test( diff --git a/Examples/Tests/nci_psatd_stability/CMakeLists.txt b/Examples/Tests/nci_psatd_stability/CMakeLists.txt index ed087fc4190..f2b4ceae8ba 100644 --- a/Examples/Tests/nci_psatd_stability/CMakeLists.txt +++ b/Examples/Tests/nci_psatd_stability/CMakeLists.txt @@ -200,7 +200,8 @@ if(WarpX_FFT) 2 # nprocs inputs_test_rz_multiJ_psatd # inputs analysis_default_regression.py # analysis - diags/diag1000050 # output + diags/diag1000025 # output OFF # dependency ) + label_warpx_test(test_rz_multiJ_psatd slow) endif() diff --git a/Examples/Tests/nci_psatd_stability/inputs_test_rz_multiJ_psatd b/Examples/Tests/nci_psatd_stability/inputs_test_rz_multiJ_psatd index 5e263856256..6350b9aee51 100644 --- a/Examples/Tests/nci_psatd_stability/inputs_test_rz_multiJ_psatd +++ b/Examples/Tests/nci_psatd_stability/inputs_test_rz_multiJ_psatd @@ -1,8 +1,8 @@ # Iterations -max_step = 50 +max_step = 25 # Domain -amr.n_cell = 64 128 +amr.n_cell = 32 64 amr.max_level = 0 warpx.numprocs = 1 2 @@ -115,7 +115,7 @@ plasma_p.do_continuous_injection = 1 # Diagnostics diagnostics.diags_names = diag1 -diag1.intervals = 50 +diag1.intervals = 25 diag1.diag_type = Full diag1.fields_to_plot = Er Ez Bt jr jz rho rho_driver rho_plasma_e rho_plasma_p diag1.species = driver plasma_e plasma_p diff --git a/Examples/Tests/nodal_electrostatic/CMakeLists.txt b/Examples/Tests/nodal_electrostatic/CMakeLists.txt index 915298f15ab..a6b3f5b0102 100644 --- a/Examples/Tests/nodal_electrostatic/CMakeLists.txt +++ b/Examples/Tests/nodal_electrostatic/CMakeLists.txt @@ -10,3 +10,4 @@ add_warpx_test( diags/diag1000010 # output OFF # dependency ) +label_warpx_test(test_3d_nodal_electrostatic_solver slow) diff --git a/Examples/Tests/ohm_solver_em_modes/CMakeLists.txt b/Examples/Tests/ohm_solver_em_modes/CMakeLists.txt index e689c83a1e4..a08c321d88d 100644 --- a/Examples/Tests/ohm_solver_em_modes/CMakeLists.txt +++ b/Examples/Tests/ohm_solver_em_modes/CMakeLists.txt @@ -20,3 +20,4 @@ add_warpx_test( diags/diag1000100 # output OFF # dependency ) +label_warpx_test(test_rz_ohm_solver_em_modes_picmi slow) diff --git a/Examples/Tests/ohm_solver_ion_Landau_damping/CMakeLists.txt b/Examples/Tests/ohm_solver_ion_Landau_damping/CMakeLists.txt index 3b2d0bb794b..501b1ce2ced 100644 --- a/Examples/Tests/ohm_solver_ion_Landau_damping/CMakeLists.txt +++ b/Examples/Tests/ohm_solver_ion_Landau_damping/CMakeLists.txt @@ -10,3 +10,4 @@ add_warpx_test( diags/diag1000100 # output OFF # dependency ) +label_warpx_test(test_2d_ohm_solver_landau_damping_picmi slow) diff --git a/Examples/Tests/ohm_solver_ion_beam_instability/CMakeLists.txt b/Examples/Tests/ohm_solver_ion_beam_instability/CMakeLists.txt index 53a9bbdeada..81c6b0d41fd 100644 --- a/Examples/Tests/ohm_solver_ion_beam_instability/CMakeLists.txt +++ b/Examples/Tests/ohm_solver_ion_beam_instability/CMakeLists.txt @@ -10,3 +10,4 @@ add_warpx_test( diags/diag1002500 # output OFF # dependency ) +label_warpx_test(test_1d_ohm_solver_ion_beam_picmi slow) diff --git a/Examples/Tests/particles_in_pml/CMakeLists.txt b/Examples/Tests/particles_in_pml/CMakeLists.txt index c1782dc4d1f..fb539461ec2 100644 --- a/Examples/Tests/particles_in_pml/CMakeLists.txt +++ b/Examples/Tests/particles_in_pml/CMakeLists.txt @@ -40,3 +40,4 @@ add_warpx_test( diags/diag1000200 # output OFF # dependency ) +label_warpx_test(test_3d_particles_in_pml_mr slow) diff --git a/Regression/Checksum/benchmarks_json/test_1d_laser_acceleration_fluid.json b/Regression/Checksum/benchmarks_json/test_1d_laser_acceleration_fluid.json deleted file mode 100644 index 2843e49ce22..00000000000 --- a/Regression/Checksum/benchmarks_json/test_1d_laser_acceleration_fluid.json +++ /dev/null @@ -1,15 +0,0 @@ -{ - "lev=0": { - "Bx": 14264658.38987597, - "By": 0.0, - "Bz": 0.0, - "Ex": 0.0, - "Ey": 4276056420659746.5, - "Ez": 762168740318568.1, - "jx": 0.0, - "jy": 7.47674123799233e+16, - "jz": 4.817762115932484e+17, - "rho": 1609691680.1267354 - } -} - diff --git a/Regression/Checksum/benchmarks_json/test_3d_beam_beam_collision.json b/Regression/Checksum/benchmarks_json/test_3d_beam_beam_collision.json index 799692135ce..ad478e96e2f 100644 --- a/Regression/Checksum/benchmarks_json/test_3d_beam_beam_collision.json +++ b/Regression/Checksum/benchmarks_json/test_3d_beam_beam_collision.json @@ -1,96 +1,96 @@ { "lev=0": { - "Bx": 958613893612.4355, - "By": 958606286084.2804, - "Bz": 50291310.73170665, - "Ex": 2.873993867215236e+20, - "Ey": 2.8740263458334176e+20, - "Ez": 1.3469564521662371e+17, - "rho_beam1": 7.96926761425663e+16, - "rho_beam2": 7.969234189119718e+16, - "rho_ele1": 325562788515568.7, - "rho_ele2": 301373974746269.7, - "rho_pos1": 314013278169396.75, - "rho_pos2": 311298336003435.56 + "Bx": 461189750202.8409, + "By": 461177227889.41614, + "Bz": 20074725.954957746, + "Ex": 1.3827102525768512e+20, + "Ey": 1.3827414851469086e+20, + "Ez": 3.5091001092648376e+16, + "rho_beam1": 3.829927705533972e+16, + "rho_beam2": 3.830015035966193e+16, + "rho_ele1": 162411820580232.94, + "rho_ele2": 150511095307150.62, + "rho_pos1": 159308102449286.12, + "rho_pos2": 156118194376152.7 }, "beam1": { - "particle_opticalDepthQSR": 104848.69019384333, - "particle_position_x": 0.0015001641452988207, - "particle_position_y": 0.0015001296473841738, - "particle_position_z": 0.004965480036291212, - "particle_momentum_x": 6.203657034942546e-15, - "particle_momentum_y": 6.161790111190829e-15, - "particle_momentum_z": 6.806194292286189e-12, - "particle_weight": 635868088.3867786 + "particle_opticalDepthQSR": 52563.18675345914, + "particle_position_x": 0.0007501135479816878, + "particle_position_y": 0.0007501095180907413, + "particle_position_z": 0.002483074586668512, + "particle_momentum_x": 3.0935661449588394e-15, + "particle_momentum_y": 3.070048977445414e-15, + "particle_momentum_z": 3.4017266957145416e-12, + "particle_weight": 636083515.3729652 }, "beam2": { - "particle_opticalDepthQSR": 104197.28107371755, - "particle_position_x": 0.0015001452558405398, - "particle_position_y": 0.0015001281739351966, - "particle_position_z": 0.0049656445643994716, - "particle_momentum_x": 6.202758467582172e-15, - "particle_momentum_y": 6.18910011814166e-15, - "particle_momentum_z": 6.7994521022372906e-12, - "particle_weight": 635874794.3085052 + "particle_opticalDepthQSR": 52275.42552501091, + "particle_position_x": 0.0007500428425956199, + "particle_position_y": 0.0007500867178448842, + "particle_position_z": 0.0024830812114940977, + "particle_momentum_x": 3.1124995623090863e-15, + "particle_momentum_y": 3.094827769550167e-15, + "particle_momentum_z": 3.4015150389915676e-12, + "particle_weight": 636114264.930704 }, "ele1": { - "particle_opticalDepthQSR": 398.7154177999122, - "particle_position_x": 5.2166787076833645e-06, - "particle_position_y": 5.005755590473258e-06, - "particle_position_z": 1.856829463647771e-05, - "particle_momentum_x": 6.0736878569270085e-18, - "particle_momentum_y": 5.735020185191748e-18, - "particle_momentum_z": 2.827581034346608e-15, - "particle_weight": 2602683.4209351614 + "particle_opticalDepthQSR": 156.022199846285, + "particle_position_x": 2.4401923319868757e-06, + "particle_position_y": 2.399150249907213e-06, + "particle_position_z": 8.791577071017722e-06, + "particle_momentum_x": 2.7299291171683895e-18, + "particle_momentum_y": 2.6551510668418435e-18, + "particle_momentum_z": 1.33445643731407e-15, + "particle_weight": 2656838.9769653436 }, "ele2": { - "particle_opticalDepthQSR": 328.6975869797729, - "particle_position_x": 4.984003903707884e-06, - "particle_position_y": 4.695016970410262e-06, - "particle_position_z": 1.606918799511055e-05, - "particle_momentum_x": 4.524294388810778e-18, - "particle_momentum_y": 4.193609622515901e-18, - "particle_momentum_z": 2.624217472737641e-15, - "particle_weight": 2432495.8168380223 + "particle_opticalDepthQSR": 163.79686010701988, + "particle_position_x": 2.724737203764692e-06, + "particle_position_y": 2.9829403746737846e-06, + "particle_position_z": 9.127382649103148e-06, + "particle_momentum_x": 2.1813342297510976e-18, + "particle_momentum_y": 2.7643067192718357e-18, + "particle_momentum_z": 1.259574872512064e-15, + "particle_weight": 2517356.358594387 }, "pho1": { - "particle_opticalDepthBW": 10028.214317531058, - "particle_position_x": 0.00014200324200040716, - "particle_position_y": 0.00014310262095706036, - "particle_position_z": 0.00047470309948487784, + "particle_opticalDepthBW": 5007.597539698783, + "particle_position_x": 7.214053121897416e-05, + "particle_position_y": 7.223804186317301e-05, + "particle_position_z": 0.00024115699590772453, "particle_momentum_x": 0.0, "particle_momentum_y": 0.0, "particle_momentum_z": 0.0, - "particle_weight": 61455533.15171491 + "particle_weight": 62860453.544321 }, "pho2": { - "particle_opticalDepthBW": 10261.48950301913, - "particle_position_x": 0.0001465092909391631, - "particle_position_y": 0.00014555115652303745, - "particle_position_z": 0.00048686081947093, + "particle_opticalDepthBW": 5113.753887045111, + "particle_position_x": 7.271625175781002e-05, + "particle_position_y": 7.311023374122331e-05, + "particle_position_z": 0.00024123464128276151, "particle_momentum_x": 0.0, "particle_momentum_y": 0.0, "particle_momentum_z": 0.0, - "particle_weight": 61924991.09906147 + "particle_weight": 59821371.52413007 }, "pos1": { - "particle_opticalDepthQSR": 380.4787933889546, - "particle_position_x": 5.59226140958729e-06, - "particle_position_y": 5.199149983019462e-06, - "particle_position_z": 1.7261766049926983e-05, - "particle_momentum_x": 5.182944941041321e-18, - "particle_momentum_y": 4.665394338329992e-18, - "particle_momentum_z": 2.565450485567441e-15, - "particle_weight": 2523696.1743423166 + "particle_opticalDepthQSR": 176.87865344045926, + "particle_position_x": 2.597218595285766e-06, + "particle_position_y": 2.5071002403671178e-06, + "particle_position_z": 8.190923176799435e-06, + "particle_momentum_x": 2.409544640420923e-18, + "particle_momentum_y": 2.5096320511498773e-18, + "particle_momentum_z": 1.3349884612525734e-15, + "particle_weight": 2604339.6419650833 }, "pos2": { - "particle_opticalDepthQSR": 378.7526306435402, - "particle_position_x": 4.812490588954386e-06, - "particle_position_y": 4.351750384371962e-06, - "particle_position_z": 1.7621416174292307e-05, - "particle_momentum_x": 4.979887438720444e-18, - "particle_momentum_y": 4.8215630209506066e-18, - "particle_momentum_z": 2.193964301475807e-15, - "particle_weight": 2513162.277112609 + "particle_opticalDepthQSR": 229.50925371797547, + "particle_position_x": 2.6205324097963396e-06, + "particle_position_y": 2.8134541282576216e-06, + "particle_position_z": 9.865542956073817e-06, + "particle_momentum_x": 2.536744632018102e-18, + "particle_momentum_y": 3.035517414633681e-18, + "particle_momentum_z": 1.3203905663185877e-15, + "particle_weight": 2570091.732557297 } } \ No newline at end of file diff --git a/Regression/Checksum/benchmarks_json/test_rz_multiJ_psatd.json b/Regression/Checksum/benchmarks_json/test_rz_multiJ_psatd.json index f30e773e7e0..d9ca66cf0c3 100644 --- a/Regression/Checksum/benchmarks_json/test_rz_multiJ_psatd.json +++ b/Regression/Checksum/benchmarks_json/test_rz_multiJ_psatd.json @@ -1,40 +1,40 @@ { "lev=0": { - "Bt": 24080.416715463354, - "Er": 4536303784778.672, - "Ez": 4298815071343.07, - "jr": 361182004529074.8, - "jz": 1802215706551553.5, - "rho": 4884623.957368025, - "rho_driver": 6288266.101815153, - "rho_plasma_e": 49568366.40537152, - "rho_plasma_p": 50769182.21072973 + "Bt": 5436.145827903481, + "Er": 1084033329144.5951, + "Ez": 1031727477244.4946, + "jr": 86689049352046.02, + "jz": 384358875752423.8, + "rho": 863406.9511533659, + "rho_driver": 1142390.635353391, + "rho_plasma_e": 12207428.74629265, + "rho_plasma_p": 12452342.269855343 }, "driver": { - "particle_momentum_x": 8.723405122353729e-16, - "particle_momentum_y": 8.719285567351437e-16, - "particle_momentum_z": 5.461771334692466e-13, - "particle_position_x": 6.269566322411488, - "particle_position_y": 17.934200805964075, - "particle_theta": 1570790.0436095877, - "particle_weight": 6241484108.424456 - }, - "plasma_p": { - "particle_momentum_x": 6.650539058831456e-20, - "particle_momentum_y": 6.776260514923786e-20, - "particle_momentum_z": 5.470216831432835e-20, - "particle_position_x": 1.1365201443471713, - "particle_position_y": 0.6152066517828133, - "particle_theta": 20286.92798337582, - "particle_weight": 1002457942911.3788 + "particle_momentum_x": 8.723365602723476e-16, + "particle_momentum_y": 8.719184082046568e-16, + "particle_momentum_z": 5.461727890375063e-13, + "particle_position_x": 6.269474429349755, + "particle_position_y": 17.933429487411857, + "particle_theta": 1570777.477238974, + "particle_weight": 6241434176.35186 }, "plasma_e": { - "particle_momentum_x": 6.655027717314839e-20, - "particle_momentum_y": 6.730480164464723e-20, - "particle_momentum_z": 2.8073811669581434e-20, - "particle_position_x": 1.1423427658689635, - "particle_position_y": 0.6140113094028803, - "particle_theta": 20188.939948727297, - "particle_weight": 1002457942911.3788 + "particle_momentum_x": 1.369719083664514e-20, + "particle_momentum_y": 1.5823095211640957e-20, + "particle_momentum_z": 7.189697105606772e-21, + "particle_position_x": 0.28510775565436686, + "particle_position_y": 0.1491507116912657, + "particle_theta": 5011.181913404598, + "particle_weight": 1001422925327.429 + }, + "plasma_p": { + "particle_momentum_x": 1.4836311665634252e-20, + "particle_momentum_y": 1.3653689459385548e-20, + "particle_momentum_z": 1.2505361981099512e-20, + "particle_position_x": 0.2838368160801851, + "particle_position_y": 0.14950442866368444, + "particle_theta": 4995.577541103406, + "particle_weight": 1001422925327.4291 } } \ No newline at end of file diff --git a/Source/BoundaryConditions/PML.H b/Source/BoundaryConditions/PML.H index f6576500c3e..9e7dbc0034c 100644 --- a/Source/BoundaryConditions/PML.H +++ b/Source/BoundaryConditions/PML.H @@ -145,7 +145,8 @@ public: amrex::Real dt, int nox_fft, int noy_fft, int noz_fft, ablastr::utils::enums::GridType grid_type, int do_moving_window, int pml_has_particles, int do_pml_in_domain, - int psatd_solution_type, int J_in_time, int rho_in_time, + PSATDSolutionType psatd_solution_type, + JInTime J_in_time, RhoInTime rho_in_time, bool do_pml_dive_cleaning, bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 1aed10ce654..655d3b59b07 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -552,7 +552,8 @@ PML::PML (const int lev, const BoxArray& grid_ba, Real dt, int nox_fft, int noy_fft, int noz_fft, ablastr::utils::enums::GridType grid_type, int do_moving_window, int /*pml_has_particles*/, int do_pml_in_domain, - const int psatd_solution_type, const int J_in_time, const int rho_in_time, + const PSATDSolutionType psatd_solution_type, + const JInTime J_in_time, const RhoInTime rho_in_time, const bool do_pml_dive_cleaning, const bool do_pml_divb_cleaning, const amrex::IntVect& fill_guards_fields, const amrex::IntVect& fill_guards_current, diff --git a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp index 9de01d20adc..6e4b1b451e2 100644 --- a/Source/BoundaryConditions/WarpXFieldBoundaries.cpp +++ b/Source/BoundaryConditions/WarpXFieldBoundaries.cpp @@ -25,8 +25,8 @@ namespace /** Returns true if any field boundary is set to FieldBoundaryType FT, else returns false.*/ template [[nodiscard]] - bool isAnyBoundary (const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi) + bool isAnyBoundary (const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi) { const auto isFT = [](const auto& b){ return b == FT;}; @@ -37,8 +37,8 @@ namespace /** Returns true if any particle boundary is set to ParticleBoundaryType PT, else returns false.*/ template [[nodiscard]] - bool isAnyBoundary (const amrex::Vector& particle_boundary_lo, - const amrex::Vector& particle_boundary_hi) + bool isAnyBoundary (const amrex::Array& particle_boundary_lo, + const amrex::Array& particle_boundary_hi) { const auto isPT = [](const auto& b){ return b == PT;}; diff --git a/Source/BoundaryConditions/WarpX_PEC.H b/Source/BoundaryConditions/WarpX_PEC.H index a6af894beb4..c387d8c1793 100644 --- a/Source/BoundaryConditions/WarpX_PEC.H +++ b/Source/BoundaryConditions/WarpX_PEC.H @@ -31,8 +31,8 @@ namespace PEC { */ void ApplyPECtoEfield ( std::array Efield, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, const amrex::IntVect& ng_fieldgather, const amrex::Geometry& geom, int lev, PatchType patch_type, const amrex::Vector& ref_ratios, bool split_pml_field = false); @@ -52,8 +52,8 @@ namespace PEC { */ void ApplyPECtoBfield ( std::array Bfield, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, const amrex::IntVect& ng_fieldgather, const amrex::Geometry& geom, int lev, PatchType patch_type, const amrex::Vector& ref_ratios); @@ -73,10 +73,10 @@ namespace PEC { */ void ApplyReflectiveBoundarytoRhofield( amrex::MultiFab* rho, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, - const amrex::Vector& particle_boundary_lo, - const amrex::Vector& particle_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, + const amrex::Array& particle_boundary_lo, + const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, int lev, PatchType patch_type, const amrex::Vector& ref_ratios); @@ -97,10 +97,10 @@ namespace PEC { void ApplyReflectiveBoundarytoJfield( amrex::MultiFab* Jx, amrex::MultiFab* Jy, amrex::MultiFab* Jz, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, - const amrex::Vector& particle_boundary_lo, - const amrex::Vector& particle_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, + const amrex::Array& particle_boundary_lo, + const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, int lev, PatchType patch_type, const amrex::Vector& ref_ratios); @@ -117,8 +117,8 @@ namespace PEC { */ void ApplyPECtoElectronPressure ( amrex::MultiFab* Pefield, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, const amrex::Geometry& geom, int lev, PatchType patch_type, const amrex::Vector& ref_ratios); } diff --git a/Source/BoundaryConditions/WarpX_PEC.cpp b/Source/BoundaryConditions/WarpX_PEC.cpp index 3ad0ab4663e..bedc5b264b7 100644 --- a/Source/BoundaryConditions/WarpX_PEC.cpp +++ b/Source/BoundaryConditions/WarpX_PEC.cpp @@ -457,8 +457,8 @@ namespace void PEC::ApplyPECtoEfield ( std::array Efield, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, const amrex::IntVect& ng_fieldgather, const amrex::Geometry& geom, const int lev, PatchType patch_type, const amrex::Vector& ref_ratios, const bool split_pml_field) @@ -540,8 +540,8 @@ PEC::ApplyPECtoEfield ( void PEC::ApplyPECtoBfield ( std::array Bfield, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, const amrex::IntVect& ng_fieldgather, const amrex::Geometry& geom, const int lev, PatchType patch_type, const amrex::Vector& ref_ratios) { @@ -627,10 +627,10 @@ PEC::ApplyPECtoBfield ( void PEC::ApplyReflectiveBoundarytoRhofield ( amrex::MultiFab* rho, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, - const amrex::Vector& particle_boundary_lo, - const amrex::Vector& particle_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, + const amrex::Array& particle_boundary_lo, + const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, const int lev, PatchType patch_type, const amrex::Vector& ref_ratios) { @@ -713,10 +713,10 @@ PEC::ApplyReflectiveBoundarytoRhofield ( void PEC::ApplyReflectiveBoundarytoJfield( amrex::MultiFab* Jx, amrex::MultiFab* Jy, amrex::MultiFab* Jz, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, - const amrex::Vector& particle_boundary_lo, - const amrex::Vector& particle_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, + const amrex::Array& particle_boundary_lo, + const amrex::Array& particle_boundary_hi, const amrex::Geometry& geom, const int lev, PatchType patch_type, const amrex::Vector& ref_ratios) { @@ -902,8 +902,8 @@ PEC::ApplyReflectiveBoundarytoJfield( void PEC::ApplyPECtoElectronPressure ( amrex::MultiFab* Pefield, - const amrex::Vector& field_boundary_lo, - const amrex::Vector& field_boundary_hi, + const amrex::Array& field_boundary_lo, + const amrex::Array& field_boundary_hi, const amrex::Geometry& geom, const int lev, PatchType patch_type, const amrex::Vector& ref_ratios) { diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index a06b5c2232f..0ddab94d95a 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -241,7 +241,7 @@ MultiParticleContainer::WriteHeader (std::ostream& os) const void storePhiOnParticles ( PinnedMemoryParticleContainer& tmp, - int electrostatic_solver_id, bool is_full_diagnostic ) { + ElectrostaticSolverAlgo electrostatic_solver_id, bool is_full_diagnostic ) { using PinnedParIter = typename PinnedMemoryParticleContainer::ParIterType; diff --git a/Source/Diagnostics/ReducedDiags/FieldReduction.H b/Source/Diagnostics/ReducedDiags/FieldReduction.H index b197ea6b767..bfc2f9aef9e 100644 --- a/Source/Diagnostics/ReducedDiags/FieldReduction.H +++ b/Source/Diagnostics/ReducedDiags/FieldReduction.H @@ -73,7 +73,7 @@ private: std::unique_ptr m_parser; // Type of reduction (e.g. Maximum, Minimum or Sum) - int m_reduction_type; + ReductionType m_reduction_type; public: diff --git a/Source/Diagnostics/ReducedDiags/FieldReduction.cpp b/Source/Diagnostics/ReducedDiags/FieldReduction.cpp index 683ae1921d6..9d0521eb181 100644 --- a/Source/Diagnostics/ReducedDiags/FieldReduction.cpp +++ b/Source/Diagnostics/ReducedDiags/FieldReduction.cpp @@ -60,9 +60,7 @@ FieldReduction::FieldReduction (const std::string& rd_name) parser_string = std::regex_replace(parser_string, std::regex("\n\\s*"), " "); // read reduction type - std::string reduction_type_string; - pp_rd_name.get("reduction_type", reduction_type_string); - m_reduction_type = GetAlgorithmInteger (pp_rd_name, "reduction_type"); + pp_rd_name.get_enum_sloppy("reduction_type", m_reduction_type, "-_"); if (amrex::ParallelDescriptor::IOProcessor()) { @@ -77,7 +75,7 @@ FieldReduction::FieldReduction (const std::string& rd_name) ofs << m_sep; ofs << "[" << c++ << "]time(s)"; ofs << m_sep; - ofs << "[" << c++ << "]" + reduction_type_string + " of " + parser_string + " (SI units)"; + ofs << "[" << c++ << "]" + amrex::getEnumNameString(m_reduction_type) + " of " + parser_string + " (SI units)"; ofs << std::endl; // close file diff --git a/Source/FieldSolver/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolver.cpp index 07dbc8c1424..69added4524 100644 --- a/Source/FieldSolver/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolver.cpp @@ -1084,7 +1084,7 @@ void ElectrostaticSolver::PoissonBoundaryHandler::definePhiBCs (const amrex::Geo else { WARPX_ABORT_WITH_MESSAGE( "Field boundary conditions have to be either periodic, PEC or neumann " - "when using the electrostatic Multigrid solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim]) + "when using the electrostatic Multigrid solver, but they are " + amrex::getEnumNameString(WarpX::field_boundary_lo[idim]) ); } @@ -1099,7 +1099,7 @@ void ElectrostaticSolver::PoissonBoundaryHandler::definePhiBCs (const amrex::Geo else { WARPX_ABORT_WITH_MESSAGE( "Field boundary conditions have to be either periodic, PEC or neumann " - "when using the electrostatic Multigrid solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_hi[idim]) + "when using the electrostatic Multigrid solver, but they are " + amrex::getEnumNameString(WarpX::field_boundary_hi[idim]) ); } } diff --git a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp index 9ae7c6fb6fd..e72260fcf4f 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/ApplySilverMuellerBoundary.cpp @@ -39,8 +39,8 @@ void FiniteDifferenceSolver::ApplySilverMuellerBoundary ( ablastr::fields::VectorField& Bfield, amrex::Box domain_box, amrex::Real const dt, - amrex::Vector field_boundary_lo, - amrex::Vector field_boundary_hi) { + amrex::Array field_boundary_lo, + amrex::Array field_boundary_hi) { using ablastr::fields::Direction; diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H index e3a113136a2..250f6280f92 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H @@ -48,7 +48,7 @@ class FiniteDifferenceSolver * \param grid_type Whether the solver is applied to a collocated or staggered grid */ FiniteDifferenceSolver ( - int fdtd_algo, + ElectromagneticSolverAlgo fdtd_algo, std::array cell_size, ablastr::utils::enums::GridType grid_type ); @@ -86,8 +86,8 @@ class FiniteDifferenceSolver ablastr::fields::VectorField & Bfield, amrex::Box domain_box, amrex::Real dt, - amrex::Vector field_boundary_lo, - amrex::Vector field_boundary_hi); + amrex::Array field_boundary_lo, + amrex::Array field_boundary_hi); void ComputeDivE ( ablastr::fields::VectorField const & Efield, @@ -180,7 +180,7 @@ class FiniteDifferenceSolver private: - int m_fdtd_algo; + ElectromagneticSolverAlgo m_fdtd_algo; ablastr::utils::enums::GridType m_grid_type; #ifdef WARPX_DIM_RZ diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp index 9af610031c0..fdd02c5249b 100644 --- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp +++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp @@ -28,7 +28,7 @@ /* This function initializes the stencil coefficients for the chosen finite-difference algorithm */ FiniteDifferenceSolver::FiniteDifferenceSolver ( - int const fdtd_algo, + ElectromagneticSolverAlgo const fdtd_algo, std::array cell_size, ablastr::utils::enums::GridType grid_type): // Register the type of finite-difference algorithm diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp index 3d5fdad264e..0e329bc4730 100644 --- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp +++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp @@ -339,7 +339,7 @@ void MagnetostaticSolver::VectorPoissonBoundaryHandler::defineVectorPotentialBCs else { WARPX_ABORT_WITH_MESSAGE( "Field boundary conditions have to be either periodic, PEC or neumann " - "when using the magnetostatic solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim]) + "when using the magnetostatic solver, but they are " + amrex::getEnumNameString(WarpX::field_boundary_lo[idim]) ); } @@ -359,7 +359,7 @@ void MagnetostaticSolver::VectorPoissonBoundaryHandler::defineVectorPotentialBCs else { WARPX_ABORT_WITH_MESSAGE( "Field boundary conditions have to be either periodic, PEC or neumann " - "when using the magnetostatic solver, but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim]) + "when using the magnetostatic solver, but they are " + amrex::getEnumNameString(WarpX::field_boundary_lo[idim]) ); } } diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H index 23c38f33a85..94a16c13ec1 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.H @@ -55,8 +55,8 @@ class PsatdAlgorithmFirstOrder : public SpectralBaseAlgorithm ablastr::utils::enums::GridType grid_type, amrex::Real dt, bool div_cleaning, - int J_in_time, - int rho_in_time); + JInTime J_in_time, + RhoInTime rho_in_time); /** * \brief Updates E, B, F, and G fields in spectral space, @@ -93,8 +93,8 @@ class PsatdAlgorithmFirstOrder : public SpectralBaseAlgorithm // Other member variables amrex::Real m_dt; bool m_div_cleaning; - int m_J_in_time; - int m_rho_in_time; + JInTime m_J_in_time; + RhoInTime m_rho_in_time; }; #endif // WARPX_USE_FFT #endif // WARPX_PSATD_ALGORITHM_FIRST_ORDER_H_ diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp index 14db3c9af48..c5f60e18d24 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmFirstOrder.cpp @@ -38,8 +38,8 @@ PsatdAlgorithmFirstOrder::PsatdAlgorithmFirstOrder ( ablastr::utils::enums::GridType grid_type, const amrex::Real dt, const bool div_cleaning, - const int J_in_time, - const int rho_in_time + const JInTime J_in_time, + const RhoInTime rho_in_time ) // Initializer list : SpectralBaseAlgorithm(spectral_kspace, dm, spectral_index, norder_x, norder_y, norder_z, grid_type), diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H index c2efbd79935..2cbb4d7402d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H @@ -27,8 +27,8 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ amrex::Real dt_step, bool update_with_rho, bool time_averaging, - int J_in_time, - int rho_in_time, + JInTime J_in_time, + RhoInTime rho_in_time, bool dive_cleaning, bool divb_cleaning); // Redefine functions from base class @@ -65,7 +65,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ amrex::Real m_dt; bool m_update_with_rho; bool m_time_averaging; - int m_J_in_time; + JInTime m_J_in_time; bool m_dive_cleaning; bool m_divb_cleaning; SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef; diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp index efff5c30a41..66d99e29715 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp @@ -25,8 +25,8 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace, amrex::Real const dt, bool const update_with_rho, const bool time_averaging, - const int J_in_time, - const int rho_in_time, + const JInTime J_in_time, + const RhoInTime rho_in_time, const bool dive_cleaning, const bool divb_cleaning): // Initialize members of base class and member variables diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H index d6c4916bdac..42bf32c4724 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H @@ -11,6 +11,7 @@ #include "SpectralFieldData_fwd.H" #include "SpectralKSpace.H" +#include "Utils/WarpXAlgorithmSelection.H" #include "Utils/WarpX_Complex.H" #include @@ -59,8 +60,8 @@ class SpectralFieldIndex */ SpectralFieldIndex (bool update_with_rho, bool time_averaging, - int J_in_time, - int rho_in_time, + JInTime J_in_time, + RhoInTime rho_in_time, bool dive_cleaning, bool divb_cleaning, bool pml, diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 8e7b9ed9ae4..a20429eea68 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -35,8 +35,8 @@ using namespace amrex; SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho, const bool time_averaging, - const int J_in_time, - const int rho_in_time, + const JInTime J_in_time, + const RhoInTime rho_in_time, const bool dive_cleaning, const bool divb_cleaning, const bool pml, diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H index 145ea317f7b..bcd80e421a8 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H @@ -11,6 +11,7 @@ #include "SpectralAlgorithms/SpectralBaseAlgorithm.H" #include "SpectralFieldData.H" +#include "Utils/WarpXAlgorithmSelection.H" #include #include @@ -87,9 +88,9 @@ class SpectralSolver bool periodic_single_box, bool update_with_rho, bool fft_do_time_averaging, - int psatd_solution_type, - int J_in_time, - int rho_in_time, + PSATDSolutionType psatd_solution_type, + JInTime J_in_time, + RhoInTime rho_in_time, bool dive_cleaning, bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp index c6192187064..fcd52597e07 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp @@ -37,9 +37,9 @@ SpectralSolver::SpectralSolver ( const bool pml, const bool periodic_single_box, const bool update_with_rho, const bool fft_do_time_averaging, - const int psatd_solution_type, - const int J_in_time, - const int rho_in_time, + const PSATDSolutionType psatd_solution_type, + const JInTime J_in_time, + const RhoInTime rho_in_time, const bool dive_cleaning, const bool divb_cleaning) { diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H index 040acd00c51..61cf64036eb 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H @@ -42,8 +42,8 @@ class SpectralSolverRZ bool with_pml, bool update_with_rho, bool fft_do_time_averaging, - int J_in_time, - int rho_in_time, + JInTime J_in_time, + RhoInTime rho_in_time, bool dive_cleaning, bool divb_cleaning); diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp index b558f3983a3..9a8cff9f1f3 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp @@ -36,8 +36,8 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev, bool const with_pml, bool const update_with_rho, const bool fft_do_time_averaging, - const int J_in_time, - const int rho_in_time, + const JInTime J_in_time, + const RhoInTime rho_in_time, const bool dive_cleaning, const bool divb_cleaning) : k_space(realspace_ba, dm, dx) diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index a095bce6ae3..66976045943 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -64,34 +64,25 @@ void BilinearFilter::ComputeStencils(){ WARPX_PROFILE("BilinearFilter::ComputeStencils()"); int i = 0; for (const auto& el : npass_each_dir ) { - stencil_length_each_dir[i++] = static_cast(el); + stencil_length_each_dir[i++] = static_cast(el) + 1; } - stencil_length_each_dir += 1.; + + m_stencil_0.resize( 1u + npass_each_dir[0] ); + compute_stencil(m_stencil_0, npass_each_dir[0]); +#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) || defined(WARPX_DIM_3D) + m_stencil_1.resize( 1u + npass_each_dir[1] ); + compute_stencil(m_stencil_1, npass_each_dir[1]); +#endif #if defined(WARPX_DIM_3D) - // npass_each_dir = npass_x npass_y npass_z - stencil_x.resize( 1u + npass_each_dir[0] ); - stencil_y.resize( 1u + npass_each_dir[1] ); - stencil_z.resize( 1u + npass_each_dir[2] ); - compute_stencil(stencil_x, npass_each_dir[0]); - compute_stencil(stencil_y, npass_each_dir[1]); - compute_stencil(stencil_z, npass_each_dir[2]); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - // npass_each_dir = npass_x npass_z - stencil_x.resize( 1u + npass_each_dir[0] ); - stencil_z.resize( 1u + npass_each_dir[1] ); - compute_stencil(stencil_x, npass_each_dir[0]); - compute_stencil(stencil_z, npass_each_dir[1]); -#elif defined(WARPX_DIM_1D_Z) - // npass_each_dir = npass_z - stencil_z.resize( 1u + npass_each_dir[0] ); - compute_stencil(stencil_z, npass_each_dir[0]); + m_stencil_2.resize( 1u + npass_each_dir[2] ); + compute_stencil(m_stencil_2, npass_each_dir[2]); #endif + slen = stencil_length_each_dir.dim3(); -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) +#if defined(WARPX_DIM_1D_Z) || defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) slen.z = 1; #endif #if defined(WARPX_DIM_1D_Z) slen.y = 1; - slen.z = 1; #endif } diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index 5c648a75c4c..584f6b151d7 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -21,9 +21,9 @@ public: // Apply stencil on MultiFab. // Guard cells are handled inside this function - void ApplyStencil(amrex::MultiFab& dstmf, - const amrex::MultiFab& srcmf, int lev, int scomp=0, - int dcomp=0, int ncomp=10000); + void ApplyStencil (amrex::MultiFab& dstmf, + const amrex::MultiFab& srcmf, int lev, int scomp=0, + int dcomp=0, int ncomp=10000); // Apply stencil on a FabArray. void ApplyStencil (amrex::FArrayBox& dstfab, @@ -31,20 +31,18 @@ public: int scomp=0, int dcomp=0, int ncomp=10000); // public for cuda - void DoFilter(const amrex::Box& tbx, - amrex::Array4 const& tmp, - amrex::Array4 const& dst, - int scomp, int dcomp, int ncomp); + void DoFilter (const amrex::Box& tbx, + amrex::Array4 const& tmp, + amrex::Array4 const& dst, + int scomp, int dcomp, int ncomp); - // In 2D, stencil_length_each_dir = {length(stencil_x), length(stencil_z)} + // Length of stencil in each included direction amrex::IntVect stencil_length_each_dir; protected: // Stencil along each direction. - // in 2D, stencil_y is not initialized. - amrex::Gpu::DeviceVector stencil_x, stencil_y, stencil_z; - // Length of each stencil. - // In 2D, slen = {length(stencil_x), length(stencil_z), 1} + amrex::Gpu::DeviceVector m_stencil_0, m_stencil_1, m_stencil_2; + // Length of each stencil, 1 for dimensions not included amrex::Dim3 slen; private: diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 6243ce4ebbf..40dfc54f8cb 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -87,23 +87,21 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, DoFilter(tbx, src, dst, scomp, dcomp, ncomp); } -/* \brief Apply stencil (2D/3D, CPU/GPU) +/* \brief Apply stencil (CPU/GPU) */ void Filter::DoFilter (const Box& tbx, Array4 const& src, Array4 const& dst, int scomp, int dcomp, int ncomp) { -#if (AMREX_SPACEDIM >= 2) - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); -#endif -#if defined(WARPX_DIM_3D) - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); -#endif - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + AMREX_D_TERM( + amrex::Real const* AMREX_RESTRICT s0 = m_stencil_0.data();, + amrex::Real const* AMREX_RESTRICT s1 = m_stencil_1.data();, + amrex::Real const* AMREX_RESTRICT s2 = m_stencil_2.data(); + ) Dim3 slen_local = slen; -#if defined(WARPX_DIM_3D) +#if AMREX_SPACEDIM == 3 AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, { Real d = 0.0; @@ -115,25 +113,25 @@ void Filter::DoFilter (const Box& tbx, return src.contains(jj,kk,ll) ? src(jj,kk,ll,nn) : 0.0_rt; }; - for (int iz=0; iz < slen_local.z; ++iz){ - for (int iy=0; iy < slen_local.y; ++iy){ - for (int ix=0; ix < slen_local.x; ++ix){ - Real sss = sx[ix]*sy[iy]*sz[iz]; - d += sss*( src_zeropad(i-ix,j-iy,k-iz,scomp+n) - +src_zeropad(i+ix,j-iy,k-iz,scomp+n) - +src_zeropad(i-ix,j+iy,k-iz,scomp+n) - +src_zeropad(i+ix,j+iy,k-iz,scomp+n) - +src_zeropad(i-ix,j-iy,k+iz,scomp+n) - +src_zeropad(i+ix,j-iy,k+iz,scomp+n) - +src_zeropad(i-ix,j+iy,k+iz,scomp+n) - +src_zeropad(i+ix,j+iy,k+iz,scomp+n)); + for (int i2=0; i2 < slen_local.z; ++i2){ + for (int i1=0; i1 < slen_local.y; ++i1){ + for (int i0=0; i0 < slen_local.x; ++i0){ + Real sss = s0[i0]*s1[i1]*s2[i2]; + d += sss*( src_zeropad(i-i0,j-i1,k-i2,scomp+n) + +src_zeropad(i+i0,j-i1,k-i2,scomp+n) + +src_zeropad(i-i0,j+i1,k-i2,scomp+n) + +src_zeropad(i+i0,j+i1,k-i2,scomp+n) + +src_zeropad(i-i0,j-i1,k+i2,scomp+n) + +src_zeropad(i+i0,j-i1,k+i2,scomp+n) + +src_zeropad(i-i0,j+i1,k+i2,scomp+n) + +src_zeropad(i+i0,j+i1,k+i2,scomp+n)); } } } dst(i,j,k,dcomp+n) = d; }); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) +#elif AMREX_SPACEDIM == 2 AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, { Real d = 0.0; @@ -145,21 +143,21 @@ void Filter::DoFilter (const Box& tbx, return src.contains(jj,kk,ll) ? src(jj,kk,ll,nn) : 0.0_rt; }; - for (int iz=0; iz < slen_local.z; ++iz){ - for (int iy=0; iy < slen_local.y; ++iy){ - for (int ix=0; ix < slen_local.x; ++ix){ - Real sss = sx[ix]*sz[iy]; - d += sss*( src_zeropad(i-ix,j-iy,k,scomp+n) - +src_zeropad(i+ix,j-iy,k,scomp+n) - +src_zeropad(i-ix,j+iy,k,scomp+n) - +src_zeropad(i+ix,j+iy,k,scomp+n)); + for (int i2=0; i2 < slen_local.z; ++i2){ + for (int i1=0; i1 < slen_local.y; ++i1){ + for (int i0=0; i0 < slen_local.x; ++i0){ + Real sss = s0[i0]*s1[i1]; + d += sss*( src_zeropad(i-i0,j-i1,k,scomp+n) + +src_zeropad(i+i0,j-i1,k,scomp+n) + +src_zeropad(i-i0,j+i1,k,scomp+n) + +src_zeropad(i+i0,j+i1,k,scomp+n)); } } } dst(i,j,k,dcomp+n) = d; }); -#elif defined(WARPX_DIM_1D_Z) +#elif AMREX_SPACEDIM == 1 AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, { Real d = 0.0; @@ -171,21 +169,18 @@ void Filter::DoFilter (const Box& tbx, return src.contains(jj,kk,ll) ? src(jj,kk,ll,nn) : 0.0_rt; }; - for (int iz=0; iz < slen_local.z; ++iz){ - for (int iy=0; iy < slen_local.y; ++iy){ - for (int ix=0; ix < slen_local.x; ++ix){ - Real sss = sz[iy]; - d += sss*( src_zeropad(i-ix,j,k,scomp+n) - +src_zeropad(i+ix,j,k,scomp+n)); + for (int i2=0; i2 < slen_local.z; ++i2){ + for (int i1=0; i1 < slen_local.y; ++i1){ + for (int i0=0; i0 < slen_local.x; ++i0){ + Real sss = s0[i0]; + d += sss*( src_zeropad(i-i0,j,k,scomp+n) + +src_zeropad(i+i0,j,k,scomp+n)); } } } dst(i,j,k,dcomp+n) = d; }); -#else - WARPX_ABORT_WITH_MESSAGE( - "Filter not implemented for the current geometry!"); #endif } @@ -278,13 +273,11 @@ void Filter::DoFilter (const Box& tbx, const auto lo = amrex::lbound(tbx); const auto hi = amrex::ubound(tbx); // tmp and dst are of type Array4 (Fortran ordering) -#if (AMREX_SPACEDIM >= 2) - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); -#endif -#if defined(WARPX_DIM_3D) - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); -#endif - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + AMREX_D_TERM( + amrex::Real const* AMREX_RESTRICT s0 = m_stencil_0.data();, + amrex::Real const* AMREX_RESTRICT s1 = m_stencil_1.data();, + amrex::Real const* AMREX_RESTRICT s2 = m_stencil_2.data(); + ) for (int n = 0; n < ncomp; ++n) { // Set dst value to 0. for (int k = lo.z; k <= hi.z; ++k) { @@ -295,41 +288,32 @@ void Filter::DoFilter (const Box& tbx, } } // 3 nested loop on 3D stencil - for (int iz=0; iz < slen.z; ++iz){ - for (int iy=0; iy < slen.y; ++iy){ - for (int ix=0; ix < slen.x; ++ix){ -#if defined(WARPX_DIM_3D) - const Real sss = sx[ix]*sy[iy]*sz[iz]; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const Real sss = sx[ix]*sz[iy]; -#else - const Real sss = sz[ix]; -#endif + for (int i2=0; i2 < slen.z; ++i2){ + for (int i1=0; i1 < slen.y; ++i1){ + for (int i0=0; i0 < slen.x; ++i0){ + const Real sss = AMREX_D_TERM(s0[i0], *s1[i1], *s2[i2]); // 3 nested loop on 3D array for (int k = lo.z; k <= hi.z; ++k) { for (int j = lo.y; j <= hi.y; ++j) { AMREX_PRAGMA_SIMD for (int i = lo.x; i <= hi.x; ++i) { -#if defined(WARPX_DIM_3D) - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n) - +tmp(i+ix,j-iy,k-iz,scomp+n) - +tmp(i-ix,j+iy,k-iz,scomp+n) - +tmp(i+ix,j+iy,k-iz,scomp+n) - +tmp(i-ix,j-iy,k+iz,scomp+n) - +tmp(i+ix,j-iy,k+iz,scomp+n) - +tmp(i-ix,j+iy,k+iz,scomp+n) - +tmp(i+ix,j+iy,k+iz,scomp+n)); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) - +tmp(i+ix,j-iy,k,scomp+n) - +tmp(i-ix,j+iy,k,scomp+n) - +tmp(i+ix,j+iy,k,scomp+n)); -#elif defined(WARPX_DIM_1D_Z) - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j,k,scomp+n) - +tmp(i+ix,j,k,scomp+n)); -#else - WARPX_ABORT_WITH_MESSAGE( - "Filter not implemented for the current geometry!"); +#if AMREX_SPACEDIM == 3 + dst(i,j,k,dcomp+n) += sss*(tmp(i-i0,j-i1,k-i2,scomp+n) + +tmp(i+i0,j-i1,k-i2,scomp+n) + +tmp(i-i0,j+i1,k-i2,scomp+n) + +tmp(i+i0,j+i1,k-i2,scomp+n) + +tmp(i-i0,j-i1,k+i2,scomp+n) + +tmp(i+i0,j-i1,k+i2,scomp+n) + +tmp(i-i0,j+i1,k+i2,scomp+n) + +tmp(i+i0,j+i1,k+i2,scomp+n)); +#elif AMREX_SPACEDIM == 2 + dst(i,j,k,dcomp+n) += sss*(tmp(i-i0,j-i1,k,scomp+n) + +tmp(i+i0,j-i1,k,scomp+n) + +tmp(i-i0,j+i1,k,scomp+n) + +tmp(i+i0,j+i1,k,scomp+n)); +#elif AMREX_SPACEDIM == 1 + dst(i,j,k,dcomp+n) += sss*(tmp(i-i0,j,k,scomp+n) + +tmp(i+i0,j,k,scomp+n)); #endif } } diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 9567bdf1bb2..a73efb0ec64 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -121,17 +121,18 @@ void NCIGodfreyFilter::ComputeStencils() # endif h_stencil_z[0] /= 2._rt; - stencil_x.resize(h_stencil_x.size()); + m_stencil_0.resize(h_stencil_x.size()); + Gpu::copyAsync(Gpu::hostToDevice,h_stencil_x.begin(),h_stencil_x.end(),m_stencil_0.begin()); # if defined(WARPX_DIM_3D) - stencil_y.resize(h_stencil_y.size()); + m_stencil_1.resize(h_stencil_y.size()); + m_stencil_2.resize(h_stencil_z.size()); + Gpu::copyAsync(Gpu::hostToDevice,h_stencil_y.begin(),h_stencil_y.end(),m_stencil_1.begin()); + Gpu::copyAsync(Gpu::hostToDevice,h_stencil_z.begin(),h_stencil_z.end(),m_stencil_2.begin()); +# elif (AMREX_SPACEDIM == 2) + // In 2D, the filter applies stencil_1 to the 2nd dimension + m_stencil_1.resize(h_stencil_z.size()); + Gpu::copyAsync(Gpu::hostToDevice,h_stencil_z.begin(),h_stencil_z.end(),m_stencil_1.begin()); # endif - stencil_z.resize(h_stencil_z.size()); - - Gpu::copyAsync(Gpu::hostToDevice,h_stencil_x.begin(),h_stencil_x.end(),stencil_x.begin()); -# if defined(WARPX_DIM_3D) - Gpu::copyAsync(Gpu::hostToDevice,h_stencil_y.begin(),h_stencil_y.end(),stencil_y.begin()); -# endif - Gpu::copyAsync(Gpu::hostToDevice,h_stencil_z.begin(),h_stencil_z.end(),stencil_z.begin()); Gpu::synchronize(); } diff --git a/Source/Parallelization/GuardCellManager.H b/Source/Parallelization/GuardCellManager.H index c70bd6d3a35..c9de3d8b250 100644 --- a/Source/Parallelization/GuardCellManager.H +++ b/Source/Parallelization/GuardCellManager.H @@ -7,6 +7,7 @@ #ifndef WARPX_GUARDCELLMANAGER_H_ #define WARPX_GUARDCELLMANAGER_H_ +#include #include #include @@ -63,7 +64,7 @@ public: int nox, int nox_fft, int noy_fft, int noz_fft, int nci_corr_stencil, - int electromagnetic_solver_id, + ElectromagneticSolverAlgo electromagnetic_solver_id, int max_level, const amrex::Vector& v_galilean, const amrex::Vector& v_comoving, diff --git a/Source/Parallelization/GuardCellManager.cpp b/Source/Parallelization/GuardCellManager.cpp index c8c0ff40acf..a36e39d9497 100644 --- a/Source/Parallelization/GuardCellManager.cpp +++ b/Source/Parallelization/GuardCellManager.cpp @@ -42,7 +42,7 @@ guardCellManager::Init ( const int nox, const int nox_fft, const int noy_fft, const int noz_fft, const int nci_corr_stencil, - const int electromagnetic_solver_id, + const ElectromagneticSolverAlgo electromagnetic_solver_id, const int max_level, const amrex::Vector& v_galilean, const amrex::Vector& v_comoving, diff --git a/Source/Parallelization/WarpXComm_K.H b/Source/Parallelization/WarpXComm_K.H index c3362087ad9..79f2b34fba0 100644 --- a/Source/Parallelization/WarpXComm_K.H +++ b/Source/Parallelization/WarpXComm_K.H @@ -44,18 +44,19 @@ void warpx_interp (int j, int k, int l, // Refinement ratio const int rj = rr[0]; - const int rk = (AMREX_SPACEDIM == 1) ? 1 : rr[1]; - const int rl = (AMREX_SPACEDIM <= 2) ? 1 : rr[2]; + const int rk = (AMREX_SPACEDIM > 1) ? rr[1] : 1; + const int rl = (AMREX_SPACEDIM > 2) ? rr[2] : 1; // Staggering (0: cell-centered; 1: nodal) + // Unused dimensions are considered nodal. const int sj = arr_stag[0]; - const int sk = (AMREX_SPACEDIM == 1) ? 0 : arr_stag[1]; - const int sl = (AMREX_SPACEDIM <= 2) ? 0 : arr_stag[2]; + const int sk = (AMREX_SPACEDIM > 1) ? arr_stag[1] : 1; + const int sl = (AMREX_SPACEDIM > 2) ? arr_stag[2] : 1; // Number of points used for interpolation from coarse grid to fine grid const int nj = 2; - const int nk = 2; - const int nl = 2; + const int nk = (AMREX_SPACEDIM > 1) ? 2 : 1; + const int nl = (AMREX_SPACEDIM > 2) ? 2 : 1; const int jc = (sj == 0) ? amrex::coarsen(j - rj/2, rj) : amrex::coarsen(j, rj); const int kc = (sk == 0) ? amrex::coarsen(k - rk/2, rk) : amrex::coarsen(k, rk); @@ -133,38 +134,20 @@ void warpx_interp (int j, int k, int l, // Refinement ratio const int rj = rr[0]; -#if defined(WARPX_DIM_1D_Z) - constexpr int rk = 1; - constexpr int rl = 1; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const int rk = rr[1]; - constexpr int rl = 1; -#else - const int rk = rr[1]; - const int rl = rr[2]; -#endif + const int rk = (AMREX_SPACEDIM > 1) ? rr[1] : 1; + const int rl = (AMREX_SPACEDIM > 2) ? rr[2] : 1; // Staggering of fine array (0: cell-centered; 1: nodal) + // Unused dimensions are considered nodal. const int sj_fp = arr_fine_stag[0]; -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const int sk_fp = arr_fine_stag[1]; -#elif defined(WARPX_DIM_3D) - const int sk_fp = arr_fine_stag[1]; - const int sl_fp = arr_fine_stag[2]; -#endif + const int sk_fp = (AMREX_SPACEDIM > 1) ? arr_fine_stag[1] : 1; + const int sl_fp = (AMREX_SPACEDIM > 2) ? arr_fine_stag[2] : 1; // Staggering of coarse array (0: cell-centered; 1: nodal) + // Unused dimensions are considered nodal. const int sj_cp = arr_coarse_stag[0]; -#if defined(WARPX_DIM_1D_Z) - constexpr int sk_cp = 0; - constexpr int sl_cp = 0; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const int sk_cp = arr_coarse_stag[1]; - constexpr int sl_cp = 0; -#else - const int sk_cp = arr_coarse_stag[1]; - const int sl_cp = arr_coarse_stag[2]; -#endif + const int sk_cp = (AMREX_SPACEDIM > 1) ? arr_coarse_stag[1] : 1; + const int sl_cp = (AMREX_SPACEDIM > 2) ? arr_coarse_stag[2] : 1; // Number of points used for interpolation from coarse grid to fine grid int nj; @@ -182,27 +165,19 @@ void warpx_interp (int j, int k, int l, // 1) Interpolation from coarse nodal to fine nodal nj = 2; -#if defined(WARPX_DIM_1D_Z) - nk = 1; - nl = 1; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - nk = 2; - nl = 1; -#else - nk = 2; - nl = 2; -#endif + nk = (AMREX_SPACEDIM > 1) ? 2 : 1; + nl = (AMREX_SPACEDIM > 2) ? 2 : 1; for (int jj = 0; jj < nj; jj++) { for (int kk = 0; kk < nk; kk++) { for (int ll = 0; ll < nl; ll++) { auto c = arr_tmp_zeropad(jc+jj,kc+kk,lc+ll); c *= (rj - amrex::Math::abs(j - (jc + jj) * rj)) / static_cast(rj); -#if (AMREX_SPACEDIM >= 2) +#if (AMREX_SPACEDIM > 1) c *= (rk - amrex::Math::abs(k - (kc + kk) * rk)) / static_cast(rk); -#endif -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM > 2) c *= (rl - amrex::Math::abs(l - (lc + ll) * rl)) / static_cast(rl); +#endif #endif tmp += c; } @@ -212,16 +187,8 @@ void warpx_interp (int j, int k, int l, // 2) Interpolation from coarse staggered to fine nodal nj = 2; -#if defined(WARPX_DIM_1D_Z) - nk = 1; - nl = 1; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - nk = 2; - nl = 1; -#else - nk = 2; - nl = 2; -#endif + nk = (AMREX_SPACEDIM > 1) ? 2 : 1; + nl = (AMREX_SPACEDIM > 2) ? 2 : 1; const int jn = (sj_cp == 1) ? j : j - rj / 2; const int kn = (sk_cp == 1) ? k : k - rk / 2; @@ -236,11 +203,11 @@ void warpx_interp (int j, int k, int l, for (int ll = 0; ll < nl; ll++) { auto c = arr_coarse_zeropad(jc+jj,kc+kk,lc+ll); c *= (rj - amrex::Math::abs(jn - (jc + jj) * rj)) / static_cast(rj); -#if (AMREX_SPACEDIM >= 2) +#if (AMREX_SPACEDIM > 1) c *= (rk - amrex::Math::abs(kn - (kc + kk) * rk)) / static_cast(rk); -#endif -#if (AMREX_SPACEDIM == 3) +#if (AMREX_SPACEDIM > 2) c *= (rl - amrex::Math::abs(ln - (lc + ll) * rl)) / static_cast(rl); +#endif #endif coarse += c; } @@ -250,28 +217,12 @@ void warpx_interp (int j, int k, int l, // 3) Interpolation from fine staggered to fine nodal nj = (sj_fp == 0) ? 2 : 1; -#if defined(WARPX_DIM_1D_Z) - nk = 1; - nl = 1; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - nk = (sk_fp == 0) ? 2 : 1; - nl = 1; -#else nk = (sk_fp == 0) ? 2 : 1; nl = (sl_fp == 0) ? 2 : 1; -#endif const int jm = (sj_fp == 0) ? j-1 : j; -#if defined(WARPX_DIM_1D_Z) - const int km = k; - const int lm = l; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const int km = (sk_fp == 0) ? k-1 : k; - const int lm = l; -#else const int km = (sk_fp == 0) ? k-1 : k; const int lm = (sl_fp == 0) ? l-1 : l; -#endif for (int jj = 0; jj < nj; jj++) { for (int kk = 0; kk < nk; kk++) { @@ -285,6 +236,7 @@ void warpx_interp (int j, int k, int l, // Final result arr_aux(j,k,l) = tmp + (fine - coarse); } + /** * \brief Interpolation function called within WarpX::UpdateAuxilaryDataStagToNodal * to interpolate data from the coarse and fine grids to the fine aux grid, @@ -320,13 +272,10 @@ void warpx_interp (int j, int k, int l, // - (x,y,z) in 3D // Staggering of fine array (0: cell-centered; 1: nodal) + // Unused dimensions are considered nodal. const int sj_fp = arr_fine_stag[0]; -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const int sk_fp = arr_fine_stag[1]; -#elif defined(WARPX_DIM_3D) - const int sk_fp = arr_fine_stag[1]; - const int sl_fp = arr_fine_stag[2]; -#endif + const int sk_fp = (AMREX_SPACEDIM > 1) ? arr_fine_stag[1] : 1; + const int sl_fp = (AMREX_SPACEDIM > 2) ? arr_fine_stag[2] : 1; // Number of points used for interpolation from coarse grid to fine grid int nj; @@ -338,28 +287,12 @@ void warpx_interp (int j, int k, int l, // 3) Interpolation from fine staggered to fine nodal nj = (sj_fp == 0) ? 2 : 1; -#if defined(WARPX_DIM_1D_Z) - nk = 1; - nl = 1; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - nk = (sk_fp == 0) ? 2 : 1; - nl = 1; -#else nk = (sk_fp == 0) ? 2 : 1; nl = (sl_fp == 0) ? 2 : 1; -#endif - const int jm = (sj_fp == 0) ? j-1 : j; -#if defined(WARPX_DIM_1D_Z) - const int km = k; - const int lm = l; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const int km = (sk_fp == 0) ? k-1 : k; - const int lm = l; -#else - const int km = (sk_fp == 0) ? k-1 : k; - const int lm = (sl_fp == 0) ? l-1 : l; -#endif + int const jm = (sj_fp == 0) ? j-1 : j; + int const km = (sk_fp == 0) ? k-1 : k; + int const lm = (sl_fp == 0) ? l-1 : l; for (int jj = 0; jj < nj; jj++) { for (int kk = 0; kk < nk; kk++) { @@ -418,11 +351,7 @@ void warpx_interp (const int j, }; // Avoid compiler warnings -#if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(nox, noy, stencil_coeffs_x, stencil_coeffs_y); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(noy, stencil_coeffs_y); -#endif // If dst_nodal = true , we are centering from a staggered grid to a nodal grid // If dst_nodal = false, we are centering from a nodal grid to a staggered grid @@ -432,70 +361,32 @@ void warpx_interp (const int j, const int shift = (dst_nodal) ? 0 : 1; // Staggering (s = 0 if cell-centered, s = 1 if nodal) + // Unused dimensions are considered nodal. const int sj = (dst_nodal) ? src_stag[0] : dst_stag[0]; -#if (AMREX_SPACEDIM >= 2) - const int sk = (dst_nodal) ? src_stag[1] : dst_stag[1]; -#endif -#if defined(WARPX_DIM_3D) - const int sl = (dst_nodal) ? src_stag[2] : dst_stag[2]; -#endif + const int sk = (AMREX_SPACEDIM > 1) ? ((dst_nodal) ? src_stag[1] : dst_stag[1]) : 1; + const int sl = (AMREX_SPACEDIM > 2) ? ((dst_nodal) ? src_stag[2] : dst_stag[2]) : 1; // Interpolate along j,k,l only if source MultiFab is staggered along j,k,l const bool interp_j = (sj == 0); -#if (AMREX_SPACEDIM >= 2) const bool interp_k = (sk == 0); -#endif -#if defined(WARPX_DIM_3D) const bool interp_l = (sl == 0); -#endif -#if defined(WARPX_DIM_1D_Z) - const int noj = noz; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const int noj = nox; - const int nok = noz; -#elif defined(WARPX_DIM_3D) - const int noj = nox; - const int nok = noy; - const int nol = noz; -#endif + const int noj = AMREX_D_PICK(noz, nox, nox); + const int nok = AMREX_D_PICK(0 , noz, noy); + const int nol = AMREX_D_PICK(0 , 0 , noz); // Additional normalization factor const amrex::Real wj = (interp_j) ? 0.5_rt : 1.0_rt; -#if defined(WARPX_DIM_1D_Z) - constexpr amrex::Real wk = 1.0_rt; - constexpr amrex::Real wl = 1.0_rt; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt; - constexpr amrex::Real wl = 1.0_rt; -#elif defined(WARPX_DIM_3D) const amrex::Real wk = (interp_k) ? 0.5_rt : 1.0_rt; const amrex::Real wl = (interp_l) ? 0.5_rt : 1.0_rt; -#endif - // Min and max for interpolation loop along j + // Min and max for interpolation loop const int jmin = (interp_j) ? j - noj/2 + shift : j; const int jmax = (interp_j) ? j + noj/2 + shift - 1 : j; - - // Min and max for interpolation loop along k -#if defined(WARPX_DIM_1D_Z) - // k = 0 always - const int kmin = k; - const int kmax = k; -#else const int kmin = (interp_k) ? k - nok/2 + shift : k; const int kmax = (interp_k) ? k + nok/2 + shift - 1 : k; -#endif - - // Min and max for interpolation loop along l -#if (AMREX_SPACEDIM <= 2) - // l = 0 always - const int lmin = l; - const int lmax = l; -#elif defined(WARPX_DIM_3D) const int lmin = (interp_l) ? l - nol/2 + shift : l; const int lmax = (interp_l) ? l + nol/2 + shift - 1 : l; -#endif // Number of interpolation points const int nj = jmax - jmin; @@ -543,31 +434,16 @@ void warpx_interp (const int j, amrex::Real res = 0.0_rt; -#if defined(WARPX_DIM_1D_Z) - amrex::Real const* scj = stencil_coeffs_z; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::Real const* scj = stencil_coeffs_x; - amrex::Real const* sck = stencil_coeffs_z; -#elif defined(WARPX_DIM_3D) - amrex::Real const* scj = stencil_coeffs_x; - amrex::Real const* sck = stencil_coeffs_y; - amrex::Real const* scl = stencil_coeffs_z; -#endif + amrex::Real const* scj = AMREX_D_PICK(stencil_coeffs_z, stencil_coeffs_x, stencil_coeffs_x); + amrex::Real const* sck = AMREX_D_PICK(nullptr , stencil_coeffs_z, stencil_coeffs_y); + amrex::Real const* scl = AMREX_D_PICK(nullptr , nullptr , stencil_coeffs_z); for (int ll = 0; ll <= nl; ll++) { -#if defined(WARPX_DIM_3D) - const amrex::Real cl = (interp_l)? scl[ll] : 1.0_rt; -#else - const amrex::Real cl = 1.0_rt; -#endif + const amrex::Real cl = (interp_l)? scl[ll] : 1.0_rt; for (int kk = 0; kk <= nk; kk++) { -#if (AMREX_SPACEDIM >= 2) const amrex::Real ck = (interp_k)? sck[kk] : 1.0_rt; -#else - const amrex::Real ck = 1.0_rt; -#endif for (int jj = 0; jj <= nj; jj++) { const amrex::Real cj = (interp_j)? scj[jj] : 1.0_rt; diff --git a/Source/Particles/AddPlasmaUtilities.H b/Source/Particles/AddPlasmaUtilities.H new file mode 100644 index 00000000000..8f0489e3921 --- /dev/null +++ b/Source/Particles/AddPlasmaUtilities.H @@ -0,0 +1,210 @@ +/* Copyright 2024 The WarpX Community + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + * Authors: Andrew Myers + */ +#ifndef WARPX_ADDPLASMAUTILITIES_H_ +#define WARPX_ADDPLASMAUTILITIES_H_ + +#include "Initialization/PlasmaInjector.H" + +#ifdef WARPX_QED +# include "Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H" +# include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H" +#endif + +#include +#include +#include +#include +#include +#include + +/* + Finds the overlap region between the given tile_realbox and part_realbox, returning true + if an overlap exists and false if otherwise. This also sets the parameters overlap_realbox, + overlap_box, and shifted to the appropriate values. + */ +bool find_overlap (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted); + +/* + Finds the overlap region between the given tile_realbox, part_realbox and the surface used for + flux injection, returning true if an overlap exists and false if otherwise. This also sets the + parameters overlap_realbox, overlap_box, and shifted to the appropriate values. + */ +bool find_overlap_flux (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + const PlasmaInjector& plasma_injector, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted); + +/* + This computes the scale_fac (used for setting the particle weights) on a volumetric basis. + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real compute_scale_fac_volume (const amrex::GpuArray& dx, + const amrex::Long pcount) { + using namespace amrex::literals; + return (pcount != 0) ? AMREX_D_TERM(dx[0],*dx[1],*dx[2])/pcount : 0.0_rt; +} + +/* + Given a refinement ratio, this computes the total increase in resolution for a plane + defined by the normal_axis. + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +int compute_area_weights (const amrex::IntVect& iv, const int normal_axis) { + int r = AMREX_D_TERM(iv[0],*iv[1],*iv[2]); +#if defined(WARPX_DIM_3D) + r /= iv[normal_axis]; +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + if (normal_axis == 0) { r /= iv[0]; } + else if (normal_axis == 2) { r /= iv[1]; } +#elif defined(WARPX_DIM_1D_Z) + if (normal_axis == 2) { r /= iv[0]; } +#endif + return r; +} + +/* + This computes the scale_fac (used for setting the particle weights) on a on area basis + (used for flux injection). + */ +AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE +amrex::Real compute_scale_fac_area (const amrex::GpuArray& dx, + const amrex::Real num_ppc_real, const int flux_normal_axis) { + using namespace amrex::literals; + amrex::Real scale_fac = AMREX_D_TERM(dx[0],*dx[1],*dx[2])/num_ppc_real; + // Scale particle weight by the area of the emitting surface, within one cell +#if defined(WARPX_DIM_3D) + scale_fac /= dx[flux_normal_axis]; +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) + // When emission is in the r direction, the emitting surface is a cylinder. + // The factor 2*pi*r is added later below. + if (flux_normal_axis == 0) { scale_fac /= dx[0]; } + // When emission is in the z direction, the emitting surface is an annulus + // The factor 2*pi*r is added later below. + if (flux_normal_axis == 2) { scale_fac /= dx[1]; } + // When emission is in the theta direction (flux_normal_axis == 1), + // the emitting surface is a rectangle, within the plane of the simulation +#elif defined(WARPX_DIM_1D_Z) + if (flux_normal_axis == 2) { scale_fac /= dx[0]; } +#endif + return scale_fac; +} + +/* + These structs encapsulates several data structures needed for using the parser during plasma + injection. + */ +struct PlasmaParserWrapper +{ + PlasmaParserWrapper (std::size_t a_num_user_int_attribs, + std::size_t a_num_user_real_attribs, + const amrex::Vector< std::unique_ptr >& a_user_int_attrib_parser, + const amrex::Vector< std::unique_ptr >& a_user_real_attrib_parser); + + amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > m_user_int_attrib_parserexec_pinned; + amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > m_user_real_attrib_parserexec_pinned; +}; + +struct PlasmaParserHelper +{ + template + PlasmaParserHelper (SoAType& a_soa, std::size_t old_size, + const std::vector& a_user_int_attribs, + const std::vector& a_user_real_attribs, + std::map& a_particle_icomps, + std::map& a_particle_comps, + const PlasmaParserWrapper& wrapper) : + m_wrapper_ptr(&wrapper) { + m_pa_user_int_pinned.resize(a_user_int_attribs.size()); + m_pa_user_real_pinned.resize(a_user_real_attribs.size()); + +#ifdef AMREX_USE_GPU + m_d_pa_user_int.resize(a_user_int_attribs.size()); + m_d_pa_user_real.resize(a_user_real_attribs.size()); + m_d_user_int_attrib_parserexec.resize(a_user_int_attribs.size()); + m_d_user_real_attrib_parserexec.resize(a_user_real_attribs.size()); +#endif + + for (std::size_t ia = 0; ia < a_user_int_attribs.size(); ++ia) { + m_pa_user_int_pinned[ia] = a_soa.GetIntData(a_particle_icomps[a_user_int_attribs[ia]]).data() + old_size; + } + for (std::size_t ia = 0; ia < a_user_real_attribs.size(); ++ia) { + m_pa_user_real_pinned[ia] = a_soa.GetRealData(a_particle_comps[a_user_real_attribs[ia]]).data() + old_size; + } + +#ifdef AMREX_USE_GPU + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_pa_user_int_pinned.begin(), + m_pa_user_int_pinned.end(), m_d_pa_user_int.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, m_pa_user_real_pinned.begin(), + m_pa_user_real_pinned.end(), m_d_pa_user_real.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, wrapper.m_user_int_attrib_parserexec_pinned.begin(), + wrapper.m_user_int_attrib_parserexec_pinned.end(), m_d_user_int_attrib_parserexec.begin()); + amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, wrapper.m_user_real_attrib_parserexec_pinned.begin(), + wrapper.m_user_real_attrib_parserexec_pinned.end(), m_d_user_real_attrib_parserexec.begin()); +#endif + } + + int** getUserIntDataPtrs (); + amrex::ParticleReal** getUserRealDataPtrs (); + [[nodiscard]] amrex::ParserExecutor<7> const* getUserIntParserExecData () const; + [[nodiscard]] amrex::ParserExecutor<7> const* getUserRealParserExecData () const; + + amrex::Gpu::PinnedVector m_pa_user_int_pinned; + amrex::Gpu::PinnedVector m_pa_user_real_pinned; + +#ifdef AMREX_USE_GPU + // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, + // and them memcpy to device from host + amrex::Gpu::DeviceVector m_d_pa_user_int; + amrex::Gpu::DeviceVector m_d_pa_user_real; + amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > m_d_user_int_attrib_parserexec; + amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > m_d_user_real_attrib_parserexec; +#endif + const PlasmaParserWrapper* m_wrapper_ptr; +}; + +#ifdef WARPX_QED +struct QEDHelper +{ + template + QEDHelper (SoAType& a_soa, std::size_t old_size, + std::map& a_particle_comps, + bool a_has_quantum_sync, bool a_has_breit_wheeler, + const std::shared_ptr& a_shr_p_qs_engine, + const std::shared_ptr& a_shr_p_bw_engine) + : has_quantum_sync(a_has_quantum_sync), has_breit_wheeler(a_has_breit_wheeler) + { + if(has_quantum_sync){ + quantum_sync_get_opt = + a_shr_p_qs_engine->build_optical_depth_functor(); + p_optical_depth_QSR = a_soa.GetRealData( + a_particle_comps["opticalDepthQSR"]).data() + old_size; + } + if(has_breit_wheeler){ + breit_wheeler_get_opt = + a_shr_p_bw_engine->build_optical_depth_functor(); + p_optical_depth_BW = a_soa.GetRealData( + a_particle_comps["opticalDepthBW"]).data() + old_size; + } + } + + amrex::ParticleReal* p_optical_depth_QSR = nullptr; + amrex::ParticleReal* p_optical_depth_BW = nullptr; + + const bool has_quantum_sync; + const bool has_breit_wheeler; + + QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; + BreitWheelerGetOpticalDepth breit_wheeler_get_opt; +}; +#endif + +#endif /*WARPX_ADDPLASMAUTILITIES_H_*/ diff --git a/Source/Particles/AddPlasmaUtilities.cpp b/Source/Particles/AddPlasmaUtilities.cpp new file mode 100644 index 00000000000..31066516477 --- /dev/null +++ b/Source/Particles/AddPlasmaUtilities.cpp @@ -0,0 +1,158 @@ +/* Copyright 2024 The WarpX Community + * + * This file is part of WarpX. + * + * License: BSD-3-Clause-LBNL + * Authors: Andrew Myers + */ +#include "AddPlasmaUtilities.H" + +#include + +bool find_overlap (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted) +{ + using namespace amrex::literals; + + bool no_overlap = false; + for (int dir=0; dir= part_realbox.lo(dir) ) { + const amrex::Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); + } else { + no_overlap = true; break; + } + // Count the number of cells in this direction in overlap_realbox + overlap_box.setSmall( dir, 0 ); + overlap_box.setBig( dir, + int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) + /dx[dir] )) - 1); + shifted[dir] = + static_cast(std::round((overlap_realbox.lo(dir)-prob_lo[dir])/dx[dir])); + // shifted is exact in non-moving-window direction. That's all we care. + } + return no_overlap; +} + +bool find_overlap_flux (const amrex::RealBox& tile_realbox, const amrex::RealBox& part_realbox, + const amrex::GpuArray& dx, + const amrex::GpuArray& prob_lo, + const PlasmaInjector& plasma_injector, + amrex::RealBox& overlap_realbox, amrex::Box& overlap_box, amrex::IntVect& shifted) +{ + using namespace amrex::literals; + + bool no_overlap = false; + for (int dir=0; dir 0) { + if (plasma_injector.surface_flux_pos < tile_realbox.lo(dir) || + plasma_injector.surface_flux_pos >= tile_realbox.hi(dir)) { + no_overlap = true; + break; + } + } else { + if (plasma_injector.surface_flux_pos <= tile_realbox.lo(dir) || + plasma_injector.surface_flux_pos > tile_realbox.hi(dir)) { + no_overlap = true; + break; + } + } + overlap_realbox.setLo( dir, plasma_injector.surface_flux_pos ); + overlap_realbox.setHi( dir, plasma_injector.surface_flux_pos ); + overlap_box.setSmall( dir, 0 ); + overlap_box.setBig( dir, 0 ); + shifted[dir] = + static_cast(std::round((overlap_realbox.lo(dir)-prob_lo[dir])/dx[dir])); + } else { + if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { + const amrex::Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); + overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]); + } else { + no_overlap = true; break; + } + if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { + const amrex::Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); + overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); + } else { + no_overlap = true; break; + } + // Count the number of cells in this direction in overlap_realbox + overlap_box.setSmall( dir, 0 ); + overlap_box.setBig( dir, + int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) + /dx[dir] )) - 1); + shifted[dir] = + static_cast(std::round((overlap_realbox.lo(dir)-prob_lo[dir])/dx[dir])); + // shifted is exact in non-moving-window direction. That's all we care. + } + } + + return no_overlap; +} + +PlasmaParserWrapper::PlasmaParserWrapper (const std::size_t a_num_user_int_attribs, + const std::size_t a_num_user_real_attribs, + const amrex::Vector< std::unique_ptr >& a_user_int_attrib_parser, + const amrex::Vector< std::unique_ptr >& a_user_real_attrib_parser) + +{ + m_user_int_attrib_parserexec_pinned.resize(a_num_user_int_attribs); + m_user_real_attrib_parserexec_pinned.resize(a_num_user_real_attribs); + + for (std::size_t ia = 0; ia < a_num_user_int_attribs; ++ia) { + m_user_int_attrib_parserexec_pinned[ia] = a_user_int_attrib_parser[ia]->compile<7>(); + } + for (std::size_t ia = 0; ia < a_num_user_real_attribs; ++ia) { + m_user_real_attrib_parserexec_pinned[ia] = a_user_real_attrib_parser[ia]->compile<7>(); + } +} + +int** PlasmaParserHelper::getUserIntDataPtrs () { +#ifdef AMREX_USE_GPU + return m_d_pa_user_int.dataPtr(); +#else + return m_pa_user_int_pinned.dataPtr(); +#endif +} + +amrex::ParticleReal** PlasmaParserHelper::getUserRealDataPtrs () { +#ifdef AMREX_USE_GPU + return m_d_pa_user_real.dataPtr(); +#else + return m_pa_user_real_pinned.dataPtr(); +#endif +} + +amrex::ParserExecutor<7> const* PlasmaParserHelper::getUserIntParserExecData () const { +#ifdef AMREX_USE_GPU + return m_d_user_int_attrib_parserexec.dataPtr(); +#else + return m_wrapper_ptr->m_user_int_attrib_parserexec_pinned.dataPtr(); +#endif +} + +amrex::ParserExecutor<7> const* PlasmaParserHelper::getUserRealParserExecData () const { +#ifdef AMREX_USE_GPU + return m_d_user_real_attrib_parserexec.dataPtr(); +#else + return m_wrapper_ptr->m_user_real_attrib_parserexec_pinned.dataPtr(); +#endif +} diff --git a/Source/Particles/CMakeLists.txt b/Source/Particles/CMakeLists.txt index 67af14ef889..6b434c0a4e1 100644 --- a/Source/Particles/CMakeLists.txt +++ b/Source/Particles/CMakeLists.txt @@ -2,6 +2,7 @@ foreach(D IN LISTS WarpX_DIMS) warpx_set_suffix_dims(SD ${D}) target_sources(lib_${SD} PRIVATE + AddPlasmaUtilities.cpp MultiParticleContainer.cpp ParticleBoundaries.cpp PhotonParticleContainer.cpp diff --git a/Source/Particles/Gather/FieldGather.H b/Source/Particles/Gather/FieldGather.H index 4b4590b8642..95a010031ec 100644 --- a/Source/Particles/Gather/FieldGather.H +++ b/Source/Particles/Gather/FieldGather.H @@ -1714,9 +1714,9 @@ void doGatherShapeNImplicit ( const amrex::Dim3& lo, const int n_rz_azimuthal_modes, const int nox, - const int depos_type ) + const CurrentDepositionAlgo depos_type ) { - if (depos_type==0) { // CurrentDepositionAlgo::Esirkepov + if (depos_type == CurrentDepositionAlgo::Esirkepov) { if (nox == 1) { doGatherShapeNEsirkepovStencilImplicit<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, @@ -1743,7 +1743,7 @@ void doGatherShapeNImplicit ( dinv, xyzmin, lo, n_rz_azimuthal_modes); } } - else if (depos_type==3) { // CurrentDepositionAlgo::Villasenor + else if (depos_type == CurrentDepositionAlgo::Villasenor) { if (nox == 1) { doGatherPicnicShapeN<1>(xp_n, yp_n, zp_n, xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, @@ -1770,7 +1770,7 @@ void doGatherShapeNImplicit ( dinv, xyzmin, lo, n_rz_azimuthal_modes); } } - else if (depos_type==1) { // CurrentDepositionAlgo::Direct + else if (depos_type == CurrentDepositionAlgo::Direct) { if (nox == 1) { doGatherShapeN<1,0>(xp_nph, yp_nph, zp_nph, Exp, Eyp, Ezp, Bxp, Byp, Bzp, ex_arr, ey_arr, ez_arr, bx_arr, by_arr, bz_arr, diff --git a/Source/Particles/LaserParticleContainer.cpp b/Source/Particles/LaserParticleContainer.cpp index ee3954546c7..234b32d5359 100644 --- a/Source/Particles/LaserParticleContainer.cpp +++ b/Source/Particles/LaserParticleContainer.cpp @@ -413,12 +413,10 @@ LaserParticleContainer::InitData (int lev) #if defined(WARPX_DIM_3D) return {m_u_X[0]*(pos[0]-m_position[0])+m_u_X[1]*(pos[1]-m_position[1])+m_u_X[2]*(pos[2]-m_position[2]), m_u_Y[0]*(pos[0]-m_position[0])+m_u_Y[1]*(pos[1]-m_position[1])+m_u_Y[2]*(pos[2]-m_position[2])}; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) -# if defined(WARPX_DIM_RZ) +#elif defined(WARPX_DIM_RZ) return {pos[0]-m_position[0], 0.0_rt}; -# else +#elif defined(WARPX_DIM_XZ) return {m_u_X[0]*(pos[0]-m_position[0])+m_u_X[2]*(pos[2]-m_position[2]), 0.0_rt}; -# endif #else return {m_u_X[2]*(pos[2]-m_position[2]), 0.0_rt}; #endif @@ -743,13 +741,12 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const Sy = std::min(std::min(dx[0]/(std::abs(m_u_Y[0])+eps), dx[1]/(std::abs(m_u_Y[1])+eps)), dx[2]/(std::abs(m_u_Y[2])+eps)); -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) -# if defined(WARPX_DIM_RZ) +#elif defined(WARPX_DIM_RZ) Sx = dx[0]; -# else + Sy = 1.0; +#elif defined(WARPX_DIM_XZ) Sx = std::min(dx[0]/(std::abs(m_u_X[0])+eps), dx[2]/(std::abs(m_u_X[2])+eps)); -# endif Sy = 1.0; #else Sx = 1.0; @@ -759,7 +756,7 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const } void -LaserParticleContainer::ComputeWeightMobility (Real Sx, Real Sy) +LaserParticleContainer::ComputeWeightMobility ([[maybe_unused]] Real Sx, [[maybe_unused]] Real Sy) { // The mobility is the constant of proportionality between the field to // be emitted, and the corresponding velocity that the particles need to have. @@ -769,14 +766,7 @@ LaserParticleContainer::ComputeWeightMobility (Real Sx, Real Sy) m_mobility = eps/m_e_max; m_weight = PhysConst::ep0 / m_mobility; // Multiply by particle spacing -#if defined(WARPX_DIM_3D) - m_weight *= Sx * Sy; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - m_weight *= Sx; - amrex::ignore_unused(Sy); -#else - amrex::ignore_unused(Sx,Sy); -#endif + m_weight *= AMREX_D_TERM(1._rt, * Sx, * Sy); // When running in the boosted-frame, the input parameters (and in particular // the amplitude of the field) are given in the lab-frame. // Therefore, the mobility needs to be modified by a factor WarpX::gamma_boost. diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 58cbe11a980..69918f69940 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -1,3 +1,4 @@ +CEXE_sources += AddPlasmaUtilities.cpp CEXE_sources += MultiParticleContainer.cpp CEXE_sources += WarpXParticleContainer.cpp CEXE_sources += RigidInjectedParticleContainer.cpp diff --git a/Source/Particles/ParticleIO.H b/Source/Particles/ParticleIO.H index d5fc68f4097..8d3516e6890 100644 --- a/Source/Particles/ParticleIO.H +++ b/Source/Particles/ParticleIO.H @@ -90,6 +90,6 @@ particlesConvertUnits (ConvertDirection convert_direction, T_ParticleContainer* */ void storePhiOnParticles ( PinnedMemoryParticleContainer& tmp, - int electrostatic_solver_id, bool is_full_diagnostic ); + ElectrostaticSolverAlgo electrostatic_solver_id, bool is_full_diagnostic ); #endif /* WARPX_PARTICLEIO_H_ */ diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b3e9f258816..5f5192ba262 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -356,6 +356,14 @@ public: } protected: + + /* + Finds the box defining the region where refine injection should be used, if that + option is enabled. Currently this only works for numLevels() == 2 and static mesh + refinement. + */ + bool findRefinedInjectionBox (amrex::Box& fine_injection_box, amrex::IntVect& rrfac); + std::string species_name; std::vector> plasma_injectors; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a5deb3640c3..bf77c0b2291 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -15,6 +15,7 @@ #include "Initialization/InjectorMomentum.H" #include "Initialization/InjectorPosition.H" #include "MultiParticleContainer.H" +#include "Particles/AddPlasmaUtilities.H" #ifdef WARPX_QED # include "Particles/ElementaryProcess/QEDInternals/BreitWheelerEngineWrapper.H" # include "Particles/ElementaryProcess/QEDInternals/QuantumSyncEngineWrapper.H" @@ -178,16 +179,16 @@ namespace pos.x = lo_corner[0] + (iv[0]+r.x)*dx[0]; pos.y = lo_corner[1] + (iv[1]+r.y)*dx[1]; pos.z = lo_corner[2] + (iv[2]+r.z)*dx[2]; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) +#elif defined(WARPX_DIM_XZ) pos.x = lo_corner[0] + (iv[0]+r.x)*dx[0]; pos.y = 0.0_rt; -#if defined WARPX_DIM_XZ pos.z = lo_corner[1] + (iv[1]+r.y)*dx[1]; -#elif defined WARPX_DIM_RZ +#elif defined(WARPX_DIM_RZ) // Note that for RZ, r.y will be theta + pos.x = lo_corner[0] + (iv[0]+r.x)*dx[0]; + pos.y = 0.0_rt; pos.z = lo_corner[1] + (iv[1]+r.z)*dx[1]; -#endif -#else +#elif defined(WARPX_DIM_1D_Z) pos.x = 0.0_rt; pos.y = 0.0_rt; pos.z = lo_corner[0] + (iv[0]+r.x)*dx[0]; @@ -217,29 +218,17 @@ namespace const GpuArray& pa, long& ip, const bool& do_field_ionization, int* pi #ifdef WARPX_QED - ,const bool& has_quantum_sync, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_QSR - ,const bool& has_breit_wheeler, amrex::ParticleReal* AMREX_RESTRICT p_optical_depth_BW + ,const QEDHelper& qed_helper #endif ) noexcept { - pa[PIdx::z][ip] = 0._rt; -#if (AMREX_SPACEDIM >= 2) - pa[PIdx::x][ip] = 0._rt; -#endif -#if defined(WARPX_DIM_3D) - pa[PIdx::y][ip] = 0._rt; -#endif - pa[PIdx::w ][ip] = 0._rt; - pa[PIdx::ux][ip] = 0._rt; - pa[PIdx::uy][ip] = 0._rt; - pa[PIdx::uz][ip] = 0._rt; -#ifdef WARPX_DIM_RZ - pa[PIdx::theta][ip] = 0._rt; -#endif + for (int idx=0 ; idx < PIdx::nattribs ; idx++) { + pa[idx][ip] = 0._rt; + } if (do_field_ionization) {pi[ip] = 0;} #ifdef WARPX_QED - if (has_quantum_sync) {p_optical_depth_QSR[ip] = 0._rt;} - if (has_breit_wheeler) {p_optical_depth_BW[ip] = 0._rt;} + if (qed_helper.has_quantum_sync) {qed_helper.p_optical_depth_QSR[ip] = 0._rt;} + if (qed_helper.has_breit_wheeler) {qed_helper.p_optical_depth_BW[ip] = 0._rt;} #endif idcpu[ip] = amrex::ParticleIdCpus::Invalid; @@ -758,6 +747,7 @@ PhysicalParticleContainer::AddPlasmaFromFile(PlasmaInjector & plasma_injector, const std::shared_ptr ptr_offset_z = ps["positionOffset"]["z"].loadChunk(); auto const position_unit_z = static_cast(ps["position"]["z"].unitSI()); auto const position_offset_unit_z = static_cast(ps["positionOffset"]["z"].unitSI()); + const std::shared_ptr ptr_ux = ps["momentum"]["x"].loadChunk(); auto const momentum_unit_x = static_cast(ps["momentum"]["x"].unitSI()); const std::shared_ptr ptr_uz = ps["momentum"]["z"].loadChunk(); @@ -974,22 +964,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int amrex::LayoutData* cost = WarpX::getCosts(lev); - const int nlevs = numLevels(); - static bool refine_injection = false; - static Box fine_injection_box; - static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); - // This does not work if the mesh is dynamic. But in that case, we should - // not use refined injected either. We also assume there is only one fine level. - if (WarpX::moving_window_active(WarpX::GetInstance().getistep(0)+1) and WarpX::refine_plasma - and do_continuous_injection and nlevs == 2) - { - refine_injection = true; - fine_injection_box = ParticleBoxArray(1).minimalBox(); - fine_injection_box.setSmall(WarpX::moving_window_dir, std::numeric_limits::lowest()/2); - fine_injection_box.setBig(WarpX::moving_window_dir, std::numeric_limits::max()/2); - rrfac = m_gdb->refRatio(0); - fine_injection_box.coarsen(rrfac); - } + Box fine_injection_box; + amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + const bool refine_injection = findRefinedInjectionBox(fine_injection_box, rrfac); InjectorPosition* inj_pos = plasma_injector.getInjectorPosition(); InjectorDensity* inj_rho = plasma_injector.getInjectorDensity(); @@ -1005,18 +982,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int const bool radially_weighted = plasma_injector.radially_weighted; #endif - - // User-defined integer and real attributes: prepare parsers - const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); - const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_int_attrib_parserexec_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_real_attrib_parserexec_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - user_int_attrib_parserexec_pinned[ia] = m_user_int_attrib_parser[ia]->compile<7>(); - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - user_real_attrib_parserexec_pinned[ia] = m_user_real_attrib_parser[ia]->compile<7>(); - } + auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); + auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); + const PlasmaParserWrapper plasma_parser_wrapper (m_user_int_attribs.size(), + m_user_real_attribs.size(), + m_user_int_attrib_parser, + m_user_real_attrib_parser); MFItInfo info; if (do_tiling && Gpu::notInLaunchRegion()) { @@ -1042,30 +1013,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int RealBox overlap_realbox; Box overlap_box; IntVect shifted; - bool no_overlap = false; - - for (int dir=0; dir= part_realbox.lo(dir) ) { - const Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - // Count the number of cells in this direction in overlap_realbox - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, - int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) - /dx[dir] )) - 1); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - // shifted is exact in non-moving-window direction. That's all we care. - } + const bool no_overlap = find_overlap(tile_realbox, part_realbox, dx, problo, overlap_realbox, overlap_box, shifted); if (no_overlap) { continue; // Go to the next tile } @@ -1120,19 +1068,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } return 0; }; - const int flag_pcount = checker(); - if (flag_pcount == 1) { - pcounts[index] = num_ppc*r; - } else { - pcounts[index] = 0; - } + pcounts[index] = checker() ? num_ppc*r : 0; } -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(k); -#endif -#if defined(WARPX_DIM_1D_Z) amrex::ignore_unused(j,k); -#endif }); // Max number of new particles. All of them are created, @@ -1170,40 +1108,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int pa[ia] = soa.GetRealData(ia).data() + old_size; } uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; - // user-defined integer and real attributes - amrex::Gpu::PinnedVector pa_user_int_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector pa_user_real_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - pa_user_int_pinned[ia] = soa.GetIntData(particle_icomps[m_user_int_attribs[ia]]).data() + old_size; - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - pa_user_real_pinned[ia] = soa.GetRealData(particle_comps[m_user_real_attribs[ia]]).data() + old_size; - } -#ifdef AMREX_USE_GPU - // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, - // and them memcpy to device from host - amrex::Gpu::DeviceVector d_pa_user_int(n_user_int_attribs); - amrex::Gpu::DeviceVector d_pa_user_real(n_user_real_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_int_attrib_parserexec(n_user_int_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_real_attrib_parserexec(n_user_real_attribs); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_int_pinned.begin(), - pa_user_int_pinned.end(), d_pa_user_int.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_real_pinned.begin(), - pa_user_real_pinned.end(), d_pa_user_real.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_int_attrib_parserexec_pinned.begin(), - user_int_attrib_parserexec_pinned.end(), d_user_int_attrib_parserexec.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_real_attrib_parserexec_pinned.begin(), - user_real_attrib_parserexec_pinned.end(), d_user_real_attrib_parserexec.begin()); - int** pa_user_int_data = d_pa_user_int.dataPtr(); - ParticleReal** pa_user_real_data = d_pa_user_real.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = d_user_int_attrib_parserexec.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = d_user_real_attrib_parserexec.dataPtr(); -#else - int** pa_user_int_data = pa_user_int_pinned.dataPtr(); - ParticleReal** pa_user_real_data = pa_user_real_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = user_int_attrib_parserexec_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = user_real_attrib_parserexec_pinned.dataPtr(); -#endif + + PlasmaParserHelper plasma_parser_helper (soa, old_size, m_user_int_attribs, m_user_real_attribs, particle_icomps, particle_comps, plasma_parser_wrapper); + int** pa_user_int_data = plasma_parser_helper.getUserIntDataPtrs(); + ParticleReal** pa_user_real_data = plasma_parser_helper.getUserRealDataPtrs(); + amrex::ParserExecutor<7> const* user_int_parserexec_data = plasma_parser_helper.getUserIntParserExecData(); + amrex::ParserExecutor<7> const* user_real_parserexec_data = plasma_parser_helper.getUserRealParserExecData(); int* pi = nullptr; if (do_field_ionization) { @@ -1211,34 +1121,9 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } #ifdef WARPX_QED - //Pointer to the optical depth component - amrex::ParticleReal* p_optical_depth_QSR = nullptr; - amrex::ParticleReal* p_optical_depth_BW = nullptr; - - // If a QED effect is enabled, the corresponding optical depth - // has to be initialized - const bool loc_has_quantum_sync = has_quantum_sync(); - const bool loc_has_breit_wheeler = has_breit_wheeler(); - if (loc_has_quantum_sync) { - p_optical_depth_QSR = soa.GetRealData( - particle_comps["opticalDepthQSR"]).data() + old_size; - } - if(loc_has_breit_wheeler) { - p_optical_depth_BW = soa.GetRealData( - particle_comps["opticalDepthBW"]).data() + old_size; - } - - //If needed, get the appropriate functors from the engines - QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; - BreitWheelerGetOpticalDepth breit_wheeler_get_opt; - if(loc_has_quantum_sync){ - quantum_sync_get_opt = - m_shr_p_qs_engine->build_optical_depth_functor(); - } - if(loc_has_breit_wheeler){ - breit_wheeler_get_opt = - m_shr_p_bw_engine->build_optical_depth_functor(); - } + const QEDHelper qed_helper(soa, old_size, particle_comps, + has_quantum_sync(), has_breit_wheeler(), + m_shr_p_qs_engine, m_shr_p_bw_engine); #endif const bool loc_do_field_ionization = do_field_ionization; @@ -1256,23 +1141,14 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int [=] AMREX_GPU_DEVICE (int i, int j, int k, amrex::RandomEngine const& engine) noexcept { const IntVect iv = IntVect(AMREX_D_DECL(i, j, k)); + amrex::ignore_unused(j,k); const auto index = overlap_box.index(iv); #ifdef WARPX_DIM_RZ Real theta_offset = 0._rt; if (rz_random_theta) { theta_offset = amrex::Random(engine) * 2._rt * MathConst::pi; } #endif - Real scale_fac = 0.0_rt; - if( pcounts[index] != 0) { -#if defined(WARPX_DIM_3D) - scale_fac = dx[0]*dx[1]*dx[2]/pcounts[index]; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - scale_fac = dx[0]*dx[1]/pcounts[index]; -#elif defined(WARPX_DIM_1D_Z) - scale_fac = dx[0]/pcounts[index]; -#endif - } - + const Real scale_fac = compute_scale_fac_volume(dx, pcounts[index]); for (int i_part = 0; i_part < pcounts[index]; ++i_part) { long ip = poffset[index] + i_part; @@ -1285,38 +1161,22 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int auto pos = getCellCoords(overlap_corner, dx, r, iv); #if defined(WARPX_DIM_3D) - if (!tile_realbox.contains(XDim3{pos.x,pos.y,pos.z})) { - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi -#ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW -#endif - ); - continue; - } + bool const box_contains = tile_realbox.contains(XDim3{pos.x,pos.y,pos.z}); #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) amrex::ignore_unused(k); - if (!tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt})) { - ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi -#ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW -#endif - ); - continue; - } -#else + bool const box_contains = tile_realbox.contains(XDim3{pos.x,pos.z,0.0_rt}); +#elif defined(WARPX_DIM_1D_Z) amrex::ignore_unused(j,k); - if (!tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt})) { + bool const box_contains = tile_realbox.contains(XDim3{pos.z,0.0_rt,0.0_rt}); +#endif + if (!box_contains) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; } -#endif // Save the x and y values to use in the insideBounds checks. // This is needed with WARPX_DIM_RZ since x and y are modified. @@ -1349,8 +1209,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!inj_pos->insideBounds(xb, yb, z0)) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1363,8 +1222,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if ( dens < density_min ){ ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1381,8 +1239,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if (!inj_pos->insideBounds(xb, yb, z0_lab)) { ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1393,8 +1250,7 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int if ( dens < density_min ){ ZeroInitializeAndSetNegativeID(pa_idcpu, pa, ip, loc_do_field_ionization, pi #ifdef WARPX_QED - ,loc_has_quantum_sync, p_optical_depth_QSR - ,loc_has_breit_wheeler, p_optical_depth_BW + ,qed_helper #endif ); continue; @@ -1418,12 +1274,12 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int } #ifdef WARPX_QED - if(loc_has_quantum_sync){ - p_optical_depth_QSR[ip] = quantum_sync_get_opt(engine); + if(qed_helper.has_quantum_sync){ + qed_helper.p_optical_depth_QSR[ip] = qed_helper.quantum_sync_get_opt(engine); } - if(loc_has_breit_wheeler){ - p_optical_depth_BW[ip] = breit_wheeler_get_opt(engine); + if(qed_helper.has_breit_wheeler){ + qed_helper.p_optical_depth_BW[ip] = qed_helper.breit_wheeler_get_opt(engine); } #endif // Initialize user-defined integers with user-defined parser @@ -1461,13 +1317,14 @@ PhysicalParticleContainer::AddPlasma (PlasmaInjector const& plasma_injector, int pa[PIdx::x][ip] = pos.x; pa[PIdx::y][ip] = pos.y; pa[PIdx::z][ip] = pos.z; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) -#ifdef WARPX_DIM_RZ +#elif defined(WARPX_DIM_XZ) + pa[PIdx::x][ip] = pos.x; + pa[PIdx::z][ip] = pos.z; +#elif defined(WARPX_DIM_RZ) pa[PIdx::theta][ip] = theta; -#endif pa[PIdx::x][ip] = xb; pa[PIdx::z][ip] = pos.z; -#else +#elif defined(WARPX_DIM_1D_Z) pa[PIdx::z][ip] = pos.z; #endif } @@ -1513,25 +1370,6 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const auto dx = geom.CellSizeArray(); const auto problo = geom.ProbLoArray(); - Real scale_fac = 0._rt; - // Scale particle weight by the area of the emitting surface, within one cell -#if defined(WARPX_DIM_3D) - scale_fac = dx[0]*dx[1]*dx[2]/dx[plasma_injector.flux_normal_axis]/num_ppc_real; -#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - scale_fac = dx[0]*dx[1]/num_ppc_real; - // When emission is in the r direction, the emitting surface is a cylinder. - // The factor 2*pi*r is added later below. - if (plasma_injector.flux_normal_axis == 0) { scale_fac /= dx[0]; } - // When emission is in the z direction, the emitting surface is an annulus - // The factor 2*pi*r is added later below. - if (plasma_injector.flux_normal_axis == 2) { scale_fac /= dx[1]; } - // When emission is in the theta direction (flux_normal_axis == 1), - // the emitting surface is a rectangle, within the plane of the simulation -#elif defined(WARPX_DIM_1D_Z) - scale_fac = dx[0]/num_ppc_real; - if (plasma_injector.flux_normal_axis == 2) { scale_fac /= dx[0]; } -#endif - amrex::LayoutData* cost = WarpX::getCosts(0); // Create temporary particle container to which particles will be added; @@ -1542,19 +1380,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, for (int ic = 0; ic < NumRuntimeIntComps(); ++ic) { tmp_pc.AddIntComp(false); } tmp_pc.defineAllParticleTiles(); - const int nlevs = numLevels(); - static bool refine_injection = false; - static Box fine_injection_box; - static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); - // This does not work if the mesh is dynamic. But in that case, we should - // not use refined injected either. We also assume there is only one fine level. - if (WarpX::refine_plasma && nlevs == 2) - { - refine_injection = true; - fine_injection_box = ParticleBoxArray(1).minimalBox(); - rrfac = m_gdb->refRatio(0); - fine_injection_box.coarsen(rrfac); - } + Box fine_injection_box; + amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + const bool refine_injection = findRefinedInjectionBox(fine_injection_box, rrfac); InjectorPosition* flux_pos = plasma_injector.getInjectorFluxPosition(); InjectorFlux* inj_flux = plasma_injector.getInjectorFlux(); @@ -1568,6 +1396,13 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, const bool radially_weighted = plasma_injector.radially_weighted; #endif + auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); + auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); + const PlasmaParserWrapper plasma_parser_wrapper (m_user_int_attribs.size(), + m_user_real_attribs.size(), + m_user_int_attrib_parser, + m_user_real_attrib_parser); + MFItInfo info; if (do_tiling && Gpu::notInLaunchRegion()) { info.EnableTiling(tile_size); @@ -1592,61 +1427,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, RealBox overlap_realbox; Box overlap_box; IntVect shifted; - bool no_overlap = false; - - for (int dir=0; dir 0) { - if (plasma_injector.surface_flux_pos < tile_realbox.lo(dir) || - plasma_injector.surface_flux_pos >= tile_realbox.hi(dir)) { - no_overlap = true; - break; - } - } else { - if (plasma_injector.surface_flux_pos <= tile_realbox.lo(dir) || - plasma_injector.surface_flux_pos > tile_realbox.hi(dir)) { - no_overlap = true; - break; - } - } - overlap_realbox.setLo( dir, plasma_injector.surface_flux_pos ); - overlap_realbox.setHi( dir, plasma_injector.surface_flux_pos ); - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, 0 ); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - } else { - if ( tile_realbox.lo(dir) <= part_realbox.hi(dir) ) { - const Real ncells_adjust = std::floor( (tile_realbox.lo(dir) - part_realbox.lo(dir))/dx[dir] ); - overlap_realbox.setLo( dir, part_realbox.lo(dir) + std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - if ( tile_realbox.hi(dir) >= part_realbox.lo(dir) ) { - const Real ncells_adjust = std::floor( (part_realbox.hi(dir) - tile_realbox.hi(dir))/dx[dir] ); - overlap_realbox.setHi( dir, part_realbox.hi(dir) - std::max(ncells_adjust, 0._rt) * dx[dir]); - } else { - no_overlap = true; break; - } - // Count the number of cells in this direction in overlap_realbox - overlap_box.setSmall( dir, 0 ); - overlap_box.setBig( dir, - int( std::round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir)) - /dx[dir] )) - 1); - shifted[dir] = - static_cast(std::round((overlap_realbox.lo(dir)-problo[dir])/dx[dir])); - // shifted is exact in non-moving-window direction. That's all we care. - } - } + const bool no_overlap = find_overlap_flux(tile_realbox, part_realbox, dx, problo, plasma_injector, overlap_realbox, overlap_box, shifted); if (no_overlap) { continue; // Go to the next tile } @@ -1664,6 +1445,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, Gpu::DeviceVector offset(overlap_box.numPts()); auto *pcounts = counts.data(); const amrex::IntVect lrrfac = rrfac; + const int flux_normal_axis = plasma_injector.flux_normal_axis; Box fine_overlap_box; // default Box is NOT ok(). if (refine_injection) { fine_overlap_box = overlap_box & amrex::shift(fine_injection_box, -shifted); @@ -1681,17 +1463,13 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, auto index = overlap_box.index(iv); int r; if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { - r = AMREX_D_TERM(lrrfac[0],*lrrfac[1],*lrrfac[2]); + r = compute_area_weights(lrrfac, flux_normal_axis); } else { r = 1; } pcounts[index] = num_ppc_int*r; } -#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - amrex::ignore_unused(k); -#elif defined(WARPX_DIM_1D_Z) amrex::ignore_unused(j,k); -#endif }); // Max number of new particles. All of them are created, @@ -1726,46 +1504,11 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } uint64_t * AMREX_RESTRICT pa_idcpu = soa.GetIdCPUData().data() + old_size; - // user-defined integer and real attributes - const auto n_user_int_attribs = static_cast(m_user_int_attribs.size()); - const auto n_user_real_attribs = static_cast(m_user_real_attribs.size()); - amrex::Gpu::PinnedVector pa_user_int_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector pa_user_real_pinned(n_user_real_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_int_attrib_parserexec_pinned(n_user_int_attribs); - amrex::Gpu::PinnedVector< amrex::ParserExecutor<7> > user_real_attrib_parserexec_pinned(n_user_real_attribs); - for (int ia = 0; ia < n_user_int_attribs; ++ia) { - pa_user_int_pinned[ia] = soa.GetIntData(particle_icomps[m_user_int_attribs[ia]]).data() + old_size; - user_int_attrib_parserexec_pinned[ia] = m_user_int_attrib_parser[ia]->compile<7>(); - } - for (int ia = 0; ia < n_user_real_attribs; ++ia) { - pa_user_real_pinned[ia] = soa.GetRealData(particle_comps[m_user_real_attribs[ia]]).data() + old_size; - user_real_attrib_parserexec_pinned[ia] = m_user_real_attrib_parser[ia]->compile<7>(); - } -#ifdef AMREX_USE_GPU - // To avoid using managed memory, we first define pinned memory vector, initialize on cpu, - // and them memcpy to device from host - amrex::Gpu::DeviceVector d_pa_user_int(n_user_int_attribs); - amrex::Gpu::DeviceVector d_pa_user_real(n_user_real_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_int_attrib_parserexec(n_user_int_attribs); - amrex::Gpu::DeviceVector< amrex::ParserExecutor<7> > d_user_real_attrib_parserexec(n_user_real_attribs); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_int_pinned.begin(), - pa_user_int_pinned.end(), d_pa_user_int.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, pa_user_real_pinned.begin(), - pa_user_real_pinned.end(), d_pa_user_real.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_int_attrib_parserexec_pinned.begin(), - user_int_attrib_parserexec_pinned.end(), d_user_int_attrib_parserexec.begin()); - amrex::Gpu::copyAsync(Gpu::hostToDevice, user_real_attrib_parserexec_pinned.begin(), - user_real_attrib_parserexec_pinned.end(), d_user_real_attrib_parserexec.begin()); - int** pa_user_int_data = d_pa_user_int.dataPtr(); - ParticleReal** pa_user_real_data = d_pa_user_real.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = d_user_int_attrib_parserexec.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = d_user_real_attrib_parserexec.dataPtr(); -#else - int** pa_user_int_data = pa_user_int_pinned.dataPtr(); - ParticleReal** pa_user_real_data = pa_user_real_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_int_parserexec_data = user_int_attrib_parserexec_pinned.dataPtr(); - amrex::ParserExecutor<7> const* user_real_parserexec_data = user_real_attrib_parserexec_pinned.dataPtr(); -#endif + PlasmaParserHelper plasma_parser_helper (soa, old_size, m_user_int_attribs, m_user_real_attribs, particle_icomps, particle_comps, plasma_parser_wrapper); + int** pa_user_int_data = plasma_parser_helper.getUserIntDataPtrs(); + ParticleReal** pa_user_real_data = plasma_parser_helper.getUserRealDataPtrs(); + amrex::ParserExecutor<7> const* user_int_parserexec_data = plasma_parser_helper.getUserIntParserExecData(); + amrex::ParserExecutor<7> const* user_real_parserexec_data = plasma_parser_helper.getUserRealParserExecData(); int* p_ion_level = nullptr; if (do_field_ionization) { @@ -1773,34 +1516,9 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } #ifdef WARPX_QED - //Pointer to the optical depth component - amrex::ParticleReal* p_optical_depth_QSR = nullptr; - amrex::ParticleReal* p_optical_depth_BW = nullptr; - - // If a QED effect is enabled, the corresponding optical depth - // has to be initialized - const bool loc_has_quantum_sync = has_quantum_sync(); - const bool loc_has_breit_wheeler = has_breit_wheeler(); - if (loc_has_quantum_sync) { - p_optical_depth_QSR = soa.GetRealData( - particle_comps["opticalDepthQSR"]).data() + old_size; - } - if(loc_has_breit_wheeler) { - p_optical_depth_BW = soa.GetRealData( - particle_comps["opticalDepthBW"]).data() + old_size; - } - - //If needed, get the appropriate functors from the engines - QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; - BreitWheelerGetOpticalDepth breit_wheeler_get_opt; - if(loc_has_quantum_sync){ - quantum_sync_get_opt = - m_shr_p_qs_engine->build_optical_depth_functor(); - } - if(loc_has_breit_wheeler){ - breit_wheeler_get_opt = - m_shr_p_bw_engine->build_optical_depth_functor(); - } + const QEDHelper qed_helper(soa, old_size, particle_comps, + has_quantum_sync(), has_breit_wheeler(), + m_shr_p_qs_engine, m_shr_p_bw_engine); #endif const bool loc_do_field_ionization = do_field_ionization; @@ -1819,6 +1537,24 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, { const IntVect iv = IntVect(AMREX_D_DECL(i, j, k)); const auto index = overlap_box.index(iv); + + Real scale_fac = compute_scale_fac_area(dx, num_ppc_real, flux_normal_axis); + + auto lo = getCellCoords(overlap_corner, dx, {0._rt, 0._rt, 0._rt}, iv); + auto hi = getCellCoords(overlap_corner, dx, {1._rt, 1._rt, 1._rt}, iv); + + if (flux_pos->overlapsWith(lo, hi)) + { + int r; + if (fine_overlap_box.ok() && fine_overlap_box.contains(iv)) { + r = compute_area_weights(lrrfac, flux_normal_axis); + } else { + r = 1; + } + scale_fac /= r; + } + amrex::ignore_unused(j,k); + for (int i_part = 0; i_part < pcounts[index]; ++i_part) { const long ip = poffset[index] + i_part; @@ -1910,14 +1646,15 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, } #ifdef WARPX_QED - if (loc_has_quantum_sync) { - p_optical_depth_QSR[ip] = quantum_sync_get_opt(engine); + if(qed_helper.has_quantum_sync){ + qed_helper.p_optical_depth_QSR[ip] = qed_helper.quantum_sync_get_opt(engine); } - if(loc_has_breit_wheeler){ - p_optical_depth_BW[ip] = breit_wheeler_get_opt(engine); + if(qed_helper.has_breit_wheeler){ + qed_helper.p_optical_depth_BW[ip] = qed_helper.breit_wheeler_get_opt(engine); } #endif + // Initialize user-defined integers with user-defined parser for (int ia = 0; ia < n_user_int_attribs; ++ia) { pa_user_int_data[ia][ip] = static_cast(user_int_parserexec_data[ia](pos.x, pos.y, pos.z, u.x, u.y, u.z, t)); @@ -1970,7 +1707,7 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, #elif defined(WARPX_DIM_XZ) pa[PIdx::x][ip] = ppos.x; pa[PIdx::z][ip] = ppos.z; -#else +#elif defined(WARPX_DIM_1D_Z) pa[PIdx::z][ip] = ppos.z; #endif } @@ -2002,26 +1739,8 @@ PhysicalParticleContainer::AddPlasmaFlux (PlasmaInjector const& plasma_injector, // are in the right tile.) tmp_pc.Redistribute(); - // Add the particles to the current container, tile by tile - for (int lev=0; levaddParticles(tmp_pc, true); } void @@ -2397,13 +2116,7 @@ PhysicalParticleContainer::SplitParticles (int lev) long np_split; if(split_type==0) { - #if defined(WARPX_DIM_3D) - np_split = 8; - #elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - np_split = 4; - #else - np_split = 2; - #endif + np_split = amrex::Math::powi(2); } else { np_split = 2*AMREX_SPACEDIM; } @@ -2862,15 +2575,12 @@ PhysicalParticleContainer::PushPX (WarpXParIter& pti, if (save_previous_position) { #if (AMREX_SPACEDIM >= 2) x_old = pti.GetAttribs(particle_comps["prev_x"]).dataPtr() + offset; -#else - amrex::ignore_unused(x_old); #endif #if defined(WARPX_DIM_3D) y_old = pti.GetAttribs(particle_comps["prev_y"]).dataPtr() + offset; -#else - amrex::ignore_unused(y_old); #endif z_old = pti.GetAttribs(particle_comps["prev_z"]).dataPtr() + offset; + amrex::ignore_unused(x_old, y_old); } // Loop over the particles and update their momentum @@ -3070,7 +2780,7 @@ PhysicalParticleContainer::ImplicitPushXP (WarpXParIter& pti, const Dim3 lo = lbound(box); - const int depos_type = WarpX::current_deposition_algo; + const auto depos_type = WarpX::current_deposition_algo; const int nox = WarpX::nox; const int n_rz_azimuthal_modes = WarpX::n_rz_azimuthal_modes; @@ -3480,6 +3190,28 @@ void PhysicalParticleContainer::resample (const int timestep, const bool verbose WARPX_PROFILE_VAR_STOP(blp_resample_actual); } +bool +PhysicalParticleContainer::findRefinedInjectionBox (amrex::Box& a_fine_injection_box, amrex::IntVect& a_rrfac) +{ + WARPX_PROFILE("PhysicalParticleContainer::findRefinedInjectionBox"); + + // This does not work if the mesh is dynamic. But in that case, we should + // not use refined injected either. We also assume there is only one fine level. + static bool refine_injection = false; + static Box fine_injection_box; + static amrex::IntVect rrfac(AMREX_D_DECL(1,1,1)); + if (!refine_injection and WarpX::moving_window_active(WarpX::GetInstance().getistep(0)+1) and WarpX::refine_plasma and do_continuous_injection and numLevels() == 2) { + refine_injection = true; + fine_injection_box = ParticleBoxArray(1).minimalBox(); + fine_injection_box.setSmall(WarpX::moving_window_dir, std::numeric_limits::lowest()/2); + fine_injection_box.setBig(WarpX::moving_window_dir, std::numeric_limits::max()/2); + rrfac = m_gdb->refRatio(0); + fine_injection_box.coarsen(rrfac); + } + a_fine_injection_box = fine_injection_box; + a_rrfac = rrfac; + return refine_injection; +} #ifdef WARPX_QED diff --git a/Source/Particles/Pusher/GetAndSetPosition.H b/Source/Particles/Pusher/GetAndSetPosition.H index b9e7dc91684..ab06fe3d6cd 100644 --- a/Source/Particles/Pusher/GetAndSetPosition.H +++ b/Source/Particles/Pusher/GetAndSetPosition.H @@ -63,25 +63,16 @@ struct GetParticlePosition { using RType = amrex::ParticleReal; -#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) - const RType* AMREX_RESTRICT m_x = nullptr; - const RType* AMREX_RESTRICT m_z = nullptr; -#elif defined(WARPX_DIM_3D) const RType* AMREX_RESTRICT m_x = nullptr; const RType* AMREX_RESTRICT m_y = nullptr; const RType* AMREX_RESTRICT m_z = nullptr; -#elif defined(WARPX_DIM_1D_Z) - const RType* AMREX_RESTRICT m_z = nullptr; -#endif #if defined(WARPX_DIM_RZ) const RType* m_theta = nullptr; -#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ) - static constexpr RType m_y_default = RType(0.0); -#elif defined(WARPX_DIM_1D_Z) +#endif + static constexpr RType m_x_default = RType(0.0); static constexpr RType m_y_default = RType(0.0); -#endif GetParticlePosition () = default; @@ -118,20 +109,20 @@ struct GetParticlePosition AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator() (const long i, RType& x, RType& y, RType& z) const noexcept { -#ifdef WARPX_DIM_RZ +#if defined(WARPX_DIM_RZ) RType const r = m_x[i]; x = r*std::cos(m_theta[i]); y = r*std::sin(m_theta[i]); z = m_z[i]; -#elif WARPX_DIM_3D +#elif defined(WARPX_DIM_3D) x = m_x[i]; y = m_y[i]; z = m_z[i]; -#elif WARPX_DIM_XZ +#elif defined(WARPX_DIM_XZ) x = m_x[i]; y = m_y_default; z = m_z[i]; -#else +#elif defined(WARPX_DIM_1D_Z) x = m_x_default; y = m_y_default; z = m_z[i]; @@ -146,19 +137,19 @@ struct GetParticlePosition AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AsStored (const long i, RType& x, RType& y, RType& z) const noexcept { -#ifdef WARPX_DIM_RZ +#if defined(WARPX_DIM_RZ) x = m_x[i]; y = m_theta[i]; z = m_z[i]; -#elif WARPX_DIM_3D +#elif defined(WARPX_DIM_3D) x = m_x[i]; y = m_y[i]; z = m_z[i]; -#elif WARPX_DIM_XZ +#elif defined(WARPX_DIM_XZ) x = m_x[i]; y = m_y_default; z = m_z[i]; -#else +#elif defined(WARPX_DIM_1D_Z) x = m_x_default; y = m_y_default; z = m_z[i]; @@ -178,16 +169,17 @@ struct SetParticlePosition { using RType = amrex::ParticleReal; -#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) +#if defined(WARPX_DIM_3D) RType* AMREX_RESTRICT m_x; + RType* AMREX_RESTRICT m_y; RType* AMREX_RESTRICT m_z; -#elif defined(WARPX_DIM_3D) +#elif defined(WARPX_DIM_RZ) || defined(WARPX_DIM_XZ) RType* AMREX_RESTRICT m_x; - RType* AMREX_RESTRICT m_y; RType* AMREX_RESTRICT m_z; #elif defined(WARPX_DIM_1D_Z) RType* AMREX_RESTRICT m_z; #endif + #if defined(WARPX_DIM_RZ) RType* AMREX_RESTRICT m_theta; #endif @@ -217,18 +209,18 @@ struct SetParticlePosition void operator() (const long i, RType x, RType y, RType z) const noexcept { amrex::ignore_unused(x,y,z); -#ifdef WARPX_DIM_RZ +#if defined(WARPX_DIM_RZ) m_theta[i] = std::atan2(y, x); m_x[i] = std::sqrt(x*x + y*y); m_z[i] = z; -#elif WARPX_DIM_3D +#elif defined(WARPX_DIM_3D) m_x[i] = x; m_y[i] = y; m_z[i] = z; -#elif WARPX_DIM_XZ +#elif defined(WARPX_DIM_XZ) m_x[i] = x; m_z[i] = z; -#else +#elif defined(WARPX_DIM_1D_Z) m_z[i] = z; #endif } @@ -241,18 +233,18 @@ struct SetParticlePosition void AsStored (const long i, RType x, RType y, RType z) const noexcept { amrex::ignore_unused(x,y,z); -#ifdef WARPX_DIM_RZ +#if defined(WARPX_DIM_RZ) m_x[i] = x; m_theta[i] = y; m_z[i] = z; -#elif WARPX_DIM_3D +#elif defined(WARPX_DIM_3D) m_x[i] = x; m_y[i] = y; m_z[i] = z; -#elif WARPX_DIM_XZ +#elif defined(WARPX_DIM_XZ) m_x[i] = x; m_z[i] = z; -#else +#elif defined(WARPX_DIM_1D_Z) m_z[i] = z; #endif } diff --git a/Source/Particles/Pusher/PushSelector.H b/Source/Particles/Pusher/PushSelector.H index 4a82e582bfb..d256a1a5e40 100644 --- a/Source/Particles/Pusher/PushSelector.H +++ b/Source/Particles/Pusher/PushSelector.H @@ -49,7 +49,7 @@ void doParticleMomentumPush(amrex::ParticleReal& ux, const int ion_lev, const amrex::ParticleReal m, const amrex::ParticleReal a_q, - const int pusher_algo, + const ParticlePusherAlgo pusher_algo, const int do_crr, #ifdef WARPX_QED const amrex::Real t_chi_max, diff --git a/Source/Particles/Pusher/UpdatePosition.H b/Source/Particles/Pusher/UpdatePosition.H index 89c2de88e47..d11ba6fe21f 100644 --- a/Source/Particles/Pusher/UpdatePosition.H +++ b/Source/Particles/Pusher/UpdatePosition.H @@ -22,7 +22,9 @@ * x^{n+1} - x^{n} = dt*u^{n+1/2}/gamma^{n+1/2} */ AMREX_GPU_HOST_DEVICE AMREX_INLINE -void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, +void UpdatePosition([[maybe_unused]] amrex::ParticleReal& x, + [[maybe_unused]] amrex::ParticleReal& y, + [[maybe_unused]] amrex::ParticleReal& z, const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::Real dt ) { @@ -35,13 +37,9 @@ void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::Parti // Update positions over one time step #if (AMREX_SPACEDIM >= 2) x += ux * inv_gamma * dt; -#else - amrex::ignore_unused(x); #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) // RZ pushes particles in 3D y += uy * inv_gamma * dt; -#else - amrex::ignore_unused(y); #endif z += uz * inv_gamma * dt; } @@ -53,10 +51,12 @@ void UpdatePosition(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::Parti * See Eqs. 15 and 17 in Chen, JCP 407 (2020) 109228. */ AMREX_GPU_HOST_DEVICE AMREX_INLINE -void UpdatePositionImplicit(amrex::ParticleReal& x, amrex::ParticleReal& y, amrex::ParticleReal& z, - const amrex::ParticleReal ux_n, const amrex::ParticleReal uy_n, const amrex::ParticleReal uz_n, - const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, - const amrex::Real dt ) +void UpdatePositionImplicit ([[maybe_unused]] amrex::ParticleReal& x, + [[maybe_unused]] amrex::ParticleReal& y, + [[maybe_unused]] amrex::ParticleReal& z, + const amrex::ParticleReal ux_n, const amrex::ParticleReal uy_n, const amrex::ParticleReal uz_n, + const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, + const amrex::Real dt ) { using namespace amrex::literals; @@ -74,13 +74,9 @@ void UpdatePositionImplicit(amrex::ParticleReal& x, amrex::ParticleReal& y, amre // Update positions over one time step #if (AMREX_SPACEDIM >= 2) x += ux * inv_gamma * dt; -#else - amrex::ignore_unused(x); #endif #if defined(WARPX_DIM_3D) || defined(WARPX_DIM_RZ) // RZ pushes particles in 3D y += uy * inv_gamma * dt; -#else - amrex::ignore_unused(y); #endif z += uz * inv_gamma * dt; } @@ -90,20 +86,19 @@ void UpdatePositionImplicit(amrex::ParticleReal& x, amrex::ParticleReal& y, amre * of the particles for given electric and magnetic fields on the grid. */ AMREX_GPU_HOST_DEVICE AMREX_INLINE -void PositionNorm( amrex::ParticleReal dxp, amrex::ParticleReal dyp, amrex::ParticleReal dzp, - amrex::ParticleReal& dxp_save, amrex::ParticleReal& dyp_save, amrex::ParticleReal& dzp_save, - amrex::ParticleReal idxg2, amrex::ParticleReal idyg2, amrex::ParticleReal idzg2, +void PositionNorm ([[maybe_unused]] amrex::ParticleReal dxp, + [[maybe_unused]] amrex::ParticleReal dyp, + [[maybe_unused]] amrex::ParticleReal dzp, + [[maybe_unused]] amrex::ParticleReal& dxp_save, + [[maybe_unused]] amrex::ParticleReal& dyp_save, + [[maybe_unused]] amrex::ParticleReal& dzp_save, + [[maybe_unused]] amrex::ParticleReal idxg2, + [[maybe_unused]] amrex::ParticleReal idyg2, + [[maybe_unused]] amrex::ParticleReal idzg2, amrex::ParticleReal& step_norm, const int iter ) { using namespace amrex::literals; -#if defined(WARPX_DIM_1D_Z) - amrex::ignore_unused(dxp, dxp_save, idxg2); -#endif -#if !defined(WARPX_DIM_3D) - amrex::ignore_unused(dyp, dyp_save, idyg2); -#endif - if (iter==0) { step_norm = 1.0_prt; } else { step_norm = (dzp - dzp_save)*(dzp - dzp_save)*idzg2; diff --git a/Source/Utils/CMakeLists.txt b/Source/Utils/CMakeLists.txt index 3d804fe9cde..82053bfc88a 100644 --- a/Source/Utils/CMakeLists.txt +++ b/Source/Utils/CMakeLists.txt @@ -6,7 +6,6 @@ foreach(D IN LISTS WarpX_DIMS) ParticleUtils.cpp SpeciesUtils.cpp RelativeCellPosition.cpp - WarpXAlgorithmSelection.cpp WarpXMovingWindow.cpp WarpXTagging.cpp WarpXUtil.cpp diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index dd7e61ff4fa..dc1f1da5c4c 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -2,7 +2,6 @@ CEXE_sources += WarpXMovingWindow.cpp CEXE_sources += WarpXTagging.cpp CEXE_sources += WarpXUtil.cpp CEXE_sources += WarpXVersion.cpp -CEXE_sources += WarpXAlgorithmSelection.cpp CEXE_sources += Interpolate.cpp CEXE_sources += IntervalsParser.cpp CEXE_sources += RelativeCellPosition.cpp diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H index b9105557462..f67aeddadd0 100644 --- a/Source/Utils/WarpXAlgorithmSelection.H +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -9,6 +9,7 @@ #define WARPX_UTILS_WARPXALGORITHMSELECTION_H_ #include +#include #include #include @@ -20,23 +21,19 @@ using namespace ablastr::utils::enums; // NOLINT(google-global-names-in-headers * \brief struct to determine the computational medium, i.e., vacuum or material/macroscopic default is vacuum. */ -struct MediumForEM { - enum { - Vacuum = 0, - Macroscopic = 1 - }; -}; +AMREX_ENUM(MediumForEM, + Vacuum, + Macroscopic, + Default = Vacuum); /** * \brief struct to select the overall evolve scheme */ -struct EvolveScheme { - enum { - Explicit = 0, - ThetaImplicitEM = 1, - SemiImplicitEM = 2 - }; -}; +AMREX_ENUM(EvolveScheme, + Explicit, + ThetaImplicitEM, + SemiImplicitEM, + Default = Explicit); /** * \brief struct to select algorithm for macroscopic Maxwell solver @@ -44,158 +41,116 @@ struct EvolveScheme { Backward Euler (fully-implicit) represents sigma*E = sigma*E^(n+1) default is Backward Euler as it is more robust. */ -struct MacroscopicSolverAlgo { - enum { - BackwardEuler = 0, - LaxWendroff = 1 - }; -}; - -struct ElectromagneticSolverAlgo { - enum { - None = 0, - Yee = 1, - CKC = 2, - PSATD = 3, - ECT = 4, - HybridPIC = 5 - }; -}; - -struct ElectrostaticSolverAlgo { - enum { - None = 0, - Relativistic = 1, - LabFrameElectroMagnetostatic = 2, - LabFrame = 3 // Non relativistic - }; -}; - -struct PoissonSolverAlgo { - enum { - Multigrid = 1, - IntegratedGreenFunction = 2, - }; -}; - -struct ParticlePusherAlgo { - enum { - Boris = 0, - Vay = 1, - HigueraCary = 2 - }; -}; - -struct CurrentDepositionAlgo { - enum { - Esirkepov = 0, - Direct = 1, - Vay = 2, - Villasenor = 3 - }; -}; - -struct ChargeDepositionAlgo { - // Only the Standard algorithm is implemented - enum { - Standard = 0 - }; -}; - -struct GatheringAlgo { - enum { - EnergyConserving = 0, - MomentumConserving - }; -}; - -struct PSATDSolutionType { - enum { - FirstOrder = 0, - SecondOrder = 1 - }; -}; - -struct JInTime { - enum { - Constant = 0, - Linear = 1 - }; -}; - -struct RhoInTime { - enum { - Constant = 0, - Linear = 1 - }; -}; +AMREX_ENUM(MacroscopicSolverAlgo, + BackwardEuler, + LaxWendroff, + Default = BackwardEuler); + +AMREX_ENUM(ElectromagneticSolverAlgo, + None, + Yee, + CKC, + PSATD, + ECT, + HybridPIC, + hybrid = HybridPIC, + Default = Yee); + +AMREX_ENUM(ElectrostaticSolverAlgo, + None, + Relativistic, + LabFrameElectroMagnetostatic, + LabFrame, + Default = None); + +AMREX_ENUM(PoissonSolverAlgo, + Multigrid, + IntegratedGreenFunction, + fft = IntegratedGreenFunction, + Default = Multigrid); + +AMREX_ENUM(ParticlePusherAlgo, + Boris, + Vay, + HigueraCary, + higuera = HigueraCary, + Default = Boris); + +AMREX_ENUM(CurrentDepositionAlgo, + Esirkepov, + Direct, + Vay, + Villasenor, + Default = Esirkepov); + +AMREX_ENUM(ChargeDepositionAlgo, + Standard, + Default = Standard); + +AMREX_ENUM(GatheringAlgo, + EnergyConserving, + MomentumConserving, + Default = EnergyConserving); + +AMREX_ENUM(PSATDSolutionType, + FirstOrder, + SecondOrder, + Default = SecondOrder); + +AMREX_ENUM(JInTime, + Constant, + Linear, + Default = Constant); + +AMREX_ENUM(RhoInTime, + Constant, + Linear, + Default = Linear); /** Strategy to compute weights for use in load balance. */ -struct LoadBalanceCostsUpdateAlgo { - enum { - Timers = 0, //!< load balance according to in-code timer-based weights (i.e., with `costs`) - Heuristic = 1 /**< load balance according to weights computed from number of cells - and number of particles per box (i.e., with `costs_heuristic`) */ - }; -}; +AMREX_ENUM(LoadBalanceCostsUpdateAlgo, + Timers, //!< load balance according to in-code timer-based weights (i.e., with `costs`) + Heuristic, /**< load balance according to weights computed from number of cells + and number of particles per box (i.e., with `costs_heuristic`) */ + Default = Timers); /** Field boundary conditions at the domain boundary */ -enum struct FieldBoundaryType { - PML = 0, - Periodic = 1, - PEC = 2, //!< perfect electric conductor (PEC) with E_tangential=0 - PMC = 3, //!< perfect magnetic conductor (PMC) with B_tangential=0 - Damped = 4, // Fields in the guard cells are damped for PSATD - //in the moving window direction - Absorbing_SilverMueller = 5, // Silver-Mueller boundary condition - Neumann = 6, // For electrostatic, the normal E is set to zero - None = 7, // The fields values at the boundary are not updated. This is - // useful for RZ simulations, at r=0. - Open = 8 // Used in the Integrated Green Function Poisson solver - // Note that the solver implicitely assumes open BCs: - // no need to enforce them separately -}; +AMREX_ENUM(FieldBoundaryType, + PML, + Periodic, + PEC, //!< perfect electric conductor (PEC) with E_tangential=0 + PMC, //!< perfect magnetic conductor (PMC) with B_tangential=0 + Damped, // Fields in the guard cells are damped for PSATD + //in the moving window direction + Absorbing_SilverMueller, // Silver-Mueller boundary condition + absorbingsilvermueller = Absorbing_SilverMueller, + Neumann, // For electrostatic, the normal E is set to zero + None, // The fields values at the boundary are not updated. This is + // useful for RZ simulations, at r=0. + Open, // Used in the Integrated Green Function Poisson solver + // Note that the solver implicitely assumes open BCs: + // no need to enforce them separately + Default = PML); /** Particle boundary conditions at the domain boundary */ -enum struct ParticleBoundaryType { - Absorbing = 0, //!< particles crossing domain boundary are removed - Open = 1, //!< particles cross domain boundary leave with damped j - Reflecting = 2, //!< particles are reflected - Periodic = 3, //!< particles are introduced from the periodic boundary - Thermal = 4, - None = 5 //!< For r=0 boundary with RZ simulations -}; +AMREX_ENUM(ParticleBoundaryType, + Absorbing, //!< particles crossing domain boundary are removed + Open, //!< particles cross domain boundary leave with damped j + Reflecting, //!< particles are reflected + Periodic, //!< particles are introduced from the periodic boundary + Thermal, + None, //!< For r=0 boundary with RZ simulations + Default = Absorbing); /** MPI reductions */ -struct ReductionType { - enum { - Maximum = 0, - Minimum = 1, - Sum = 2 - }; -}; - -int -GetAlgorithmInteger(const amrex::ParmParse& pp, const char* pp_search_key ); - -/** Select BC Type for fields, if field=true - * else select BCType for particles. - */ -FieldBoundaryType -GetFieldBCTypeInteger( std::string BCType ); - -/** Select BC Type for particles. - */ -ParticleBoundaryType -GetParticleBCTypeInteger( std::string BCType ); - -/** Find the name associated with a BC type - */ -std::string -GetFieldBCTypeString( FieldBoundaryType fb_type ); +AMREX_ENUM(ReductionType, + Maximum, + Minimum, + Sum, + Integral = Sum); #endif // WARPX_UTILS_WARPXALGORITHMSELECTION_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp deleted file mode 100644 index edcf5991c71..00000000000 --- a/Source/Utils/WarpXAlgorithmSelection.cpp +++ /dev/null @@ -1,271 +0,0 @@ -/* Copyright 2019-2020 Axel Huebl, David Grote, Luca Fedeli - * Remi Lehe, Weiqun Zhang, Yinjian Zhao - * - * - * This file is part of WarpX. - * - * License: BSD-3-Clause-LBNL - */ -#include "WarpX.H" -#include "WarpXAlgorithmSelection.H" -#include "Utils/TextMsg.H" - -#include - -#include -#include - -#include -#include -#include -#include -#include - -// Define dictionary with correspondence between user-input strings, -// and corresponding integer for use inside the code - -const std::map evolve_scheme_to_int = { - {"explicit", EvolveScheme::Explicit }, - {"theta_implicit_em", EvolveScheme::ThetaImplicitEM }, - {"semi_implicit_em", EvolveScheme::SemiImplicitEM }, - {"default", EvolveScheme::Explicit } -}; - -const std::map grid_to_int = { - {"collocated", static_cast(ablastr::utils::enums::GridType::Collocated)}, - {"staggered", static_cast(ablastr::utils::enums::GridType::Staggered)}, - {"hybrid", static_cast(ablastr::utils::enums::GridType::Hybrid)}, - {"default", static_cast(ablastr::utils::enums::GridType::Staggered)} -}; - -const std::map electromagnetic_solver_algo_to_int = { - {"none", ElectromagneticSolverAlgo::None }, - {"yee", ElectromagneticSolverAlgo::Yee }, - {"ckc", ElectromagneticSolverAlgo::CKC }, - {"psatd", ElectromagneticSolverAlgo::PSATD }, - {"ect", ElectromagneticSolverAlgo::ECT }, - {"hybrid", ElectromagneticSolverAlgo::HybridPIC }, - {"default", ElectromagneticSolverAlgo::Yee } -}; - -const std::map electrostatic_solver_algo_to_int = { - {"none", ElectrostaticSolverAlgo::None }, - {"relativistic", ElectrostaticSolverAlgo::Relativistic}, - {"labframe-electromagnetostatic", ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic}, - {"labframe", ElectrostaticSolverAlgo::LabFrame}, - {"default", ElectrostaticSolverAlgo::None } -}; - -const std::map poisson_solver_algo_to_int = { - {"multigrid", PoissonSolverAlgo::Multigrid}, - {"fft", PoissonSolverAlgo::IntegratedGreenFunction}, - {"default", PoissonSolverAlgo::Multigrid } -}; - -const std::map particle_pusher_algo_to_int = { - {"boris", ParticlePusherAlgo::Boris }, - {"vay", ParticlePusherAlgo::Vay }, - {"higuera", ParticlePusherAlgo::HigueraCary }, - {"default", ParticlePusherAlgo::Boris } -}; - -const std::map current_deposition_algo_to_int = { - {"esirkepov", CurrentDepositionAlgo::Esirkepov }, - {"direct", CurrentDepositionAlgo::Direct }, - {"vay", CurrentDepositionAlgo::Vay }, - {"villasenor", CurrentDepositionAlgo::Villasenor }, - {"default", CurrentDepositionAlgo::Esirkepov } // NOTE: overwritten for PSATD and Hybrid-PIC below -}; - -const std::map charge_deposition_algo_to_int = { - {"standard", ChargeDepositionAlgo::Standard }, - {"default", ChargeDepositionAlgo::Standard } -}; - -const std::map gathering_algo_to_int = { - {"energy-conserving", GatheringAlgo::EnergyConserving }, - {"momentum-conserving", GatheringAlgo::MomentumConserving }, - {"default", GatheringAlgo::EnergyConserving } -}; - -const std::map psatd_solution_type_to_int = { - {"first-order", PSATDSolutionType::FirstOrder}, - {"second-order", PSATDSolutionType::SecondOrder}, - {"default", PSATDSolutionType::SecondOrder} -}; - -const std::map J_in_time_to_int = { - {"constant", JInTime::Constant}, - {"linear", JInTime::Linear}, - {"default", JInTime::Constant} -}; - -const std::map rho_in_time_to_int = { - {"constant", RhoInTime::Constant}, - {"linear", RhoInTime::Linear}, - {"default", RhoInTime::Linear} -}; - -const std::map load_balance_costs_update_algo_to_int = { - {"timers", LoadBalanceCostsUpdateAlgo::Timers }, - {"heuristic", LoadBalanceCostsUpdateAlgo::Heuristic }, - {"default", LoadBalanceCostsUpdateAlgo::Timers } -}; - -const std::map MaxwellSolver_medium_algo_to_int = { - {"vacuum", MediumForEM::Vacuum}, - {"macroscopic", MediumForEM::Macroscopic}, - {"default", MediumForEM::Vacuum} -}; - -const std::map MacroscopicSolver_algo_to_int = { - {"backwardeuler", MacroscopicSolverAlgo::BackwardEuler}, - {"laxwendroff", MacroscopicSolverAlgo::LaxWendroff}, - {"default", MacroscopicSolverAlgo::BackwardEuler} -}; - -const std::map FieldBCType_algo_to_enum = { - {"pml", FieldBoundaryType::PML}, - {"periodic", FieldBoundaryType::Periodic}, - {"pec", FieldBoundaryType::PEC}, - {"pmc", FieldBoundaryType::PMC}, - {"damped", FieldBoundaryType::Damped}, - {"absorbing_silver_mueller", FieldBoundaryType::Absorbing_SilverMueller}, - {"neumann", FieldBoundaryType::Neumann}, - {"open", FieldBoundaryType::Open}, - {"none", FieldBoundaryType::None}, - {"default", FieldBoundaryType::PML} -}; - -const std::map ParticleBCType_algo_to_enum = { - {"absorbing", ParticleBoundaryType::Absorbing}, - {"open", ParticleBoundaryType::Open}, - {"reflecting", ParticleBoundaryType::Reflecting}, - {"periodic", ParticleBoundaryType::Periodic}, - {"thermal", ParticleBoundaryType::Thermal}, - {"none", ParticleBoundaryType::None}, - {"default", ParticleBoundaryType::Absorbing} -}; - -const std::map ReductionType_algo_to_int = { - {"maximum", ReductionType::Maximum}, - {"minimum", ReductionType::Minimum}, - {"integral", ReductionType::Sum} -}; - -int -GetAlgorithmInteger (const amrex::ParmParse& pp, const char* pp_search_key ) -{ - // Read user input ; use "default" if it is not found - std::string algo = "default"; - pp.query( pp_search_key, algo ); - // Convert to lower case - std::transform(algo.begin(), algo.end(), algo.begin(), ::tolower); - - // Pick the right dictionary - std::map algo_to_int; - if (0 == std::strcmp(pp_search_key, "evolve_scheme")) { - algo_to_int = evolve_scheme_to_int; - } else if (0 == std::strcmp(pp_search_key, "maxwell_solver")) { - algo_to_int = electromagnetic_solver_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "grid_type")) { - algo_to_int = grid_to_int; - } else if (0 == std::strcmp(pp_search_key, "do_electrostatic")) { - algo_to_int = electrostatic_solver_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "poisson_solver")) { - algo_to_int = poisson_solver_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "particle_pusher")) { - algo_to_int = particle_pusher_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "current_deposition")) { - algo_to_int = current_deposition_algo_to_int; - if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD || - WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC || - WarpX::electrostatic_solver_id != ElectrostaticSolverAlgo::None) { - algo_to_int["default"] = CurrentDepositionAlgo::Direct; - } - } else if (0 == std::strcmp(pp_search_key, "charge_deposition")) { - algo_to_int = charge_deposition_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "field_gathering")) { - algo_to_int = gathering_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "solution_type")) { - algo_to_int = psatd_solution_type_to_int; - } else if (0 == std::strcmp(pp_search_key, "J_in_time")) { - algo_to_int = J_in_time_to_int; - } else if (0 == std::strcmp(pp_search_key, "rho_in_time")) { - algo_to_int = rho_in_time_to_int; - } else if (0 == std::strcmp(pp_search_key, "load_balance_costs_update")) { - algo_to_int = load_balance_costs_update_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "em_solver_medium")) { - algo_to_int = MaxwellSolver_medium_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "macroscopic_sigma_method")) { - algo_to_int = MacroscopicSolver_algo_to_int; - } else if (0 == std::strcmp(pp_search_key, "reduction_type")) { - algo_to_int = ReductionType_algo_to_int; - } else { - std::string const pp_search_string = pp_search_key; - WARPX_ABORT_WITH_MESSAGE("Unknown algorithm type: " + pp_search_string); - } - - // Check if the user-input is a valid key for the dictionary - if (algo_to_int.count(algo) == 0) { - // Not a valid key ; print error message - const std::string pp_search_string = pp_search_key; - std::string error_message = "Invalid string for algo." + pp_search_string - + ": " + algo + ".\nThe valid values are:\n"; - for ( const auto &valid_pair : algo_to_int ) { - if (valid_pair.first != "default"){ - error_message += " - " + valid_pair.first + "\n"; - } - } - WARPX_ABORT_WITH_MESSAGE(error_message); - } - - // If the input is a valid key, return the value - return algo_to_int[algo]; -} - -FieldBoundaryType -GetFieldBCTypeInteger( std::string BCType ){ - std::transform(BCType.begin(), BCType.end(), BCType.begin(), ::tolower); - - if (FieldBCType_algo_to_enum.count(BCType) == 0) { - std::string error_message = "Invalid string for field/particle BC. : " + BCType + "\nThe valid values are : \n"; - for (const auto &valid_pair : FieldBCType_algo_to_enum) { - if (valid_pair.first != "default"){ - error_message += " - " + valid_pair.first + "\n"; - } - } - WARPX_ABORT_WITH_MESSAGE(error_message); - } - // return FieldBCType_algo_to_enum[BCType]; // This operator cannot be used for a const map - return FieldBCType_algo_to_enum.at(BCType); -} - -ParticleBoundaryType -GetParticleBCTypeInteger( std::string BCType ){ - std::transform(BCType.begin(), BCType.end(), BCType.begin(), ::tolower); - - if (ParticleBCType_algo_to_enum.count(BCType) == 0) { - std::string error_message = "Invalid string for particle BC. : " + BCType + "\nThe valid values are : \n"; - for (const auto &valid_pair : ParticleBCType_algo_to_enum) { - if (valid_pair.first != "default"){ - error_message += " - " + valid_pair.first + "\n"; - } - } - WARPX_ABORT_WITH_MESSAGE(error_message); - } - // return ParticleBCType_algo_to_enum[BCType]; // This operator cannot be used for a const map - return ParticleBCType_algo_to_enum.at(BCType); -} - -std::string -GetFieldBCTypeString( FieldBoundaryType fb_type ) { - std::string boundary_name; - for (const auto &valid_pair : FieldBCType_algo_to_enum) { - if ((valid_pair.second == fb_type)&&(valid_pair.first != "default")){ - boundary_name = valid_pair.first; - break; - } - } - return boundary_name; -} diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index 44952091257..856e021abb3 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -77,7 +77,8 @@ void ParseGeometryInput() #ifdef WARPX_DIM_RZ const ParmParse pp_algo("algo"); - const int electromagnetic_solver_id = GetAlgorithmInteger(pp_algo, "maxwell_solver"); + auto electromagnetic_solver_id = ElectromagneticSolverAlgo::Default; + pp_algo.query_enum_sloppy("maxwell_solver", electromagnetic_solver_id, "-_"); if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD) { WARPX_ALWAYS_ASSERT_WITH_MESSAGE(prob_lo[0] == 0., @@ -346,7 +347,8 @@ void CheckGriddingForRZSpectral () CheckDims(); const ParmParse pp_algo("algo"); - const int electromagnetic_solver_id = GetAlgorithmInteger(pp_algo, "maxwell_solver"); + auto electromagnetic_solver_id = ElectromagneticSolverAlgo::Default; + pp_algo.query_enum_sloppy("maxwell_solver", electromagnetic_solver_id, "-_"); // only check for PSATD in RZ if (electromagnetic_solver_id != ElectromagneticSolverAlgo::PSATD) { @@ -431,16 +433,14 @@ void CheckGriddingForRZSpectral () void ReadBCParams () { - amrex::Vector field_BC_lo(AMREX_SPACEDIM,"default"); - amrex::Vector field_BC_hi(AMREX_SPACEDIM,"default"); - amrex::Vector particle_BC_lo(AMREX_SPACEDIM,"default"); - amrex::Vector particle_BC_hi(AMREX_SPACEDIM,"default"); amrex::Vector geom_periodicity(AMREX_SPACEDIM,0); ParmParse pp_geometry("geometry"); const ParmParse pp_warpx("warpx"); const ParmParse pp_algo("algo"); - const int electromagnetic_solver_id = GetAlgorithmInteger(pp_algo, "maxwell_solver"); - const int poisson_solver_id = GetAlgorithmInteger(pp_warpx, "poisson_solver"); + auto electromagnetic_solver_id = ElectromagneticSolverAlgo::Default; + pp_algo.query_enum_sloppy("maxwell_solver", electromagnetic_solver_id, "-_"); + auto poisson_solver_id = PoissonSolverAlgo::Default; + pp_warpx.query_enum_sloppy("poisson_solver", poisson_solver_id, "-_"); if (pp_geometry.queryarr("is_periodic", geom_periodicity)) { @@ -455,26 +455,21 @@ void ReadBCParams () // particle boundary may not be explicitly specified for some applications bool particle_boundary_specified = false; const ParmParse pp_boundary("boundary"); - pp_boundary.queryarr("field_lo", field_BC_lo, 0, AMREX_SPACEDIM); - pp_boundary.queryarr("field_hi", field_BC_hi, 0, AMREX_SPACEDIM); - if (pp_boundary.queryarr("particle_lo", particle_BC_lo, 0, AMREX_SPACEDIM)) { - particle_boundary_specified = true; - } - if (pp_boundary.queryarr("particle_hi", particle_BC_hi, 0, AMREX_SPACEDIM)) { - particle_boundary_specified = true; - } - AMREX_ALWAYS_ASSERT(field_BC_lo.size() == AMREX_SPACEDIM); - AMREX_ALWAYS_ASSERT(field_BC_hi.size() == AMREX_SPACEDIM); - AMREX_ALWAYS_ASSERT(particle_BC_lo.size() == AMREX_SPACEDIM); - AMREX_ALWAYS_ASSERT(particle_BC_hi.size() == AMREX_SPACEDIM); - for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { // Get field boundary type - WarpX::field_boundary_lo[idim] = GetFieldBCTypeInteger(field_BC_lo[idim]); - WarpX::field_boundary_hi[idim] = GetFieldBCTypeInteger(field_BC_hi[idim]); + pp_boundary.query_enum_sloppy("field_lo", + WarpX::field_boundary_lo[idim], "-_", idim); + pp_boundary.query_enum_sloppy("field_hi", + WarpX::field_boundary_hi[idim], "-_", idim); // Get particle boundary type - WarpX::particle_boundary_lo[idim] = GetParticleBCTypeInteger(particle_BC_lo[idim]); - WarpX::particle_boundary_hi[idim] = GetParticleBCTypeInteger(particle_BC_hi[idim]); + if (pp_boundary.query_enum_sloppy("particle_lo", + WarpX::particle_boundary_lo[idim], "-_", idim)) { + particle_boundary_specified = true; + } + if (pp_boundary.query_enum_sloppy("particle_hi", + WarpX::particle_boundary_hi[idim], "-_", idim)) { + particle_boundary_specified = true; + } if (WarpX::field_boundary_lo[idim] == FieldBoundaryType::Periodic || WarpX::field_boundary_hi[idim] == FieldBoundaryType::Periodic || diff --git a/Source/WarpX.H b/Source/WarpX.H index 1fc7bc2b1fb..581ef20eb42 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -172,17 +172,17 @@ public: // Algorithms //! Integer that corresponds to the current deposition algorithm (Esirkepov, direct, Vay, Villasenor) - static short current_deposition_algo; + static inline auto current_deposition_algo = CurrentDepositionAlgo::Default; //! Integer that corresponds to the charge deposition algorithm (only standard deposition) - static short charge_deposition_algo; + static inline auto charge_deposition_algo = ChargeDepositionAlgo::Default; //! Integer that corresponds to the field gathering algorithm (energy-conserving, momentum-conserving) - static short field_gathering_algo; + static inline auto field_gathering_algo = GatheringAlgo::Default; //! Integer that corresponds to the particle push algorithm (Boris, Vay, Higuera-Cary) - static short particle_pusher_algo; + static inline auto particle_pusher_algo = ParticlePusherAlgo::Default; //! Integer that corresponds to the type of Maxwell solver (Yee, CKC, PSATD, ECT) - static short electromagnetic_solver_id; + static inline auto electromagnetic_solver_id = ElectromagneticSolverAlgo::Default; //! Integer that corresponds to the evolve scheme (explicit, semi_implicit_em, theta_implicit_em) - static short evolve_scheme; + static inline auto evolve_scheme = EvolveScheme::Default; //! Maximum iterations used for self-consistent particle update in implicit particle-suppressed evolve schemes static int max_particle_its_in_implicit_scheme; //! Relative tolerance used for self-consistent particle update in implicit particle-suppressed evolve schemes @@ -190,43 +190,55 @@ public: /** Records a number corresponding to the load balance cost update strategy * being used (0 or 1 corresponding to timers or heuristic). */ - static short load_balance_costs_update_algo; + static inline auto load_balance_costs_update_algo = LoadBalanceCostsUpdateAlgo::Default; //! Integer that corresponds to electromagnetic Maxwell solver (vacuum - 0, macroscopic - 1) - static int em_solver_medium; + static inline auto em_solver_medium = MediumForEM::Default; /** Integer that correspond to macroscopic Maxwell solver algorithm * (BackwardEuler - 0, Lax-Wendroff - 1) */ - static int macroscopic_solver_algo; + static inline auto macroscopic_solver_algo = MacroscopicSolverAlgo::Default; /** Integers that correspond to boundary condition applied to fields at the * lower domain boundaries * (0 to 6 correspond to PML, Periodic, PEC, PMC, Damped, Absorbing Silver-Mueller, None) */ - static amrex::Vector field_boundary_lo; + static inline amrex::Array + field_boundary_lo {AMREX_D_DECL(FieldBoundaryType::Default, + FieldBoundaryType::Default, + FieldBoundaryType::Default)}; /** Integers that correspond to boundary condition applied to fields at the * upper domain boundaries * (0 to 6 correspond to PML, Periodic, PEC, PMC, Damped, Absorbing Silver-Mueller, None) */ - static amrex::Vector field_boundary_hi; + static inline amrex::Array + field_boundary_hi {AMREX_D_DECL(FieldBoundaryType::Default, + FieldBoundaryType::Default, + FieldBoundaryType::Default)}; /** Integers that correspond to boundary condition applied to particles at the * lower domain boundaries * (0 to 4 correspond to Absorbing, Open, Reflecting, Periodic, Thermal) */ - static amrex::Vector particle_boundary_lo; + static inline amrex::Array + particle_boundary_lo {AMREX_D_DECL(ParticleBoundaryType::Default, + ParticleBoundaryType::Default, + ParticleBoundaryType::Default)}; /** Integers that correspond to boundary condition applied to particles at the * upper domain boundaries * (0 to 4 correspond to Absorbing, Open, Reflecting, Periodic, Thermal) */ - static amrex::Vector particle_boundary_hi; + static inline amrex::Array + particle_boundary_hi {AMREX_D_DECL(ParticleBoundaryType::Default, + ParticleBoundaryType::Default, + ParticleBoundaryType::Default)}; //! Integer that corresponds to the order of the PSATD solution //! (whether the PSATD equations are derived from first-order or //! second-order solution) - static short psatd_solution_type; + static inline auto psatd_solution_type = PSATDSolutionType::Default; //! Integers that correspond to the time dependency of J (constant, linear) //! and rho (linear, quadratic) for the PSATD algorithm - static short J_in_time; - static short rho_in_time; + static inline auto J_in_time = JInTime::Default; + static inline auto rho_in_time = RhoInTime::Default; //! If true, the current is deposited on a nodal grid and then centered onto a staggered grid //! using finite centering of order given by #current_centering_nox, #current_centering_noy, @@ -379,7 +391,7 @@ public: //! Integer that corresponds to the type of grid used in the simulation //! (collocated, staggered, hybrid) - static ablastr::utils::enums::GridType grid_type; + static inline auto grid_type = ablastr::utils::enums::GridType::Default; // Global rho nodal flag to know about rho index type when rho MultiFab is not allocated amrex::IntVect m_rho_nodal_flag; @@ -858,8 +870,8 @@ public: static const amrex::iMultiFab* CurrentBufferMasks (int lev); static const amrex::iMultiFab* GatherBufferMasks (int lev); - static int electrostatic_solver_id; - static int poisson_solver_id; + static inline auto electrostatic_solver_id = ElectrostaticSolverAlgo::Default; + static inline auto poisson_solver_id = PoissonSolverAlgo::Default; // Parameters for lab frame electrostatic static amrex::Real self_fields_required_precision; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index f62dcde1e0d..a30381b7c5f 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -110,23 +110,11 @@ bool WarpX::compute_max_step_from_btd = false; Real WarpX::zmax_plasma_to_compute_max_step = 0._rt; Real WarpX::zmin_domain_boost_step_0 = 0._rt; -short WarpX::current_deposition_algo; -short WarpX::charge_deposition_algo; -short WarpX::field_gathering_algo; -short WarpX::particle_pusher_algo; -short WarpX::electromagnetic_solver_id; -short WarpX::evolve_scheme; int WarpX::max_particle_its_in_implicit_scheme = 21; ParticleReal WarpX::particle_tol_in_implicit_scheme = 1.e-10; -short WarpX::psatd_solution_type; -short WarpX::J_in_time; -short WarpX::rho_in_time; -short WarpX::load_balance_costs_update_algo; bool WarpX::do_dive_cleaning = false; bool WarpX::do_divb_cleaning = false; bool WarpX::do_divb_cleaning_external = false; -int WarpX::em_solver_medium; -int WarpX::macroscopic_solver_algo; bool WarpX::do_single_precision_comms = false; bool WarpX::do_shared_mem_charge_deposition = false; @@ -141,11 +129,6 @@ amrex::IntVect WarpX::shared_tilesize(AMREX_D_DECL(1,1,1)); #endif int WarpX::shared_mem_current_tpb = 128; -amrex::Vector WarpX::field_boundary_lo(AMREX_SPACEDIM,FieldBoundaryType::PML); -amrex::Vector WarpX::field_boundary_hi(AMREX_SPACEDIM,FieldBoundaryType::PML); -amrex::Vector WarpX::particle_boundary_lo(AMREX_SPACEDIM,ParticleBoundaryType::Absorbing); -amrex::Vector WarpX::particle_boundary_hi(AMREX_SPACEDIM,ParticleBoundaryType::Absorbing); - int WarpX::n_rz_azimuthal_modes = 1; int WarpX::ncomps = 1; @@ -191,8 +174,6 @@ amrex::IntVect WarpX::sort_idx_type(AMREX_D_DECL(0,0,0)); bool WarpX::do_dynamic_scheduling = true; -int WarpX::electrostatic_solver_id; -int WarpX::poisson_solver_id; Real WarpX::self_fields_required_precision = 1.e-11_rt; Real WarpX::self_fields_absolute_tolerance = 0.0_rt; int WarpX::self_fields_max_iters = 200; @@ -211,7 +192,6 @@ IntVect WarpX::filter_npass_each_dir(1); int WarpX::n_field_gather_buffer = -1; int WarpX::n_current_deposition_buffer = -1; -ablastr::utils::enums::GridType WarpX::grid_type; amrex::IntVect m_rho_nodal_flag; WarpX* WarpX::m_instance = nullptr; @@ -473,8 +453,7 @@ WarpX::ReadParameters () { const ParmParse pp_algo("algo"); - electromagnetic_solver_id = static_cast(GetAlgorithmInteger(pp_algo, "maxwell_solver")); - + pp_algo.query_enum_sloppy("maxwell_solver", electromagnetic_solver_id, "-_"); if (electromagnetic_solver_id == ElectromagneticSolverAlgo::ECT && !EB::enabled()) { throw std::runtime_error("ECP Solver requires to enable embedded boundaries at runtime."); } @@ -681,7 +660,7 @@ WarpX::ReadParameters () maxlevel_extEMfield_init = maxLevel(); pp_warpx.query("maxlevel_extEMfield_init", maxlevel_extEMfield_init); - electrostatic_solver_id = GetAlgorithmInteger(pp_warpx, "do_electrostatic"); + pp_warpx.query_enum_sloppy("do_electrostatic", electrostatic_solver_id, "-_"); // if an electrostatic solver is used, set the Maxwell solver to None if (electrostatic_solver_id != ElectrostaticSolverAlgo::None) { electromagnetic_solver_id = ElectromagneticSolverAlgo::None; @@ -701,7 +680,7 @@ WarpX::ReadParameters () pp_warpx.query("self_fields_verbosity", self_fields_verbosity); } - poisson_solver_id = GetAlgorithmInteger(pp_warpx, "poisson_solver"); + pp_warpx.query_enum_sloppy("poisson_solver", poisson_solver_id, "-_"); #ifndef WARPX_DIM_3D WARPX_ALWAYS_ASSERT_WITH_MESSAGE( poisson_solver_id!=PoissonSolverAlgo::IntegratedGreenFunction, @@ -1022,7 +1001,7 @@ WarpX::ReadParameters () // Integer that corresponds to the type of grid used in the simulation // (collocated, staggered, hybrid) - grid_type = static_cast(GetAlgorithmInteger(pp_warpx, "grid_type")); + pp_warpx.query_enum_sloppy("grid_type", grid_type, "-_"); // Use same shape factors in all directions, for gathering if (grid_type == GridType::Collocated) { galerkin_interpolation = false; } @@ -1165,10 +1144,15 @@ WarpX::ReadParameters () // note: current_deposition must be set after maxwell_solver (electromagnetic_solver_id) or // do_electrostatic (electrostatic_solver_id) are already determined, // because its default depends on the solver selection - current_deposition_algo = static_cast(GetAlgorithmInteger(pp_algo, "current_deposition")); - charge_deposition_algo = static_cast(GetAlgorithmInteger(pp_algo, "charge_deposition")); - particle_pusher_algo = static_cast(GetAlgorithmInteger(pp_algo, "particle_pusher")); - evolve_scheme = static_cast(GetAlgorithmInteger(pp_algo, "evolve_scheme")); + if (electromagnetic_solver_id == ElectromagneticSolverAlgo::PSATD || + electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC || + electrostatic_solver_id != ElectrostaticSolverAlgo::None) { + current_deposition_algo = CurrentDepositionAlgo::Direct; + } + pp_algo.query_enum_sloppy("current_deposition", current_deposition_algo, "-_"); + pp_algo.query_enum_sloppy("charge_deposition", charge_deposition_algo, "-_"); + pp_algo.query_enum_sloppy("particle_pusher", particle_pusher_algo, "-_"); + pp_algo.query_enum_sloppy("evolve_scheme", evolve_scheme, "-_"); // check for implicit evolve scheme if (evolve_scheme == EvolveScheme::SemiImplicitEM) { @@ -1229,7 +1213,7 @@ WarpX::ReadParameters () // Query algo.field_gathering from input, set field_gathering_algo to // "default" if not found (default defined in Utils/WarpXAlgorithmSelection.cpp) - field_gathering_algo = static_cast(GetAlgorithmInteger(pp_algo, "field_gathering")); + pp_algo.query_enum_sloppy("field_gathering", field_gathering_algo, "-_"); // Set default field gathering algorithm for hybrid grids (momentum-conserving) std::string tmp_algo; @@ -1286,9 +1270,10 @@ WarpX::ReadParameters () " combined with mesh refinement is currently not implemented"); } - em_solver_medium = GetAlgorithmInteger(pp_algo, "em_solver_medium"); + pp_algo.query_enum_sloppy("em_solver_medium", em_solver_medium, "-_"); if (em_solver_medium == MediumForEM::Macroscopic ) { - macroscopic_solver_algo = GetAlgorithmInteger(pp_algo,"macroscopic_sigma_method"); + pp_algo.query_enum_sloppy("macroscopic_sigma_method", + macroscopic_solver_algo, "-_"); } if (evolve_scheme == EvolveScheme::SemiImplicitEM || @@ -1327,7 +1312,7 @@ WarpX::ReadParameters () } utils::parser::queryWithParser(pp_algo, "load_balance_efficiency_ratio_threshold", load_balance_efficiency_ratio_threshold); - load_balance_costs_update_algo = static_cast(GetAlgorithmInteger(pp_algo, "load_balance_costs_update")); + pp_algo.query_enum_sloppy("load_balance_costs_update", load_balance_costs_update_algo, "-_"); if (WarpX::load_balance_costs_update_algo==LoadBalanceCostsUpdateAlgo::Heuristic) { utils::parser::queryWithParser( pp_algo, "costs_heuristic_cells_wt", costs_heuristic_cells_wt); @@ -1498,12 +1483,12 @@ WarpX::ReadParameters () // Integer that corresponds to the order of the PSATD solution // (whether the PSATD equations are derived from first-order or // second-order solution) - psatd_solution_type = static_cast(GetAlgorithmInteger(pp_psatd, "solution_type")); + pp_psatd.query_enum_sloppy("solution_type", psatd_solution_type, "-_"); // Integers that correspond to the time dependency of J (constant, linear) // and rho (linear, quadratic) for the PSATD algorithm - J_in_time = static_cast(GetAlgorithmInteger(pp_psatd, "J_in_time")); - rho_in_time = static_cast(GetAlgorithmInteger(pp_psatd, "rho_in_time")); + pp_psatd.query_enum_sloppy("J_in_time", J_in_time, "-_"); + pp_psatd.query_enum_sloppy("rho_in_time", rho_in_time, "-_"); if (psatd_solution_type != PSATDSolutionType::FirstOrder || !do_multi_J) { diff --git a/Source/ablastr/utils/Enums.H b/Source/ablastr/utils/Enums.H index 1f89bede9e4..7c7129cae77 100644 --- a/Source/ablastr/utils/Enums.H +++ b/Source/ablastr/utils/Enums.H @@ -8,17 +8,19 @@ #ifndef ABLASTR_UTILS_ENUMS_H_ #define ABLASTR_UTILS_ENUMS_H_ +#include + namespace ablastr::utils::enums { /** Type of grids used in a simulation: * * Collocated at the same location (AMReX: all "NODAL"), staggered (Yee-style), or hybrid. */ - enum struct GridType { - Collocated = 0, - Staggered = 1, - Hybrid = 2 - }; + AMREX_ENUM(GridType, + Collocated, + Staggered, + Hybrid, + Default = Staggered); /** Mesh-refinement patch * diff --git a/Tools/Release/updateAMReX.py b/Tools/Release/updateAMReX.py index 7ba3bca8357..beeb12e85ff 100755 --- a/Tools/Release/updateAMReX.py +++ b/Tools/Release/updateAMReX.py @@ -15,15 +15,6 @@ import requests -try: - from configupdater import ConfigUpdater -except ImportError: - print("Warning: Cannot update .ini files without 'configupdater'") - print("Consider running 'python -m pip install configupdater'") - ConfigUpdater = None - sys.exit(1) - - # Maintainer Inputs ########################################################### print("""Hi there, this is a WarpX maintainer tool to update the source @@ -110,22 +101,6 @@ # Updates ##################################################################### -# run_test.sh (used also for Azure Pipelines) -run_test_path = str(REPO_DIR.joinpath("run_test.sh")) -with open(run_test_path, encoding="utf-8") as f: - run_test_content = f.read() - # branch/commit/tag (git fetcher) version - # cd amrex && git checkout COMMIT_TAG_OR_BRANCH && cd - - run_test_content = re.sub( - r"(.*cd\s+amrex.+git checkout\s+--detach\s+)(.+)(\s+&&\s.*)", - r"\g<1>{}\g<3>".format(amrex_new_branch), - run_test_content, - flags=re.MULTILINE, - ) - -with open(run_test_path, "w", encoding="utf-8") as f: - f.write(run_test_content) - # CI: legacy build check in .github/workflows/cuda.yml ci_gnumake_path = str(REPO_DIR.joinpath(".github/workflows/cuda.yml")) with open(ci_gnumake_path, encoding="utf-8") as f: @@ -142,23 +117,6 @@ with open(ci_gnumake_path, "w", encoding="utf-8") as f: f.write(ci_gnumake_content) -if ConfigUpdater is not None: - # WarpX-tests.ini - tests_ini_path = str(REPO_DIR.joinpath("Regression/WarpX-tests.ini")) - cp = ConfigUpdater() - cp.optionxform = str - cp.read(tests_ini_path) - cp["AMReX"]["branch"].value = amrex_new_branch - cp.update_file() - - # WarpX-GPU-tests.ini - tests_gpu_ini_path = str(REPO_DIR.joinpath("Regression/WarpX-GPU-tests.ini")) - cp = ConfigUpdater() - cp.optionxform = str - cp.read(tests_gpu_ini_path) - cp["AMReX"]["branch"].value = amrex_new_branch - cp.update_file() - # WarpX references to AMReX: cmake/dependencies/AMReX.cmake with open(amrex_cmake_path, encoding="utf-8") as f: amrex_cmake_content = f.read() diff --git a/cmake/dependencies/AMReX.cmake b/cmake/dependencies/AMReX.cmake index f65f5d36cce..e3682b69ff5 100644 --- a/cmake/dependencies/AMReX.cmake +++ b/cmake/dependencies/AMReX.cmake @@ -279,7 +279,7 @@ set(WarpX_amrex_src "" set(WarpX_amrex_repo "https://github.com/AMReX-Codes/amrex.git" CACHE STRING "Repository URI to pull and build AMReX from if(WarpX_amrex_internal)") -set(WarpX_amrex_branch "216ce6f37de4b65be57fc1006b3457b4fc318e03" +set(WarpX_amrex_branch "028638564f7be0694b9898f8d4088cdbf9a6f9f5" CACHE STRING "Repository branch for WarpX_amrex_repo if(WarpX_amrex_internal)") diff --git a/cmake/dependencies/pyAMReX.cmake b/cmake/dependencies/pyAMReX.cmake index 4c92ffa99ba..e93851443c0 100644 --- a/cmake/dependencies/pyAMReX.cmake +++ b/cmake/dependencies/pyAMReX.cmake @@ -74,7 +74,7 @@ option(WarpX_pyamrex_internal "Download & build pyAMReX" ON) set(WarpX_pyamrex_repo "https://github.com/AMReX-Codes/pyamrex.git" CACHE STRING "Repository URI to pull and build pyamrex from if(WarpX_pyamrex_internal)") -set(WarpX_pyamrex_branch "da2d5a000330395b3fcbcb43a519b3c8a318c584" +set(WarpX_pyamrex_branch "41c856b8a588c3c8b04bb35d2d05b56f6ce0dd7f" CACHE STRING "Repository branch for WarpX_pyamrex_repo if(WarpX_pyamrex_internal)")