diff options
52 files changed, 1659 insertions, 86 deletions
diff --git a/Docs/requirements.txt b/Docs/requirements.txt index 478b0ba60..62755cc38 100644 --- a/Docs/requirements.txt +++ b/Docs/requirements.txt @@ -2,3 +2,5 @@ sphinx_rtd_theme>=0.3.1 recommonmark sphinx pygments +breathe +exhale diff --git a/Docs/source/conf.py b/Docs/source/conf.py index 475504d7c..4a4716b66 100644 --- a/Docs/source/conf.py +++ b/Docs/source/conf.py @@ -17,11 +17,10 @@ # add these directories to sys.path here. If the directory is relative to the # documentation root, use os.path.abspath to make it absolute, like shown here. # -import os, sys +import os, sys, subprocess import sphinx_rtd_theme sys.path.insert(0, os.path.join( os.path.abspath(__file__), '../Python') ) - # -- General configuration ------------------------------------------------ # If your documentation needs a minimal Sphinx version, state it here. @@ -33,7 +32,10 @@ sys.path.insert(0, os.path.join( os.path.abspath(__file__), '../Python') ) # ones. extensions = ['sphinx.ext.autodoc', 'sphinx.ext.mathjax', - 'sphinx.ext.viewcode' ] + 'sphinx.ext.viewcode', + 'breathe', + 'exhale' + ] # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -162,3 +164,34 @@ texinfo_documents = [ # Example configuration for intersphinx: refer to the Python standard library. intersphinx_mapping = {'https://amrex-codes.github.io/': None} + +# Setup the breathe extension +breathe_projects = { + "WarpX": "../doxyxml/" +} +breathe_default_project = "WarpX" + +# Setup the exhale extension +exhale_args = { + # These arguments are required + "containmentFolder": "./api", + "rootFileName": "library_root.rst", + "rootFileTitle": "Doxygen", + "doxygenStripFromPath": "..", + # Suggested optional arguments + "createTreeView": True, + # TIP: if using the sphinx-bootstrap-theme, you need + # "treeViewIsBootstrap": True, + "exhaleExecutesDoxygen": True, + "exhaleDoxygenStdin": "INPUT = ../../Source/" +} + +# Tell sphinx what the primary language being documented is. +primary_domain = 'cpp' +# Tell sphinx what the pygments highlight language should be. +highlight_language = 'cpp' + + +read_the_docs_build = os.environ.get('READTHEDOCS', None) == 'True' +if read_the_docs_build: + subprocess.call('cd ../; doxygen', shell=True) diff --git a/Docs/source/developers/amrex_basics.rst b/Docs/source/developers/amrex_basics.rst new file mode 100644 index 000000000..577a6547b --- /dev/null +++ b/Docs/source/developers/amrex_basics.rst @@ -0,0 +1,32 @@ +.. _developers-amrex-basics: + +AMReX basics (excessively basic) +================================ + +WarpX is built on the Adaptive Mesh Refinement (AMR) library `AMReX <https://github.com/AMReX-Codes/amrex>`__. This section provides a very sporadic description of the main AMReX classes and concepts relevant for WarpX, that can serve as a reminder. Please read the AMReX basics `doc page <https://amrex-codes.github.io/amrex/docs_html/Basics.html>`__, of which this section is largely inspired. + +* ``amrex::Box``: Dimension-dependent lower and upper indices defining a rectangular volume in 3D (or surface in 2D) in the index space. ``Box`` is a lightweight meta-data class, with useful member functions. + +* ``amrex::BoxArray``: Collection of ``Box`` on a single AMR level. The information of which MPI rank owns which ``Box`` in a ``BoxArray`` is in ``DistributionMapping``. + +* ``amrex::FArrayBox``: Fortran-ordered array of floating-point ``amrex::Real`` elements defined on a ``Box``. A ``FArrayBox`` can represent scalar data or vector data, with ``ncomp`` components. + +* ``amrex::MultiFab``: Collection of `FAB` (= ``FArrayBox``) on a single AMR level, distributed over MPI ranks. The concept of `ghost cells` is defined at the ``MultiFab`` level. + +* ``amrex::ParticleContainer``: A collection of particles, typically for particles of a physical species. Particles in a ``ParticleContainer`` are organized per ``Box``. Particles in a ``Box`` are organized per tile (this feature is off when running on GPU). Particles within a tile are stored in several structures, each being contiguous in memory: (i) an Array-Of-Struct (AoS) (often called `data`, they are the 3D position, the particle ID and the index of the CPU owning the particle), where the Struct is an ``amrex::Particle`` and (ii) Struct-Of-Arrays (SoA) for extra variables (often called ``attribs``, in WarpX they are the momentum, field on particle etc.). + +The simulation domain is decomposed in several ``Box``, and each MPI rank owns (and performs operations on) the fields and particles defined on a few of these ``Box``, but has the metadata of all of them. For convenience, AMReX provides iterators, to easily iterate over all ``FArrayBox`` (or even tile-by-tile, optionally) in a ``MultiFab`` own by the MPI rank (``MFIter``), or over all particles in a ``ParticleContainer`` on a per-box basis (``ParIter``, or its derived class ``WarpXParIter``). These are respectively done in loops like: + +.. code-block:: cpp + + // mf is a pointer to MultiFab + for ( amrex::MFIter mfi(mf, false); mfi.isValid(); ++mfi ) { ... } + +and + +.. code-block:: cpp + + // *this is a pointer to a ParticleContainer + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { ... } + +When looping over ``FArrayBox`` in a ``MultiFab``, the iterator provides functions to retrieve the metadata of the ``Box`` on which the ``FAB`` is defined (``MFIter::box()``, ``MFIter::tilebox()`` or variations) or the particles defined on this ``Box`` (``ParIter::GetParticles()``). diff --git a/Docs/source/developers/developers.rst b/Docs/source/developers/developers.rst new file mode 100644 index 000000000..4eab0102a --- /dev/null +++ b/Docs/source/developers/developers.rst @@ -0,0 +1,20 @@ +.. _developers: + +Developers documentation +======================== + +For general information on how to contribute to WarpX, including our ``git`` workflow and code practices, have a look at our `CONTRIBUTING.md <https://github.com/ECP-WarpX/WarpX/blob/dev/CONTRIBUTING.md>`__! + +.. toctree:: + :maxdepth: 1 + + amrex_basics + repo_organization + fields + particles + initialization + diagnostics + moving_window + portability + profiling + documentation diff --git a/Docs/source/developers/diagnostics.rst b/Docs/source/developers/diagnostics.rst new file mode 100644 index 000000000..6945670b5 --- /dev/null +++ b/Docs/source/developers/diagnostics.rst @@ -0,0 +1,18 @@ +.. _developers-diagnostics: + +Diagnostics +=========== + +Regular Diagnostics (plotfiles) +------------------------------- + +.. note:: + + Section empty! + +Back-Transformed Diagnostics +---------------------------- + +.. note:: + + Section empty! diff --git a/Docs/source/developers/documentation.rst b/Docs/source/developers/documentation.rst new file mode 100644 index 000000000..a6f06c7bd --- /dev/null +++ b/Docs/source/developers/documentation.rst @@ -0,0 +1,48 @@ +.. _developers-docs: + +Documentation +============= + +.. warning:: + + Needs info on how to install all pieces to compile the doc! Some basic info is in the ``CONTRIBUTING.md`` . + +Doxygen documentation +--------------------- + +WarpX uses a Doxygen documentation. Whenever you create a new class, please document it where it is declared (typically in the header file): + +.. code-block:: cpp + + /** + * few-line description explaining the purpose of my_class. + * + * If you are kind enough, also quickly explain how things in my_class work. + * (typically a few more lines) + */ + class my_class + { ... } + +Doxygen reads this docstring, so please be accurate with the syntax! See `Doxygen manual <http://www.doxygen.nl/manual/docblocks.html>`__ for more information. Similarly, please document functions when you declare them (typically in a header file) like: + +.. code-block:: cpp + + /** + * few-line description explaining the purpose of my_function. + * + * \param[in,out] my_int a pointer to an integer variable on which + * my_function will operate. + */ + void my_class::my_function(int* my_int); + +Breathe documentation +--------------------- + +Your Doxygen documentation is not only useful for people looking into the code, it is also part of the `WarpX online documentation <https://ecp-warpx.github.io>`__ based on `Sphinx <http://www.sphinx-doc.org/en/master/>`__! This is done using the Python module `Breathe <http://breathe.readthedocs.org>`__, that allows you to read Doxygen documentation dorectly in the source and include it in your Sphinx documentation, by calling Breathe functions. For instance, the following line will get the Doxygen documentation for ``WarpXParticleContainer`` in ``Source/Particles/WarpXParticleContainer.H`` and include it to the html page generated by Sphinx: + + .. doxygenclass:: WarpXParticleContainer + +Exhale documentation +-------------------- + +Very similar to Breathe, the Python module `exhale <https://exhale.readthedocs.io/en/latest/>`__ reads the full Doxygen documentation and renders it in ` rst <https://en.wikipedia.org/wiki/ReStructuredText>`__ format, and is accessible from the main WarpX ReadTheDocs page. diff --git a/Docs/source/developers/fields.rst b/Docs/source/developers/fields.rst new file mode 100644 index 000000000..6f43246b8 --- /dev/null +++ b/Docs/source/developers/fields.rst @@ -0,0 +1,127 @@ +.. _developers-fields: + +Fields +====== + +.. note:: + + Add info on staggering and domain decomposition. Synchronize with section ``initialization``. + +The main fields are the electric field ``Efield``, the magnetic field ``Bfield``, the current density ``current`` and the charge density ``rho``. When a divergence-cleaner is used, we add another field ``F`` (containing :math:`\vec \nabla \cdot \vec E - \rho`). + +Due the AMR strategy used in WarpX (see section :ref:`Theory: AMR <theory-amr>` for a complete description), each field on a given refinement level ``lev`` (except for the coarsest ``0``) is defined on: + +* **the fine patch** (suffix ``_fp``, the actual resolution on ``lev``). + +* **the coarse patch** (suffix ``_cp``, same physical domain with the resolution of MR level ``lev-1``). + +* **the auxiliary grid** (suffix ``_aux``, same resolution as ``_fp``), from which the fields are gathered from the grids to particle positions. For this reason. only ``E`` and ``B`` are defined on this ``_aux`` grid (not the current density or charge density). + +* In some conditions, i.e., when buffers are used for the field gather (for numerical reasons), a **copy of E and B on the auxiliary grid** ``_aux`` **of the level below** ``lev-1`` is stored in fields with suffix ``_cax`` (for `coarse aux`). + +As an example, the structures for the electric field are ``Efield_fp``, ``Efield_cp``, ``Efield_aux`` (and optionally ``Efield_cax``). + +Declaration +----------- + +All the fields described above are public members of class ``WarpX``, defined in ``WarpX.H``. They are defined as an ``amrex::Vector`` (over MR levels) of ``std::array`` (for the 3 spatial components :math:`E_x`, :math:`E_y`, :math:`E_z`) of ``std::unique_ptr`` of ``amrex::MultiFab``, i.e.: + + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_fp; + +Hence, ``Ex`` on MR level ``lev`` is a pointer to an ``amrex::MultiFab``. The other fields are organized in the same way. + +Allocation and initialization +----------------------------- + +The ``MultiFab`` constructor (for, e.g., ``Ex`` on level ``lev``) is called in ``WarpX::AllocLevelMFs``. + +By default, the ``MultiFab`` are set to ``0`` at initialization. They can be assigned a different value in ``WarpX::InitLevelData``. + +Field solver +------------ + +The field solver is performed in ``WarpX::EvolveE`` for the electric field and ``WarpX::EvolveB`` for the magnetic field, called from ``WarpX::OneStep_nosub`` in ``WarpX::EvolveEM``. This section describes the FDTD field push (the PSATD field push should be added later). It is implemented in ``Source/FieldSolver/WarpXPushFieldsEM.cpp``. + +As all cell-wise operation, the field push is done as follows (this is split in multiple functions in the actual implementation to aboid code duplication) +: + +.. code-block:: cpp + + // Loop over MR levels + for (int lev = 0; lev <= finest_level; ++lev) { + // Get pointer to MultiFab Ex on level lev + MultiFab* Ex = Efield_fp[lev][0].get(); + // Loop over boxes (or tiles if not on GPU) + for ( MFIter mfi(*Ex, TilingIfNotGPU()); mfi.isValid(); ++mfi ) { + // Apply field solver on the FAB + } + } + +The innermost step ``// Apply field solver on the FAB`` could be done with 3 nested ``for`` loops for the 3 dimensions (in 3D). However, for portability reasons (see section :ref:`Developers: Portability <developers-portability>`), this is done in two steps: (i) extract AMReX data structures into plain-old-data simple structures, and (ii) call a general ``ParallelFor`` function (translated into nested loops on CPU or a kernel launch on GPU, for instance): + +.. code-block:: cpp + + // Get Box corresponding to the current MFIter + const Box& tex = mfi.tilebox(Ex_nodal_flag); + // Extract the FArrayBox into a simple structure, for portability + Array4<Real> const& Exfab = Ex->array(mfi); + // Loop over cells and perform stencil operation + amrex::ParallelFor(tex, + [=] AMREX_GPU_DEVICE (int j, int k, int l) + { + warpx_push_ex_yee(...); + } + ); + +Function ``warpx_push_ex_yee`` performs the FDTD stencil operation on a single cell. It is implemented in ``Source/FieldSolver/WarpX_K.H`` (where ``_K`` stands for kernel). + +Guard cells exchanges +--------------------- + +Communications are mostly handled in ``Source/Parallelization/``. + +For E and B guard cell **exchanges**, the main functions are variants of ``amrex::FillBoundary(amrex::MultiFab, ...)`` (or ``amrex::MultiFab::FillBoundary(...)``) that fill guard cells of all ``amrex::FArrayBox`` in an ``amrex::MultiFab`` with valid cells of corresponding ``amrex::FArrayBox`` neighbors of the same ``amrex::MultiFab``. There are a number of ``FillBoundaryE``, ``FillBoundaryB`` etc. Under the hood, ``amrex::FillBoundary`` calls ``amrex::ParallelCopy``, which is also sometimes directly called in WarpX. Most calls a + +For the current density, the valid cells of neighboring ``MultiFabs`` are accumulated (added) rather than just copied. This is done using ``amrex::MultiFab::SumBoundary``, and mostly located in ``Source/Parallelization/WarpXSumGuardCells.H``. + +Interpolations for MR +--------------------- + +This is mostly implemented in ``Source/Parallelization``, see the following functions (you may complain to the authors if the documentation is empty) + +.. doxygenfunction:: WarpX::SyncCurrent + +.. doxygenfunction:: interpolateCurrentFineToCoarse + +.. doxygenfunction:: WarpX::RestrictCurrentFromFineToCoarsePatch + +.. doxygenfunction:: WarpX::AddCurrentFromFineLevelandSumBoundary + +Filter +------ + +General functions for filtering can be found in ``Source/Filter/``, where the main ``Filter`` class is defined (see below). All filters (so far there are two of them) in WarpX derive from this class. + +.. doxygenclass:: Filter + +Bilinear filter +~~~~~~~~~~~~~~~ + +The multi-pass bilinear filter (applied on the current density) is implemented in ``Source/Filter/``, and class ``WarpX`` holds an instance of this class in member variable ``WarpX::bilinear_filter``. For performance reasons (to avoid creating too many guard cells), this filter is directly applied in communication routines, see + +.. doxygenfunction:: WarpX::AddCurrentFromFineLevelandSumBoundary + +and + +.. doxygenfunction:: WarpX::ApplyFilterandSumBoundaryJ + +Godfrey's anti-NCI filter for FDTD simulations +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This filter is applied on the electric and magnetic field (on the auxiliary grid) to suppress the Numerical Cherenkov Instability when running FDTD. It is implemented in ``Source/Filter/``, and there are two different stencils, one for ``Ex``, ``Ey`` and ``Bz`` and the other for ``Ez``, ``Bx`` and ``By``. + +.. doxygenclass:: NCIGodfreyFilter + +The class ``WarpX`` holds two corresponding instances of this class in member variables ``WarpX::nci_godfrey_filter_exeybz`` and ``WarpX::nci_godfrey_filter_bxbyez``. It is a 9-point stencil (is the ``z`` direction only), for which the coefficients are computed using tabulated values (depending on dz/dx) in ``Source/Utils/NCIGodfreyTables.H``, see variable ``table_nci_godfrey_galerkin_Ex_Ey_Bz``. The filter is applied in ``PhysicalParticleContainer::Evolve``, right after field gather and before particle push, see + +.. doxygenfunction:: PhysicalParticleContainer::applyNCIFilter diff --git a/Docs/source/developers/initialization.rst b/Docs/source/developers/initialization.rst new file mode 100644 index 000000000..79e94d4d2 --- /dev/null +++ b/Docs/source/developers/initialization.rst @@ -0,0 +1,22 @@ +.. _developers-initialization: + +Initialization +============== + +.. note:: + Section almost empty!! + +General simulation initialization +--------------------------------- + +Regular simulation +~~~~~~~~~~~~~~~~~~ + +Running in a boosted frame +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Field initialization +-------------------- + +Particle initialization +----------------------- diff --git a/Docs/source/developers/moving_window.rst b/Docs/source/developers/moving_window.rst new file mode 100644 index 000000000..38e17c962 --- /dev/null +++ b/Docs/source/developers/moving_window.rst @@ -0,0 +1,8 @@ +.. _developers-moving-window: + +Moving Window +============= + +.. note:: + + Section empty! diff --git a/Docs/source/developers/particles.rst b/Docs/source/developers/particles.rst new file mode 100644 index 000000000..edd817b29 --- /dev/null +++ b/Docs/source/developers/particles.rst @@ -0,0 +1,108 @@ +.. _developers-particles: + +Particles +========= + +Particle containers +------------------- + +Particle structures and functions are defined in ``Source/Particles/``. WarpX uses the ``Particle`` class from AMReX for single particles. An ensemble of particles (e.g., a plasma species, or laser particles) is stored as a ``WarpXParticleContainer`` (see description below) in a per-box (and even per-tile on CPU) basis. + +.. doxygenclass:: WarpXParticleContainer + +Physical species are stored in ``PhysicalParticleContainer``, that derives from ``WarpXParticleContainer``. In particular, the main function to advance all particles in a physical species is ``PhysicalParticleContainer::Evolve`` (see below). + +.. doxygenfunction:: PhysicalParticleContainer::Evolve + :outline: + +Finally, all particle species (physical plasma species ``PhysicalParticleContainer``, photon species ``PhotonParticleContainer`` or non-physical species ``LaserParticleContainer``) are stored in ``MultiParticleContainer``. The class ``WarpX`` holds one instance of ``MultiParticleContainer`` as a member variable, called ``WarpX::mypc`` (where `mypc` stands for "my particle containers"): + +.. doxygenclass:: MultiParticleContainer + +Loop over particles +------------------- + +A typical loop over particles reads: + +.. code-block:: cpp + + // pc is a std::unique_ptr<WarpXParticleContainer> + // Loop over MR levels + for (int lev = 0; lev <= finest_level; ++lev) { + // Loop over particles, box by box + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { + // Do something on particles + // [MY INNER LOOP] + } + } + +The innermost step ``[MY INNER LOOP]`` typically calls ``amrex::ParallelFor`` to perform operations on all particles in a portable way. For this reasons, the particle data needs to be converted in plain-old-data structures. The innermost loop in the code snippet above could look like: + +.. code-block:: cpp + + // Get Array-Of-Struct particle data, also called data + // (x, y, z, id, cpu) + const auto& particles = pti.GetArrayOfStructs(); + // Get Struct-Of-Array particle data, also called attribs + // (ux, uy, uz, w, Exp, Ey, Ez, Bx, By, Bz) + auto& attribs = pti.GetAttribs(); + auto& Exp = attribs[PIdx::Ex]; + // [...] + // Number of particles in this box + const long np = pti.numParticles(); + +Link fields and particles? +-------------------------- + +In WarpX, the loop over boxes through a ``MultiFab`` iterator ``MFIter`` and the loop over boxes through a ``ParticleContainer`` iterator ``WarpXParIter`` are consistent. + +On a loop over boxes in a ``MultiFab`` (``MFIter``), it can be useful to access particle data on a GPU-friendly way. This can be done by: + +.. code-block:: cpp + + // Index of grid (= box) + const int grid_id = mfi.index(); + // Index of tile within the grid + const int tile_id = mfi.LocalTileIndex(); + // Get GPU-friendly arrays of particle data + auto& ptile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + ParticleType* pp = particle_tile.GetArrayOfStructs()().data(); + // Only need attribs (i.e., SoA data) + auto& soa = ptile.GetStructOfArrays(); + // As an example, let's get the ux momentum + const ParticleReal * const AMREX_RESTRICT ux = soa.GetRealData(PIdx::ux).data(); + +On a loop over particles it can be useful to access the fields on the box we are looping over (typically when we use both field and particle data on the same box, for field gather or current deposition for instance). This is done for instance by adding this snippet in ``[MY INNER LOOP]``: + +.. code-block:: cpp + + // E is a reference to, say, WarpX::Efield_aux + // Get the Ex field on the grid + const FArrayBox& exfab = (*E[lev][0])[pti]; + // Let's be generous and also get the underlying box (i.e., index info) + const Box& box = pti.validbox(); + +Main functions +-------------- + +.. doxygenfunction:: PhysicalParticleContainer::FieldGather + +.. doxygenfunction:: PhysicalParticleContainer::PushPX + +.. doxygenfunction:: WarpXParticleContainer::DepositCurrent + +.. note:: + The current deposition is used both by ``PhysicalParticleContainer`` and ``LaserParticleContainer``, so it is in the parent class ``WarpXParticleContainer``. + +Buffers +------- + +To reduce numerical artifacts at the boundary of a mesh-refinement patch, WarpX has an option to use buffers: When particles evolve on the fine level, they gather from the coarse level (e.g., ``Efield_cax``, a copy of the ``aux`` data from the level below) if they are located on the fine level but fewer than ``WarpX::n_field_gather_buffer`` cells away from the coarse-patch boundary. Similarly, when particles evolve on the fine level, they deposit on the coarse level (e.g., ``Efield_cp``) if they are located on the fine level but fewer than ``WarpX::n_current_deposition_buffer`` cells away from the coarse-patch boundary. + +``WarpX::gather_buffer_masks`` and ``WarpX::current_buffer_masks`` contain masks indicating if a cell is in the interior of the fine-resolution patch or in the buffers. Then, particles depending on this mask in + +.. doxygenfunction:: PhysicalParticleContainer::PartitionParticlesInBuffers + +.. note:: + + Buffers are complex! diff --git a/Docs/source/developers/portability.rst b/Docs/source/developers/portability.rst new file mode 100644 index 000000000..c20d6f41f --- /dev/null +++ b/Docs/source/developers/portability.rst @@ -0,0 +1,7 @@ +.. _developers-portability: + +Portability +=========== + +.. note:: + Section empty! diff --git a/Docs/source/developers/profiling.rst b/Docs/source/developers/profiling.rst new file mode 100644 index 000000000..86c838e8d --- /dev/null +++ b/Docs/source/developers/profiling.rst @@ -0,0 +1,7 @@ +.. _developers-profiling: + +Profiling +========= + +.. note:: + Section empty! diff --git a/Docs/source/developers/repo_organization.rst b/Docs/source/developers/repo_organization.rst new file mode 100644 index 000000000..495cc8fc7 --- /dev/null +++ b/Docs/source/developers/repo_organization.rst @@ -0,0 +1,24 @@ +.. _developers-repo-structure: + +WarpX structure +=============== + +Repo organization +----------------- + +All the WarpX source code is located in ``Source/``. All sub-directories have a pretty straigtforward name. The PIC loop is part of the WarpX class, in function ``WarpX::EvolveEM`` implemented in ``Source/WarpXEvolveEM.cpp``. The core of the PIC loop (i.e., without diagnostics etc.) is in ``WarpX::OneStep_nosub`` (when subcycling is OFF) or ``WarpX::OneStep_sub1`` (when subcycling is ON, with method 1). + +Code organization +----------------- + +The main WarpX class is WarpX, implemented in ``Source/WarpX.cpp``. + +Build system +------------ + +WarpX uses the AMReX build system. Each sub-folder contains a file ``Make.package`` with the names of header files and source files to include for the build. + +WarpX-specific vocabulary +------------------------- + +- ``Evolve`` is a generic term to advance a quantity (this comes from AMReX). For instance, ``WarpX::EvolveE(dt)`` advances the electric field for duration ``dt``, ``PhysicalParticleContainer::Evolve(...)`` does field gather + particle push + current deposition for all particles in ``PhysicalParticleContainer``, and ``WarpX::EvolveEM`` is the central ``WarpX`` function that performs 1 PIC iteration. diff --git a/Docs/source/index.rst b/Docs/source/index.rst index 46b695526..b2a8b05a7 100644 --- a/Docs/source/index.rst +++ b/Docs/source/index.rst @@ -37,3 +37,5 @@ In order to learn to use the code, please see the sections below: running_python/running_python visualization/visualization theory/theory + developers/developers + api/library_root diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 2c1d31990..2e6acfdef 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -641,6 +641,34 @@ Laser initialization the field solver. In particular, do not use any other boundary condition than periodic. +Collision initialization +------------------------ + +WarpX provides a relativistic elastic Monte Carlo binary collision model, +following the algorithm given by `Perez et al. (Phys. Plasmas 19, 083104, 2012) <https://doi.org/10.1063/1.4742167>`_. + +* ``collisions.ncollisions`` (`int`) optional (default `0`) + Number of collision types. + +* ``collisions.collision_names`` (`strings`, separated by spaces) + The name of each collision type. It must be provided if ``collisions.ncollisions`` is not zero. + This is then used in the rest of the input deck; + in this documentation we use ``<collision_name>`` as a placeholder. + The number of strings provided should match the number of collision types, + i.e. ``collisions.ncollisions``. + +* ``<collision_name>.species`` (`strings`, two species names separated by spaces) + The names of two species, between which the collision will be considered. + It must be provided if ``collisions.ncollisions`` is not zero, and + the number of provided ``<collision_name>.species`` should match + the number of collision types, i.e. ``collisions.ncollisions``. + +* ``<collision_name>.CoulombLog`` (`float`) optional + A provided fixed Coulomb logarithm of the collision type + ``<collision_name>``. + If this is not provided, or if a non-positive value is provided, + a Coulomb logarithm will be computed automatically according to the algorithm. + Numerics and algorithms ----------------------- diff --git a/Examples/Tests/collision/analysis_collision.py b/Examples/Tests/collision/analysis_collision.py new file mode 100755 index 000000000..25fa26e4b --- /dev/null +++ b/Examples/Tests/collision/analysis_collision.py @@ -0,0 +1,67 @@ +#! /usr/bin/env python + +# This script tests the collision module +# using electron-ion temperature relaxation. +# Initially, electrons and ions are both in equilibrium +# (gaussian) distributions, but have different temperatures. +# Relaxation occurs to bring the two temeratures to be +# a final same temperature through collisions. +# The code was tested to be valid, more detailed results +# were used to obtian an exponential fit with +# coefficients a and b. +# This automated test compares the results with the fit. + +# Possible errors: +# tolerance: 0.001 +# Possible running time: ~ 30.0 s + +import sys +import yt +import re +import math +import statistics +from glob import glob + +tolerance = 0.001 + +ng = 512 +ne = ng * 200 +ni = ng * 200 +np = ne + ni + +c = 299792458.0 +me = 9.10938356e-31 +mi = me * 5.0 + +# exponential fit coefficients +a = 0.041817463099883 +b = -0.083851393560288 + +last_fn = sys.argv[1] +temp = re.compile("([a-zA-Z_]+)([0-9]+)") +res = temp.match(last_fn).groups() +fn_list = glob(res[0] + "?????") + +error = 0.0 +nt = 0 +for fn in fn_list: + # load file + ds = yt.load( fn ) + ad = ds.all_data() + px = ad['particle_momentum_x'].to_ndarray() + # get time index j + buf = temp.match(fn).groups() + j = int(buf[1]) + # compute error + vxe = statistics.mean(px[ 0:ne])/me/c + vxi = statistics.mean(px[ne:np])/mi/c + vxd = vxe - vxi + fit = a*math.exp(b*j) + error = error + abs(fit-vxd) + nt = nt + 1 + +error = error / nt + +print('error = ', error) +print('tolerance = ', tolerance) +assert(error < tolerance) diff --git a/Examples/Tests/collision/inputs b/Examples/Tests/collision/inputs new file mode 100644 index 000000000..906f50b21 --- /dev/null +++ b/Examples/Tests/collision/inputs @@ -0,0 +1,65 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 150 +amr.n_cell = 8 8 8 +amr.max_grid_size = 8 +amr.blocking_factor = 8 +amr.max_level = 0 +geometry.coord_sys = 0 +geometry.is_periodic = 1 1 1 +geometry.prob_lo = 0. 0. 0. +geometry.prob_hi = 4.154046151855669e2 4.154046151855669e2 4.154046151855669e2 +warpx.do_pml = 0 + +################################# +############ NUMERICS ########### +################################# +warpx.serialize_ics = 1 +warpx.verbose = 1 +warpx.cfl = 1.0 +amr.plot_int = 10 +warpx.plot_raw_fields = 0 + +################################# +############ PLASMA ############# +################################# +particles.nspecies = 2 +particles.species_names = electron ion + +electron.charge = -q_e +electron.mass = m_e +electron.injection_style = "NRandomPerCell" +electron.num_particles_per_cell = 200 +electron.profile = constant +electron.density = 1.0e21 +electron.momentum_distribution_type = "gaussian" +electron.ux_th = 0.044237441120300 +electron.uy_th = 0.044237441120300 +electron.uz_th = 0.044237441120300 +electron.ux_m = 0.044237441120300 +electron.do_not_deposit = 1 + +ion.charge = q_e +ion.mass = 4.554691780000000e-30 +ion.injection_style = "NRandomPerCell" +ion.num_particles_per_cell = 200 +ion.profile = constant +ion.density = 1.0e21 +ion.momentum_distribution_type = "gaussian" +ion.ux_th = 0.006256118919701 +ion.uy_th = 0.006256118919701 +ion.uz_th = 0.006256118919701 +ion.do_not_deposit = 1 + +################################# +############ COLLISION ########## +################################# +collisions.ncollisions = 3 +collisions.collision_names = collision1 collision2 collision3 +collision1.species = electron ion +collision2.species = electron electron +collision3.species = ion ion +collision1.CoulombLog = 15.9 +collision2.CoulombLog = 15.9 +collision3.CoulombLog = 15.9 diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 0b78ddaf6..d1190b08c 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -53,7 +53,7 @@ branch = dev [extra-PICSAR] dir = /home/regtester/AMReX_RegTesting/picsar/ -branch = master +branch = QED # individual problems follow @@ -713,8 +713,7 @@ numthreads = 2 compileTest = 0 doVis = 0 compareParticles = 0 -analysisRoutine = Examples/Tests/photon_pusher/analysis_photon_pusher.py - +analysisRoutine = Examples/Tests/photon_pusher/analysis_photon_pusher.py [radiation_reaction] buildDir = . @@ -855,7 +854,7 @@ compileTest = 0 doVis = 0 compareParticles = 1 particleTypes = beam -outputFile = diags/plotfiles/plt00010 +outputFile = diags/plotfiles/plt00002 [PlasmaAccelerationBoost3d] buildDir = . @@ -864,6 +863,7 @@ runtime_params = warpx.do_dynamic_scheduling=0 warpx.serialize_ics=1 amr.n_cell= dim = 3 addToCompileString = restartTest = 0 +tolerance = 1e-14 useMPI = 1 numprocs = 2 useOMP = 1 @@ -1006,3 +1006,17 @@ useOMP = 1 numthreads = 2 compileTest = 0 doVis = 0 + +[collision] +buildDir = . +inputFile = Examples/Tests/collision/inputs +dim = 3 +restartTest = 0 +useMPI = 1 +numprocs = 1 +useOMP = 1 +numthreads = 1 +compileTest = 0 +doVis = 0 +compareParticles = 0 +analysisRoutine = Examples/Tests/collision/analysis_collision.py diff --git a/Source/BoundaryConditions/PML.cpp b/Source/BoundaryConditions/PML.cpp index 572030f73..51439430d 100644 --- a/Source/BoundaryConditions/PML.cpp +++ b/Source/BoundaryConditions/PML.cpp @@ -752,14 +752,14 @@ PML::CopyJtoPMLs (const std::array<amrex::MultiFab*,3>& j_fp, void -PML::ExchangeF (MultiFab* F_fp, MultiFab* F_cp, int do_pml_in_domain) +PML::ExchangeF (amrex::MultiFab* F_fp, amrex::MultiFab* F_cp, int do_pml_in_domain) { ExchangeF(PatchType::fine, F_fp, do_pml_in_domain); ExchangeF(PatchType::coarse, F_cp, do_pml_in_domain); } void -PML::ExchangeF (PatchType patch_type, MultiFab* Fp, int do_pml_in_domain) +PML::ExchangeF (PatchType patch_type, amrex::MultiFab* Fp, int do_pml_in_domain) { if (patch_type == PatchType::fine && pml_F_fp && Fp) { Exchange(*pml_F_fp, *Fp, *m_geom, do_pml_in_domain); @@ -844,11 +844,10 @@ PML::Exchange (MultiFab& pml, MultiFab& reg, const Geometry& geom, void PML::CopyToPML (MultiFab& pml, MultiFab& reg, const Geometry& geom) { - const IntVect& ngr = reg.nGrowVect(); const IntVect& ngp = pml.nGrowVect(); const auto& period = geom.periodicity(); - pml.ParallelCopy(reg, 0, 0, 1, ngr, ngp, period); + pml.ParallelCopy(reg, 0, 0, 1, IntVect(0), ngp, period); } void diff --git a/Source/Diagnostics/WarpXOpenPMD.H b/Source/Diagnostics/WarpXOpenPMD.H index b90df6945..c79a12066 100644 --- a/Source/Diagnostics/WarpXOpenPMD.H +++ b/Source/Diagnostics/WarpXOpenPMD.H @@ -119,6 +119,7 @@ private: /** This function saves the plot file * * @param[in] pc WarpX particle container + * @param[in] name species name * @param[in] iteration timestep * @param[in] write_real_comp The real attribute ids, from WarpX * @param[in] real_comp_names The real attribute names, from WarpX diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index b5fd52bdc..f5491ffe3 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -280,6 +280,7 @@ WarpX::OneStep_nosub (Real cur_time) // Loop over species. For each ionizable species, create particles in // product species. mypc->doFieldIonization(); + mypc->doCoulombCollisions(); // Push particle from x^{n} to x^{n+1} // from p^{n-1/2} to p^{n+1/2} // Deposit current j^{n+1/2} @@ -475,7 +476,7 @@ WarpX::OneStep_sub1 (Real curtime) } void -WarpX::PushParticlesandDepose (Real cur_time) +WarpX::PushParticlesandDepose (amrex::Real cur_time) { // Evolve particles to p^{n+1/2} and x^{n+1} // Depose current, j^{n+1/2} @@ -485,7 +486,7 @@ WarpX::PushParticlesandDepose (Real cur_time) } void -WarpX::PushParticlesandDepose (int lev, Real cur_time, DtType a_dt_type) +WarpX::PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type) { mypc->Evolve(lev, *Efield_aux[lev][0],*Efield_aux[lev][1],*Efield_aux[lev][2], @@ -583,7 +584,7 @@ WarpX::ComputeDt () * simulation box passes input parameter zmax_plasma_to_compute_max_step. */ void -WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ +WarpX::computeMaxStepBoostAccelerator(const amrex::Geometry& a_geom){ // Sanity checks: can use zmax_plasma_to_compute_max_step only if // the moving window and the boost are all in z direction. AMREX_ALWAYS_ASSERT_WITH_MESSAGE( diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp index 8f0853484..edd4df34d 100644 --- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp +++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp @@ -3,10 +3,10 @@ using namespace amrex; /* \brief Initialize fields in spectral space, and FFT plans */ -SpectralFieldData::SpectralFieldData( const BoxArray& realspace_ba, - const SpectralKSpace& k_space, - const DistributionMapping& dm, - const int n_field_required ) +SpectralFieldData::SpectralFieldData( const amrex::BoxArray& realspace_ba, + const SpectralKSpace& k_space, + const amrex::DistributionMapping& dm, + const int n_field_required ) { const BoxArray& spectralspace_ba = k_space.spectralspace_ba; diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index 9807665c6..4848b051e 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -90,7 +90,7 @@ WarpX::PushPSATD_localFFT (int lev, amrex::Real /* dt */) #endif void -WarpX::EvolveB (Real a_dt) +WarpX::EvolveB (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { EvolveB(lev, a_dt); @@ -98,7 +98,7 @@ WarpX::EvolveB (Real a_dt) } void -WarpX::EvolveB (int lev, Real a_dt) +WarpX::EvolveB (int lev, amrex::Real a_dt) { BL_PROFILE("WarpX::EvolveB()"); EvolveB(lev, PatchType::fine, a_dt); @@ -303,7 +303,7 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) } void -WarpX::EvolveE (Real a_dt) +WarpX::EvolveE (amrex::Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { @@ -312,7 +312,7 @@ WarpX::EvolveE (Real a_dt) } void -WarpX::EvolveE (int lev, Real a_dt) +WarpX::EvolveE (int lev, amrex::Real a_dt) { BL_PROFILE("WarpX::EvolveE()"); EvolveE(lev, PatchType::fine, a_dt); @@ -611,7 +611,7 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) } void -WarpX::EvolveF (Real a_dt, DtType a_dt_type) +WarpX::EvolveF (amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -622,7 +622,7 @@ WarpX::EvolveF (Real a_dt, DtType a_dt_type) } void -WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) +WarpX::EvolveF (int lev, amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -631,7 +631,7 @@ WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) } void -WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) +WarpX::EvolveF (int lev, PatchType patch_type, amrex::Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 5d3c14c14..93729a0ae 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -144,7 +144,7 @@ void Filter::DoFilter (const Box& tbx, * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil()"); ncomp = std::min(ncomp, srcmf.nComp()); @@ -179,8 +179,8 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco * \param ncomp Number of components on which the filter is applied. */ void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, - const Box& tbx, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (amrex::FArrayBox& dstfab, const amrex::FArrayBox& srcfab, + const amrex::Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.nComp()); diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index e32719adf..48aaac275 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -1,3 +1,6 @@ +#ifndef WARPX_F_H_ +#define WARPX_F_H_ + #include <AMReX_BLFort.H> #ifdef __cplusplus @@ -130,3 +133,5 @@ extern "C" #ifdef __cplusplus } #endif + +#endif //WARPX_F_H_ diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 300c41f7b..261a14fc5 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -21,7 +21,7 @@ namespace { amrex::Abort(string.c_str()); } - Real parseChargeName(const ParmParse pp, const std::string& name) { + Real parseChargeName(const ParmParse& pp, const std::string& name) { Real result; if (name == "q_e") { return PhysConst::q_e; @@ -33,13 +33,13 @@ namespace { } } - Real parseChargeString(const ParmParse pp, const std::string& name) { + Real parseChargeString(const ParmParse& pp, const std::string& name) { if(name.substr(0, 1) == "-") return -1.0 * parseChargeName(pp, name.substr(1, name.size() - 1)); return parseChargeName(pp, name); } - Real parseMassString(const ParmParse pp, const std::string& name) { + Real parseMassString(const ParmParse& pp, const std::string& name) { Real result; if (name == "m_e") { return PhysConst::m_e; @@ -109,6 +109,7 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name) part_pos_s.end(), part_pos_s.begin(), ::tolower); + num_particles_per_cell_each_dim.assign(3, 0); if (part_pos_s == "python") { return; } else if (part_pos_s == "singleparticle") { diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index 63ace31fb..1e83ca02f 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -10,6 +10,17 @@ #include <memory> #include <limits> +/** + * The main method to inject a laser pulse in WarpX is to use an artificial + * antenna: particles evenly distributed in a given plane (one particle per + * cell) move at each iteration and deposit a current J onto the grid, which + * in turns creates an electromagnetic field on the grid. The particles' + * displacements are prescribed to create the field requested by the user. + * + * These artificial particles are contained in the LaserParticleContainer. + * LaserParticleContainer derives directly from WarpXParticleContainer. It + * requires a DepositCurrent function, but no FieldGather function. + */ class LaserParticleContainer : public WarpXParticleContainer { diff --git a/Source/Laser/LaserProfiles.H b/Source/Laser/LaserProfiles.H index 528309492..e0ec0dc28 100644 --- a/Source/Laser/LaserProfiles.H +++ b/Source/Laser/LaserProfiles.H @@ -45,7 +45,7 @@ public: /** Initialize Laser Profile * * Reads the section of the inputfile relative to the laser beam - * (e.g. <laser_name>.profile_t_peak, <laser_name>.profile_duration...) + * (e.g. laser_name.profile_t_peak, laser_name.profile_duration...) * and the "my_constants" section. It also receives some common * laser profile parameters. It uses these data to initialize the * member variables of the laser profile class. @@ -64,6 +64,7 @@ public: * * Xp, Yp and amplitude must be arrays with the same length * + * @param[in] np number of antenna particles * @param[in] Xp X coordinate of the particles of the antenna * @param[in] Yp Y coordinate of the particles of the antenna * @param[in] t time (seconds) diff --git a/Source/Particles/Collision/CollisionType.H b/Source/Particles/Collision/CollisionType.H new file mode 100644 index 000000000..d020f47e8 --- /dev/null +++ b/Source/Particles/Collision/CollisionType.H @@ -0,0 +1,39 @@ +#ifndef WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ +#define WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ + +#include "WarpXParticleContainer.H" +#include <AMReX_DenseBins.H> +#include <AMReX_REAL.H> +#include <AMReX_ParmParse.H> + +class CollisionType +{ +public: + int m_species1_index; + int m_species2_index; + bool m_isSameSpecies; + amrex::Real m_CoulombLog; + + CollisionType( + const std::vector<std::string>& species_names, + std::string const collision_name); + + /** Perform all binary collisions within a tile + * + * @param lev AMR level of the tile + * @param mfi iterator for multifab + * @param species1/2 pointer to species container + * @param isSameSpecies true if collision is between same species + * @param CoulombLog user input Coulomb logrithm + * + */ + + static void doCoulombCollisionsWithinTile ( + int const lev, amrex::MFIter const& mfi, + std::unique_ptr<WarpXParticleContainer>& species1, + std::unique_ptr<WarpXParticleContainer>& species2, + bool const isSameSpecies, amrex::Real const CoulombLog ); + +}; + +#endif // WARPX_PARTICLES_COLLISION_COLLISIONTYPE_H_ diff --git a/Source/Particles/Collision/CollisionType.cpp b/Source/Particles/Collision/CollisionType.cpp new file mode 100644 index 000000000..b8014579d --- /dev/null +++ b/Source/Particles/Collision/CollisionType.cpp @@ -0,0 +1,241 @@ +#include "CollisionType.H" +#include "ShuffleFisherYates.H" +#include "ElasticCollisionPerez.H" +#include <WarpX.H> + +CollisionType::CollisionType( + const std::vector<std::string>& species_names, + std::string const collision_name) +{ + +#if defined WARPX_DIM_XZ + amrex::Abort("Collisions only work in 3D geometry for now."); +#elif defined WARPX_DIM_RZ + amrex::Abort("Collisions only work in Cartesian geometry for now."); +#endif + + // read collision species + std::vector<std::string> collision_species; + amrex::ParmParse pp(collision_name); + pp.getarr("species", collision_species); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(collision_species.size() == 2, + "Collision species must name exactly two species."); + + // default Coulomb log, if < 0, will be computed automatically + m_CoulombLog = -1.0; + pp.query("CoulombLog", m_CoulombLog); + + for (int i=0; i<species_names.size(); i++) + { + if (species_names[i] == collision_species[0]) + { m_species1_index = i; } + if (species_names[i] == collision_species[1]) + { m_species2_index = i; } + } + + if (collision_species[0] == collision_species[1]) + m_isSameSpecies = true; + else + m_isSameSpecies = false; + +} + +using namespace amrex; +// Define shortcuts for frequently-used type names +using ParticleType = WarpXParticleContainer::ParticleType; +using ParticleTileType = WarpXParticleContainer::ParticleTileType; +using ParticleBins = DenseBins<ParticleType>; +using index_type = ParticleBins::index_type; + +namespace { + + /* Find the particles and count the particles that are in each cell. + Note that this does *not* rearrange particle arrays */ + ParticleBins + findParticlesInEachCell( int const lev, MFIter const& mfi, + ParticleTileType const& ptile) { + + // Extract particle structures for this tile + int const np = ptile.numParticles(); + ParticleType const* particle_ptr = ptile.GetArrayOfStructs()().data(); + + // Extract box properties + Geometry const& geom = WarpX::GetInstance().Geom(lev); + Box const& cbx = mfi.tilebox(IntVect::TheZeroVector()); //Cell-centered box + const auto lo = lbound(cbx); + const auto dxi = geom.InvCellSizeArray(); + const auto plo = geom.ProbLoArray(); + + // Find particles that are in each cell ; + // results are stored in the object `bins`. + ParticleBins bins; + bins.build(np, particle_ptr, cbx, + // Pass lambda function that returns the cell index + [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p) noexcept -> IntVect + { + return IntVect(AMREX_D_DECL((p.pos(0)-plo[0])*dxi[0] - lo.x, + (p.pos(1)-plo[1])*dxi[1] - lo.y, + (p.pos(2)-plo[2])*dxi[2] - lo.z)); + }); + + return bins; + } + +} + +/** Perform all binary collisions within a tile + * + * @param lev AMR level of the tile + * @param mfi iterator for multifab + * @param species1/2 pointer to species container + * @param isSameSpecies true if collision is between same species + * @param CoulombLog user input Coulomb logrithm + * + */ +void CollisionType::doCoulombCollisionsWithinTile + ( int const lev, MFIter const& mfi, + std::unique_ptr<WarpXParticleContainer>& species_1, + std::unique_ptr<WarpXParticleContainer>& species_2, + bool const isSameSpecies, Real const CoulombLog ) +{ + + if ( isSameSpecies ) // species_1 == species_2 + { + // Extract particles in the tile that `mfi` points to + ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi); + + // Find the particles that are in each cell of this tile + ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 ); + + // Loop over cells, and collide the particles in each cell + + // Extract low-level data + int const n_cells = bins_1.numBins(); + // - Species 1 + auto& soa_1 = ptile_1.GetStructOfArrays(); + ParticleReal * const AMREX_RESTRICT ux_1 = + soa_1.GetRealData(PIdx::ux).data(); + ParticleReal * const AMREX_RESTRICT uy_1 = + soa_1.GetRealData(PIdx::uy).data(); + ParticleReal * const AMREX_RESTRICT uz_1 = + soa_1.GetRealData(PIdx::uz).data(); + ParticleReal const * const AMREX_RESTRICT w_1 = + soa_1.GetRealData(PIdx::w).data(); + index_type* indices_1 = bins_1.permutationPtr(); + index_type const* cell_offsets_1 = bins_1.offsetsPtr(); + Real q1 = species_1->getCharge(); + Real m1 = species_1->getMass(); + + const Real dt = WarpX::GetInstance().getdt(lev); + Geometry const& geom = WarpX::GetInstance().Geom(lev); + const Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2); + + // Loop over cells + amrex::ParallelFor( n_cells, + [=] AMREX_GPU_DEVICE (int i_cell) noexcept + { + // The particles from species1 that are in the cell `i_cell` are + // given by the `indices_1[cell_start_1:cell_stop_1]` + index_type const cell_start_1 = cell_offsets_1[i_cell]; + index_type const cell_stop_1 = cell_offsets_1[i_cell+1]; + index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2; + + // Do not collide if there is only one particle in the cell + if ( cell_stop_1 - cell_start_1 >= 2 ) + { + // shuffle + ShuffleFisherYates( + indices_1, cell_start_1, cell_half_1 ); + + // Call the function in order to perform collisions + ElasticCollisionPerez( + cell_start_1, cell_half_1, + cell_half_1, cell_stop_1, + indices_1, indices_1, + ux_1, uy_1, uz_1, ux_1, uy_1, uz_1, w_1, w_1, + q1, q1, m1, m1, Real(-1.0), Real(-1.0), + dt, CoulombLog, dV ); + } + } + ); + } + else // species_1 != species_2 + { + // Extract particles in the tile that `mfi` points to + ParticleTileType& ptile_1 = species_1->ParticlesAt(lev, mfi); + ParticleTileType& ptile_2 = species_2->ParticlesAt(lev, mfi); + + // Find the particles that are in each cell of this tile + ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 ); + ParticleBins bins_2 = findParticlesInEachCell( lev, mfi, ptile_2 ); + + // Loop over cells, and collide the particles in each cell + + // Extract low-level data + int const n_cells = bins_1.numBins(); + // - Species 1 + auto& soa_1 = ptile_1.GetStructOfArrays(); + ParticleReal * const AMREX_RESTRICT ux_1 = + soa_1.GetRealData(PIdx::ux).data(); + ParticleReal * const AMREX_RESTRICT uy_1 = + soa_1.GetRealData(PIdx::uy).data(); + ParticleReal * const AMREX_RESTRICT uz_1 = + soa_1.GetRealData(PIdx::uz).data(); + ParticleReal const * const AMREX_RESTRICT w_1 = + soa_1.GetRealData(PIdx::w).data(); + index_type* indices_1 = bins_1.permutationPtr(); + index_type const* cell_offsets_1 = bins_1.offsetsPtr(); + Real q1 = species_1->getCharge(); + Real m1 = species_1->getMass(); + // - Species 2 + auto& soa_2 = ptile_2.GetStructOfArrays(); + Real* ux_2 = soa_2.GetRealData(PIdx::ux).data(); + Real* uy_2 = soa_2.GetRealData(PIdx::uy).data(); + Real* uz_2 = soa_2.GetRealData(PIdx::uz).data(); + Real* w_2 = soa_2.GetRealData(PIdx::w).data(); + index_type* indices_2 = bins_2.permutationPtr(); + index_type const* cell_offsets_2 = bins_2.offsetsPtr(); + Real q2 = species_2->getCharge(); + Real m2 = species_2->getMass(); + + const Real dt = WarpX::GetInstance().getdt(lev); + Geometry const& geom = WarpX::GetInstance().Geom(lev); + const Real dV = geom.CellSize(0)*geom.CellSize(1)*geom.CellSize(2); + + // Loop over cells + amrex::ParallelFor( n_cells, + [=] AMREX_GPU_DEVICE (int i_cell) noexcept + { + // The particles from species1 that are in the cell `i_cell` are + // given by the `indices_1[cell_start_1:cell_stop_1]` + index_type const cell_start_1 = cell_offsets_1[i_cell]; + index_type const cell_stop_1 = cell_offsets_1[i_cell+1]; + // Same for species 2 + index_type const cell_start_2 = cell_offsets_2[i_cell]; + index_type const cell_stop_2 = cell_offsets_2[i_cell+1]; + + // ux from species1 can be accessed like this: + // ux_1[ indices_1[i] ], where i is between + // cell_start_1 (inclusive) and cell_start_2 (exclusive) + + // Do not collide if one species is missing in the cell + if ( cell_stop_1 - cell_start_1 >= 1 && + cell_stop_2 - cell_start_2 >= 1 ) + { + // shuffle + ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1); + ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2); + + // Call the function in order to perform collisions + ElasticCollisionPerez( + cell_start_1, cell_stop_1, cell_start_2, cell_stop_2, + indices_1, indices_2, + ux_1, uy_1, uz_1, ux_2, uy_2, uz_2, w_1, w_2, + q1, q2, m1, m2, Real(-1.0), Real(-1.0), + dt, CoulombLog, dV ); + } + } + ); + } // end if ( isSameSpecies) + +} diff --git a/Source/Particles/Collision/ComputeTemperature.H b/Source/Particles/Collision/ComputeTemperature.H new file mode 100644 index 000000000..3cc96fb52 --- /dev/null +++ b/Source/Particles/Collision/ComputeTemperature.H @@ -0,0 +1,39 @@ +#ifndef WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_ +#define WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_ + +#include <WarpXConst.H> + +template <typename T_index, typename T_R> +T_R ComputeTemperature ( + T_index const Is, T_index const Ie, T_index const *I, + T_R const *ux, T_R const *uy, T_R const *uz, T_R const m ) +{ + + T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c ); + + int N = Ie - Is; + if ( N == 0 ) { return T_R(0.0); } + + T_R vx = T_R(0.0); T_R vy = T_R(0.0); + T_R vz = T_R(0.0); T_R vs = T_R(0.0); + T_R gm = T_R(0.0); T_R us = T_R(0.0); + + for (int i = Is; i < Ie; ++i) + { + us = ( ux[ I[i] ] * ux[ I[i] ] + + uy[ I[i] ] * uy[ I[i] ] + + uz[ I[i] ] * uz[ I[i] ] ); + gm = std::sqrt( T_R(1.0) + us*inv_c2 ); + vx += ux[ I[i] ] / gm; + vy += uy[ I[i] ] / gm; + vz += uz[ I[i] ] / gm; + vs += us / gm / gm; + } + + vx = vx / N; vy = vy / N; + vz = vz / N; vs = vs / N; + + return m/T_R(3.0)*(vs-(vx*vx+vy*vy+vz*vz)); +} + +#endif // WARPX_PARTICLES_COLLISION_COMPUTE_TEMPERATURE_H_ diff --git a/Source/Particles/Collision/ElasticCollisionPerez.H b/Source/Particles/Collision/ElasticCollisionPerez.H new file mode 100644 index 000000000..8e16d95cc --- /dev/null +++ b/Source/Particles/Collision/ElasticCollisionPerez.H @@ -0,0 +1,105 @@ +#ifndef WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ +#define WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ + +#include "UpdateMomentumPerezElastic.H" +#include "ComputeTemperature.H" +#include <WarpXConst.H> +#include <AMReX_Random.H> + +/** \brief Prepare information for and call + * UpdateMomentumPerezElastic(). + * @param[in] I1s,I2s is the start index for I1,I2 (inclusive). + * @param[in] I1e,I2e is the start index for I1,I2 (exclusive). + * @param[in] I1 and I2 are the index arrays. + * @param[in,out] u1 and u2 are the velocity arrays (u=v*gamma), + * they could be either different or the same, + * their lengths are not needed, + * @param[in] I1 and I2 determine all elements that will be used. + * @param[in] w1 and w2 are arrays of weights. + * @param[in] q1 and q2 are charges. m1 and m2 are masses. + * @param[in] T1 and T2 are temperatures (Joule) + * and will be used if greater than zero, + * otherwise will be computed. + * @param[in] dt is the time step length between two collision calls. + * @param[in] L is the Coulomb log and will be used if greater than zero, + * otherwise will be computed. + * @param[in] dV is the volume of the corresponding cell. +*/ + +template <typename T_index, typename T_R> +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ElasticCollisionPerez ( + T_index const I1s, T_index const I1e, + T_index const I2s, T_index const I2e, + T_index *I1, T_index *I2, + T_R *u1x, T_R *u1y, T_R *u1z, + T_R *u2x, T_R *u2y, T_R *u2z, + T_R const *w1, T_R const *w2, + T_R const q1, T_R const q2, + T_R const m1, T_R const m2, + T_R const T1, T_R const T2, + T_R const dt, T_R const L, T_R const dV) +{ + + T_R constexpr inv_c2 = T_R(1.0)/(PhysConst::c*PhysConst::c); + int NI1 = I1e - I1s; + int NI2 = I2e - I2s; + + // get local T1t and T2t + T_R T1t; T_R T2t; + if ( T1 <= T_R(0.0) && L <= T_R(0.0) ) + { + T1t = ComputeTemperature(I1s,I1e,I1,u1x,u1y,u1z,m1); + } + else { T1t = T1; } + if ( T2 <= T_R(0.0) && L <= T_R(0.0) ) + { + T2t = ComputeTemperature(I2s,I2e,I2,u2x,u2y,u2z,m2); + } + else { T2t = T2; } + + // local density + T_R n1 = T_R(0.0); + T_R n2 = T_R(0.0); + T_R n12 = T_R(0.0); + for (int i1=I1s; i1<I1e; ++i1) { n1 += w1[ I1[i1] ]; } + for (int i2=I2s; i2<I2e; ++i2) { n2 += w2[ I2[i2] ]; } + n1 = n1 / dV; n2 = n2 / dV; + { + int i1 = I1s; int i2 = I2s; + for (int k = 0; k < amrex::max(NI1,NI2); ++k) + { + n12 += amrex::min( w1[ I1[i1] ], w2[ I2[i2] ] ); + ++i1; if ( i1 == I1e ) { i1 = I1s; } + ++i2; if ( i2 == I2e ) { i2 = I2s; } + } + n12 = n12 / dV; + } + + // compute Debye length lmdD + T_R lmdD; + lmdD = T_R(1.0)/std::sqrt( n1*q1*q1/(T1t*PhysConst::ep0) + + n2*q2*q2/(T2t*PhysConst::ep0) ); + T_R rmin = std::pow( T_R(4.0) * MathConst::pi / T_R(3.0) * + amrex::max(n1,n2), T_R(-1.0/3.0) ); + lmdD = amrex::max(lmdD, rmin); + + // call UpdateMomentumPerezElastic() + { + int i1 = I1s; int i2 = I2s; + for (int k = 0; k < amrex::max(NI1,NI2); ++k) + { + UpdateMomentumPerezElastic( + u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ], + u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ], + n1, n2, n12, + q1, m1, w1[ I1[i1] ], q2, m2, w2[ I2[i2] ], + dt, L, lmdD); + ++i1; if ( i1 == I1e ) { i1 = I1s; } + ++i2; if ( i2 == I2e ) { i2 = I2s; } + } + } + +} + +#endif // WARPX_PARTICLES_COLLISION_ELASTIC_COLLISION_PEREZ_H_ diff --git a/Source/Particles/Collision/Make.package b/Source/Particles/Collision/Make.package new file mode 100644 index 000000000..163508fb7 --- /dev/null +++ b/Source/Particles/Collision/Make.package @@ -0,0 +1,10 @@ +CEXE_headers += CollisionType.H +CEXE_headers += ElasticCollisionPerez.H +CEXE_headers += ShuffleFisherYates.H +CEXE_headers += UpdateMomentumPerezElastic.H +CEXE_headers += ComputeTemperature.H + +CEXE_sources += CollisionType.cpp + +INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision +VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles/Collision diff --git a/Source/Particles/Collision/ShuffleFisherYates.H b/Source/Particles/Collision/ShuffleFisherYates.H new file mode 100644 index 000000000..621e654d6 --- /dev/null +++ b/Source/Particles/Collision/ShuffleFisherYates.H @@ -0,0 +1,29 @@ +#ifndef WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ +#define WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ + +#include <AMReX_Random.H> + +/* \brief Shuffle array according to Fisher-Yates algorithm. + * Only shuffle the part between is <= i < ie, n = ie-is. + * T_index shall be + * amrex::DenseBins<WarpXParticleContainer::ParticleType>::index_type +*/ + +template <typename T_index> +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void ShuffleFisherYates (T_index *array, T_index const is, T_index const ie) +{ + int j; + T_index buf; + for (int i = ie-1; i >= is+1; --i) + { + // get random number j: is <= j <= i + j = amrex::Random_int(i-is+1) + is; + // swop the ith array element with the jth + buf = array[i]; + array[i] = array[j]; + array[j] = buf; + } +} + +#endif // WARPX_PARTICLES_COLLISION_SHUFFLE_FISHER_YATES_H_ diff --git a/Source/Particles/Collision/UpdateMomentumPerezElastic.H b/Source/Particles/Collision/UpdateMomentumPerezElastic.H new file mode 100644 index 000000000..948e8b075 --- /dev/null +++ b/Source/Particles/Collision/UpdateMomentumPerezElastic.H @@ -0,0 +1,252 @@ +#ifndef WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ +#define WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ + +#include <WarpXConst.H> +#include <AMReX_Random.H> +#include <cmath> // isnan() isinf() +#include <limits> // numeric_limits<float>::min() + +/* \brief Update particle velocities according to + * F. Perez et al., Phys.Plasmas.19.083104 (2012), + * which is based on Nanbu's method, PhysRevE.55.4642 (1997). + * @param[in] LmdD is max(Debye length, minimal interparticle distance). + * @param[in] L is the Coulomb log. A fixed L will be used if L > 0, + * otherwise L will be calculated based on the algorithm. + * To see if there are nan or inf updated velocities, + * compile with USE_ASSERTION=TRUE. +*/ + +template <typename T_R> +AMREX_GPU_HOST_DEVICE AMREX_INLINE +void UpdateMomentumPerezElastic ( + T_R& u1x, T_R& u1y, T_R& u1z, T_R& u2x, T_R& u2y, T_R& u2z, + T_R const n1, T_R const n2, T_R const n12, + T_R const q1, T_R const m1, T_R const w1, + T_R const q2, T_R const m2, T_R const w2, + T_R const dt, T_R const L, T_R const lmdD) +{ + + // If g = u1 - u2 = 0, do not collide. + if ( std::abs(u1x-u2x) < std::numeric_limits<T_R>::min() && + std::abs(u1y-u2y) < std::numeric_limits<T_R>::min() && + std::abs(u1z-u2z) < std::numeric_limits<T_R>::min() ) + { return; } + + T_R constexpr inv_c2 = T_R(1.0) / ( PhysConst::c * PhysConst::c ); + + // Compute Lorentz factor gamma + T_R const g1 = std::sqrt( T_R(1.0) + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_c2 ); + T_R const g2 = std::sqrt( T_R(1.0) + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_c2 ); + + // Compute momenta + T_R const p1x = u1x * m1; + T_R const p1y = u1y * m1; + T_R const p1z = u1z * m1; + T_R const p2x = u2x * m2; + T_R const p2y = u2y * m2; + T_R const p2z = u2z * m2; + + // Compute center-of-mass (COM) velocity and gamma + T_R const mass_g = m1 * g1 + m2 * g2; + T_R const vcx = (p1x+p2x) / mass_g; + T_R const vcy = (p1y+p2y) / mass_g; + T_R const vcz = (p1z+p2z) / mass_g; + T_R const vcms = vcx*vcx + vcy*vcy + vcz*vcz; + T_R const gc = T_R(1.0) / std::sqrt( T_R(1.0) - vcms*inv_c2 ); + + // Compute vc dot v1 and v2 + T_R const vcDv1 = (vcx*u1x + vcy*u1y + vcz*u1z) / g1; + T_R const vcDv2 = (vcx*u2x + vcy*u2y + vcz*u2z) / g2; + + // Compute p1 star + T_R p1sx; + T_R p1sy; + T_R p1sz; + if ( vcms > std::numeric_limits<T_R>::min() ) + { + T_R const lorentz_tansform_factor = + ( (gc-T_R(1.0))/vcms*vcDv1 - gc )*m1*g1; + p1sx = p1x + vcx*lorentz_tansform_factor; + p1sy = p1y + vcy*lorentz_tansform_factor; + p1sz = p1z + vcz*lorentz_tansform_factor; + } + else // If vcms = 0, don't do Lorentz-transform. + { + p1sx = p1x; + p1sy = p1y; + p1sz = p1z; + } + T_R const p1sm = std::sqrt( p1sx*p1sx + p1sy*p1sy + p1sz*p1sz ); + + // Compute gamma star + T_R const g1s = ( T_R(1.0) - vcDv1*inv_c2 )*gc*g1; + T_R const g2s = ( T_R(1.0) - vcDv2*inv_c2 )*gc*g2; + + // Compute the Coulomb log lnLmd + T_R lnLmd; + if ( L > T_R(0.0) ) { lnLmd = L; } + else + { + // Compute b0 + T_R const b0 = std::abs(q1*q2) * inv_c2 / + (T_R(4.0)*MathConst::pi*PhysConst::ep0) * gc/mass_g * + ( m1*g1s*m2*g2s/(p1sm*p1sm*inv_c2) + T_R(1.0) ); + + // Compute the minimal impact parameter + T_R bmin = amrex::max(PhysConst::hbar*MathConst::pi/p1sm,b0); + + // Compute the Coulomb log lnLmd + lnLmd = amrex::max( T_R(2.0), + T_R(0.5)*std::log(T_R(1.0)+lmdD*lmdD/(bmin*bmin)) ); + } + + // Compute s + T_R s = n1*n2/n12 * dt*lnLmd*q1*q1*q2*q2 / + ( T_R(4.0) * MathConst::pi * PhysConst::ep0 * PhysConst::ep0 * + m1*g1*m2*g2/(inv_c2*inv_c2) ) * gc*p1sm/mass_g * + std::pow(m1*g1s*m2*g2s/(inv_c2*p1sm*p1sm) + T_R(1.0), 2.0); + + // Compute s' + T_R const vrel = mass_g*p1sm/(m1*g1s*m2*g2s*gc); + T_R const sp = std::pow(T_R(4.0)*MathConst::pi/T_R(3.0),T_R(1.0/3.0)) * + n1*n2/n12 * dt * vrel * (m1+m2) / + amrex::max( m1*std::pow(n1,T_R(2.0/3.0)), + m2*std::pow(n2,T_R(2.0/3.0)) ); + + // Determine s + s = amrex::min(s,sp); + + // Get random numbers + T_R r = amrex::Random(); + + // Compute scattering angle + T_R cosXs; + T_R sinXs; + if ( s <= T_R(0.1) ) + { + while ( true ) + { + cosXs = T_R(1.0) + s * std::log(r); + // Avoid the bug when r is too small such that cosXs < -1 + if ( cosXs >= T_R(-1.0) ) { break; } + r = amrex::Random(); + } + } + else if ( s > T_R(0.1) && s <= T_R(3.0) ) + { + T_R const Ainv = 0.0056958 + 0.9560202*s - 0.508139*s*s + + 0.47913906*s*s*s - 0.12788975*s*s*s*s + 0.02389567*s*s*s*s*s; + cosXs = Ainv * std::log( std::exp(T_R(-1.0)/Ainv) + + T_R(2.0) * r * std::sinh(T_R(1.0)/Ainv) ); + } + else if ( s > T_R(3.0) && s <= T_R(6.0) ) + { + T_R const A = T_R(3.0) * std::exp(-s); + cosXs = T_R(1.0)/A * std::log( std::exp(-A) + + T_R(2.0) * r * std::sinh(A) ); + } + else + { + cosXs = T_R(2.0) * r - T_R(1.0); + } + sinXs = std::sqrt(T_R(1.0) - cosXs*cosXs); + + // Get random azimuthal angle + T_R const phis = amrex::Random() * T_R(2.0) * MathConst::pi; + T_R const cosphis = std::cos(phis); + T_R const sinphis = std::sin(phis); + + // Compute post-collision momenta pfs in COM + T_R p1fsx; + T_R p1fsy; + T_R p1fsz; + // p1sp is the p1s perpendicular + T_R p1sp = std::sqrt( p1sx*p1sx + p1sy*p1sy ); + // Make sure p1sp is not almost zero + if ( p1sp > std::numeric_limits<T_R>::min() ) + { + p1fsx = ( p1sx*p1sz/p1sp ) * sinXs*cosphis + + ( p1sy*p1sm/p1sp ) * sinXs*sinphis + + ( p1sx ) * cosXs; + p1fsy = ( p1sy*p1sz/p1sp ) * sinXs*cosphis + + (-p1sx*p1sm/p1sp ) * sinXs*sinphis + + ( p1sy ) * cosXs; + p1fsz = (-p1sp ) * sinXs*cosphis + + ( T_R(0.0) ) * sinXs*sinphis + + ( p1sz ) * cosXs; + // Note a negative sign is different from + // Eq. (12) in Perez's paper, + // but they are the same due to the random nature of phis. + } + else + { + // If the previous p1sp is almost zero + // x->y y->z z->x + // This set is equivalent to the one in Nanbu's paper + p1sp = std::sqrt( p1sy*p1sy + p1sz*p1sz ); + p1fsy = ( p1sy*p1sx/p1sp ) * sinXs*cosphis + + ( p1sz*p1sm/p1sp ) * sinXs*sinphis + + ( p1sy ) * cosXs; + p1fsz = ( p1sz*p1sx/p1sp ) * sinXs*cosphis + + (-p1sy*p1sm/p1sp ) * sinXs*sinphis + + ( p1sz ) * cosXs; + p1fsx = (-p1sp ) * sinXs*cosphis + + ( T_R(0.0) ) * sinXs*sinphis + + ( p1sx ) * cosXs; + } + + T_R const p2fsx = -p1fsx; + T_R const p2fsy = -p1fsy; + T_R const p2fsz = -p1fsz; + + // Transform from COM to lab frame + T_R p1fx; T_R p2fx; + T_R p1fy; T_R p2fy; + T_R p1fz; T_R p2fz; + if ( vcms > std::numeric_limits<T_R>::min() ) + { + T_R const vcDp1fs = vcx*p1fsx + vcy*p1fsy + vcz*p1fsz; + T_R const vcDp2fs = vcx*p2fsx + vcy*p2fsy + vcz*p2fsz; + T_R const factor = (gc-T_R(1.0))/vcms; + T_R const factor1 = factor*vcDp1fs + m1*g1s*gc; + T_R const factor2 = factor*vcDp2fs + m2*g2s*gc; + p1fx = p1fsx + vcx * factor1; + p1fy = p1fsy + vcy * factor1; + p1fz = p1fsz + vcz * factor1; + p2fx = p2fsx + vcx * factor2; + p2fy = p2fsy + vcy * factor2; + p2fz = p2fsz + vcz * factor2; + } + else // If vcms = 0, don't do Lorentz-transform. + { + p1fx = p1fsx; + p1fy = p1fsy; + p1fz = p1fsz; + p2fx = p2fsx; + p2fy = p2fsy; + p2fz = p2fsz; + } + + // Rejection method + r = amrex::Random(); + if ( w2 > r*amrex::max(w1, w2) ) + { + u1x = p1fx / m1; + u1y = p1fy / m1; + u1z = p1fz / m1; + AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z)); + AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z)); + } + r = amrex::Random(); + if ( w1 > r*amrex::max(w1, w2) ) + { + u2x = p2fx / m2; + u2y = p2fy / m2; + u2z = p2fz / m2; + AMREX_ASSERT(!std::isnan(u1x+u1y+u1z+u2x+u2y+u2z)); + AMREX_ASSERT(!std::isinf(u1x+u1y+u1z+u2x+u2y+u2z)); + } + +} + +#endif // WARPX_PARTICLES_COLLISION_UPDATE_MOMENTUM_PEREZ_ELASTIC_H_ diff --git a/Source/Particles/Make.package b/Source/Particles/Make.package index 6d34fa0ef..8e336ba69 100644 --- a/Source/Particles/Make.package +++ b/Source/Particles/Make.package @@ -19,6 +19,7 @@ include $(WARPX_HOME)/Source/Particles/Deposition/Make.package include $(WARPX_HOME)/Source/Particles/Gather/Make.package include $(WARPX_HOME)/Source/Particles/Sorting/Make.package include $(WARPX_HOME)/Source/Particles/ParticleCreation/Make.package +include $(WARPX_HOME)/Source/Particles/Collision/Make.package INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Particles VPATH_LOCATIONS += $(WARPX_HOME)/Source/Particles diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 72064ce8c..9db129b05 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -16,19 +16,27 @@ #include <QuantumSyncEngineWrapper.H> #endif +#include "CollisionType.H" + #include <memory> #include <map> #include <string> #include <algorithm> -// -// MultiParticleContainer holds multiple (nspecies or npsecies+1 when -// using laser) WarpXParticleContainer. WarpXParticleContainer is -// derived from amrex::ParticleContainer, and it is -// polymorphic. PhysicalParticleContainer and LaserParticleContainer are -// concrete class derived from WarpXParticlContainer. -// - +/** + * The class MultiParticleContainer holds multiple instances of the polymorphic + * class WarpXParticleContainer, stored in its member variable "allcontainers". + * The class WarpX typically has a single (pointer to an) instance of + * MultiParticleContainer. + * + * MultiParticleContainer typically has two types of functions: + * - Functions that loop over all instances of WarpXParticleContainer in + * allcontainers and calls the corresponding function (for instance, + * MultiParticleContainer::Evolve loops over all particles containers and + * calls the corresponding WarpXParticleContainer::Evolve function). + * - Functions that specifically handle multiple species (for instance + * ReadParameters or mapSpeciesProduct). + */ class MultiParticleContainer { @@ -142,6 +150,8 @@ public: void doFieldIonization (); + void doCoulombCollisions (); + void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; @@ -214,6 +224,10 @@ protected: std::vector<std::string> lasers_names; + std::vector<std::string> collision_names; + + amrex::Vector<std::unique_ptr<CollisionType> > allcollisions; + //! instead of depositing (current, charge) on the finest patch level, deposit to the coarsest grid std::vector<bool> m_deposit_on_main_grid; @@ -295,5 +309,6 @@ private: // runtime parameters int nlasers = 0; int nspecies = 1; // physical particles only. nspecies+nlasers == allcontainers.size(). + int ncollisions = 0; }; #endif /*WARPX_ParticleContainer_H_*/ diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index bf1ffb74e..d84bc1afa 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -54,6 +54,14 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) } } ionization_process = IonizationProcess(); + + // collision + allcollisions.resize(ncollisions); + for (int i = 0; i < ncollisions; ++i) { + allcollisions[i].reset + (new CollisionType(species_names, collision_names[i])); + } + } void @@ -120,6 +128,15 @@ MultiParticleContainer::ReadParameters () } } + // collision + ParmParse pc("collisions"); + pc.query("ncollisions", ncollisions); + BL_ASSERT(ncollisions >= 0); + if (ncollisions > 0) { + pc.getarr("collision_names", collision_names); + BL_ASSERT(collision_names.size() == ncollisions); + } + } pp.query("use_fdtd_nci_corr", WarpX::use_fdtd_nci_corr); @@ -609,6 +626,40 @@ MultiParticleContainer::doFieldIonization () } // pc_source } +void +MultiParticleContainer::doCoulombCollisions () +{ + BL_PROFILE("MPC::doCoulombCollisions"); + + for (int i = 0; i < ncollisions; ++i) + { + auto& species1 = allcontainers[ allcollisions[i]->m_species1_index ]; + auto& species2 = allcontainers[ allcollisions[i]->m_species2_index ]; + + // Enable tiling + MFItInfo info; + if (Gpu::notInLaunchRegion()) info.EnableTiling(species1->tile_size); + + // Loop over refinement levels + for (int lev = 0; lev <= species1->finestLevel(); ++lev){ + + // Loop over all grids/tiles at this level +#ifdef _OPENMP + info.SetDynamic(true); + #pragma omp parallel if (Gpu::notInLaunchRegion()) +#endif + for (MFIter mfi = species1->MakeMFIter(lev, info); mfi.isValid(); ++mfi){ + + CollisionType::doCoulombCollisionsWithinTile + ( lev, mfi, species1, species2, + allcollisions[i]->m_isSameSpecies, + allcollisions[i]->m_CoulombLog ); + + } + } + } +} + #ifdef WARPX_QED void MultiParticleContainer::InitQED () { diff --git a/Source/Particles/ParticleCreation/TransformParticle.H b/Source/Particles/ParticleCreation/TransformParticle.H index 1b71df723..c0158db78 100644 --- a/Source/Particles/ParticleCreation/TransformParticle.H +++ b/Source/Particles/ParticleCreation/TransformParticle.H @@ -21,7 +21,7 @@ void transformSourceParticle( /** * \brief small modifications on product particle * \param i index of particle - * \param particle particle struct + * \param v_particle pointer to vector of particles */ template < ElementaryProcessType ProcessT > AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 580f13f86..9b688cc59 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -4,6 +4,13 @@ #include <PhysicalParticleContainer.H> #include <AMReX_Vector.H> +/** + * Photon particles have no mass, they deposit no charge, and see specific QED + * effects. For these reasons, they are stored in the separate particle + * container PhotonParticleContainer, that inherts from + * PhysicalParticleContainer. The particle pusher and current deposition, in + * particular, are overriden in this container. + */ class PhotonParticleContainer : public PhysicalParticleContainer { diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index e70b470b8..0192ffb37 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -15,6 +15,13 @@ #include <map> +/** + * PhysicalParticleContainer is the ParticleContainer class containing plasma + * particles (if a simulation has 2 plasma species, say "electrons" and + * "ions"), they will be two instances of PhysicalParticleContainer. + * + * PhysicalParticleContainer inherits from WarpXParticleContainer. + */ class PhysicalParticleContainer : public WarpXParticleContainer { @@ -68,6 +75,39 @@ public: int lev, int depos_lev); + /** + * \brief Evolve is the central function PhysicalParticleContainer that + * advances plasma particles for a time dt (typically one timestep). + * + * \param lev level on which particles are living + * \param Ex MultiFab from which field Ex is gathered + * \param Ey MultiFab from which field Ey is gathered + * \param Ez MultiFab from which field Ez is gathered + * \param Bx MultiFab from which field Bx is gathered + * \param By MultiFab from which field By is gathered + * \param Bz MultiFab from which field Bz is gathered + * \param jx MultiFab to which the particles' current jx is deposited + * \param jy MultiFab to which the particles' current jy is deposited + * \param jz MultiFab to which the particles' current jz is deposited + * \param cjx Same as jx (coarser, from lev-1), when using deposition buffers + * \param cjy Same as jy (coarser, from lev-1), when using deposition buffers + * \param cjz Same as jz (coarser, from lev-1), when using deposition buffers + * \param rho MultiFab to which the particles' charge is deposited + * \param crho Same as rho (coarser, from lev-1), when using deposition buffers + * \param cEx Same as Ex (coarser, from lev-1), when using gather buffers + * \param cEy Same as Ey (coarser, from lev-1), when using gather buffers + * \param cEz Same as Ez (coarser, from lev-1), when using gather buffers + * \param cBx Same as Bx (coarser, from lev-1), when using gather buffers + * \param cBy Same as By (coarser, from lev-1), when using gather buffers + * \param cBz Same as Bz (coarser, from lev-1), when using gather buffers + * \param t current physical time + * \param dt time step by which particles are advanced + * \param a_dt_type type of time step (used for sub-cycling) + * + * Evolve iterates over particle iterator (each box) and performs filtering, + * field gather, particle push and current deposition for all particles + * in the box. + */ virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -163,13 +203,39 @@ public: * \brief Apply NCI Godfrey filter to all components of E and B before gather * \param lev MR level * \param box box onto which the filter is applied - * \param exeli safeguard to avoid destructing arrays between ParIter iterations on GPU + * \param exeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Ex + * \param eyeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Ey + * \param ezeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Ez + * \param bxeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Bx + * \param byeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field By + * \param bzeli safeguard Elixir object (to avoid de-allocating too early + --between ParIter iterations-- on GPU) for field Bz * \param filtered_Ex Array containing filtered value + * \param filtered_Ey Array containing filtered value + * \param filtered_Ez Array containing filtered value + * \param filtered_Bx Array containing filtered value + * \param filtered_By Array containing filtered value + * \param filtered_Bz Array containing filtered value * \param Ex Field array before filtering (not modified) - * \param ex_ptr pointer to the array used for field gather. + * \param Ey Field array before filtering (not modified) + * \param Ez Field array before filtering (not modified) + * \param Bx Field array before filtering (not modified) + * \param By Field array before filtering (not modified) + * \param Bz Field array before filtering (not modified) + * \param exfab pointer to the Ex field (modified) + * \param eyfab pointer to the Ey field (modified) + * \param ezfab pointer to the Ez field (modified) + * \param bxfab pointer to the Bx field (modified) + * \param byfab pointer to the By field (modified) + * \param bzfab pointer to the Bz field (modified) * * The NCI Godfrey filter is applied on Ex, the result is stored in filtered_Ex - * and the pointer is modified (before this function is called, it points to Ex + * and the pointer exfab is modified (before this function is called, it points to Ex * and after this function is called, it points to Ex_filtered) */ void applyNCIFilter ( diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 95e8ce4fe..94d9bc363 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -963,8 +963,12 @@ PhysicalParticleContainer::EvolveES (const Vector<std::array<std::unique_ptr<Mul void PhysicalParticleContainer::FieldGather (int lev, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz) + const amrex::MultiFab& Ex, + const amrex::MultiFab& Ey, + const amrex::MultiFab& Ez, + const amrex::MultiFab& Bx, + const amrex::MultiFab& By, + const amrex::MultiFab& Bz) { const std::array<Real,3>& dx = WarpX::CellSize(lev); @@ -1078,11 +1082,11 @@ PhysicalParticleContainer::Evolve (int lev, for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const auto np = pti.numParticles(); - const auto lev = pti.GetLevel(); + const auto t_lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); tmp_particle_data.resize(finestLevel()+1); for (int i = 0; i < TmpIdx::nattribs; ++i) - tmp_particle_data[lev][index][i].resize(np); + tmp_particle_data[t_lev][index][i].resize(np); } } @@ -1427,15 +1431,19 @@ PhysicalParticleContainer::SplitParticles(int lev) { pti.GetPosition(xp, yp, zp); - // offset for split particles is computed as a function of cell size - // and number of particles per cell, so that a uniform distribution - // before splitting results in a uniform distribution after splitting const amrex::Vector<int> ppc_nd = plasma_injector->num_particles_per_cell_each_dim; const std::array<Real,3>& dx = WarpX::CellSize(lev); - amrex::Vector<Real> split_offset = {dx[0]/2._rt/ppc_nd[0], - dx[1]/2._rt/ppc_nd[1], - dx[2]/2._rt/ppc_nd[2]}; - + amrex::Vector<Real> split_offset = {dx[0]/2._rt, + dx[1]/2._rt, + dx[2]/2._rt}; + if (ppc_nd[0] > 0){ + // offset for split particles is computed as a function of cell size + // and number of particles per cell, so that a uniform distribution + // before splitting results in a uniform distribution after splitting + split_offset[0] /= ppc_nd[0]; + split_offset[1] /= ppc_nd[1]; + split_offset[2] /= ppc_nd[2]; + } // particle Array Of Structs data auto& particles = pti.GetArrayOfStructs(); // particle Struct Of Arrays data @@ -2148,12 +2156,12 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, RealVector& Bxp, RealVector& Byp, RealVector& Bzp, - FArrayBox const * exfab, - FArrayBox const * eyfab, - FArrayBox const * ezfab, - FArrayBox const * bxfab, - FArrayBox const * byfab, - FArrayBox const * bzfab, + amrex::FArrayBox const * exfab, + amrex::FArrayBox const * eyfab, + amrex::FArrayBox const * ezfab, + amrex::FArrayBox const * bxfab, + amrex::FArrayBox const * byfab, + amrex::FArrayBox const * bzfab, const int ngE, const int e_is_nodal, const long offset, const long np_to_gather, diff --git a/Source/Particles/RigidInjectedParticleContainer.H b/Source/Particles/RigidInjectedParticleContainer.H index a2473c5ad..fecb9c48e 100644 --- a/Source/Particles/RigidInjectedParticleContainer.H +++ b/Source/Particles/RigidInjectedParticleContainer.H @@ -4,6 +4,24 @@ #include <PhysicalParticleContainer.H> #include <AMReX_Vector.H> +/** + * When injecting a particle beam (typically for a plasma wakefield + * acceleration simulation), say propagating in the z direction, it can + * be necessary to make particles propagate in a straight line up to a given + * location z=z0. This is of particular importance when running in a boosted + * frame, where the beam may evolve due to its space charge fields before + * entering the plasma, causing the actual injected beam, and hence the whole + * simulation result, to depend on the Lorentz factor of the boost. + * + * This feature is implemented in RigidInjectedParticleContainer: At each + * iteration, for each particle, if z<z0 the particle moves in a straight line, + * and if z>z0 the particle evolves as a regular PhysicalParticleContainer. + * + * Note: This option is also useful to build self-consistent space charge + * fields for the particle beam. + * + * RigidInjectedParticleContainer derives from PhysicalParticleContainer. + */ class RigidInjectedParticleContainer : public PhysicalParticleContainer { diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 924c58382..fa5c586da 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -101,8 +101,32 @@ public: } }; +// Forward-declaration needed by WarpXParticleContainer below class MultiParticleContainer; +/** + * WarpXParticleContainer is the base polymorphic class from which all concrete + * particle container classes (that store a collection of particles) derive. Derived + * classes can be used for plasma particles, photon particles, or non-physical + * particles (e.g., for the laser antenna). + * It derives from amrex::ParticleContainer<0,0,PIdx::nattribs>, where the + * template arguments stand for the number of int and amrex::Real SoA and AoS + * data in amrex::Particle. + * - AoS amrex::Real: x, y, z (default), 0 additional (first template + * parameter) + * - AoS int: id, cpu (default), 0 additional (second template parameter) + * - SoA amrex::Real: PIdx::nattribs (third template parameter), see PIdx for + * the list. + * + * WarpXParticleContainer contains the main functions for initialization, + * interaction with the grid (field gather and current deposition) and particle + * push. + * + * Note: many functions are pure virtual (meaning they MUST be defined in + * derived classes, e.g., Evolve) or empty function (meaning they + * do not do anything, e.g., FieldGather, meant to be overriden by derived + * function) or actual functions (e.g. CurrentDeposition). + */ class WarpXParticleContainer : public amrex::ParticleContainer<0,0,PIdx::nattribs> { @@ -140,6 +164,11 @@ public: amrex::Real t, amrex::Real dt) = 0; #endif // WARPX_DO_ELECTROSTATIC + /** + * Evolve is the central WarpXParticleContainer function that advances + * particles for a time dt (typically one timestep). It is a pure virtual + * function for flexibility. + */ virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, @@ -288,6 +317,11 @@ public: std::map<std::string, int> getParticleComps () { return particle_comps;} std::map<std::string, int> getParticleiComps () { return particle_icomps;} + //amrex::Real getCharge () {return charge;} + amrex::ParticleReal getCharge () const {return charge;} + //amrex::Real getMass () {return mass;} + amrex::ParticleReal getMass () const {return mass;} + protected: std::map<std::string, int> particle_comps; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 9f02da338..e6c6d31bd 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -420,7 +420,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, void WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, const int * const ion_lev, - MultiFab* rho, int icomp, + amrex::MultiFab* rho, int icomp, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev) { @@ -507,7 +507,7 @@ WarpXParticleContainer::DepositCharge (WarpXParIter& pti, RealVector& wp, } void -WarpXParticleContainer::DepositCharge (Vector<std::unique_ptr<MultiFab> >& rho, +WarpXParticleContainer::DepositCharge (amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, bool local, bool reset, bool do_rz_volume_scaling) { @@ -760,7 +760,7 @@ WarpXParticleContainer::PushXES (Real dt) } void -WarpXParticleContainer::PushX (Real dt) +WarpXParticleContainer::PushX (amrex::Real dt) { const int nLevels = finestLevel(); for (int lev = 0; lev <= nLevels; ++lev) { @@ -769,7 +769,7 @@ WarpXParticleContainer::PushX (Real dt) } void -WarpXParticleContainer::PushX (int lev, Real dt) +WarpXParticleContainer::PushX (int lev, amrex::Real dt) { BL_PROFILE("WPC::PushX()"); @@ -846,8 +846,9 @@ WarpXParticleContainer::particlePostLocate(ParticleType& p, { p.m_idata.id = DoSplitParticleID; } - // For the moment, do not do anything if particles goes - // to lower level. + if (pld.m_lev == lev-1){ + // For the moment, do not do anything if particles goes + // to lower level. } } diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 31315f64f..54dfe1955 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -14,8 +14,9 @@ namespace *num_boxes = mf.local_size(); int shapesize = AMREX_SPACEDIM; if (mf.nComp() > 1) shapesize += 1; - *shapes = (int*) malloc(shapesize * (*num_boxes) * sizeof(int)); - amrex::Real** data = (amrex::Real**) malloc((*num_boxes) * sizeof(amrex::Real*)); + *shapes = static_cast<int*>(malloc(sizeof(int)*shapesize * (*num_boxes))); + auto data = + static_cast<amrex::Real**>(malloc((*num_boxes) * sizeof(amrex::Real*))); #ifdef _OPENMP #pragma omp parallel diff --git a/Source/QED/QedTableParserHelperFunctions.H b/Source/QED/QedTableParserHelperFunctions.H index 9f9f37017..528613727 100644 --- a/Source/QED/QedTableParserHelperFunctions.H +++ b/Source/QED/QedTableParserHelperFunctions.H @@ -15,7 +15,7 @@ namespace QedUtils{ * This function safely extracts an amrex::Vector<T> from raw binary data. * T must be a simple datatype (e.g. an int, a float, a double...). * - * @param[in] p_char a pointer to the binary stream + * @param[in] p_data a pointer to the binary stream * @param[in] how_many how many T should be read from stream * @param[in] p_last a pointer to the last element of the char* array * @return {a tuple containing @@ -43,7 +43,7 @@ namespace QedUtils{ * This function safely extracts a T from raw binary data. * T must be a simple datatype (e.g. an int, a float, a double...). * - * @param[in] p_char a pointer to the binary stream + * @param[in] p_data a pointer to the binary stream * @param[in] p_last a pointer to the last element of the char* array * @return {a tuple containing * 1) flag (which is false if p_last is exceeded) diff --git a/Source/Utils/WarpXConst.H b/Source/Utils/WarpXConst.H index 7734d9011..70436cb72 100644 --- a/Source/Utils/WarpXConst.H +++ b/Source/Utils/WarpXConst.H @@ -13,7 +13,7 @@ namespace MathConst // https://physics.nist.gov/cuu/Constants/index.html namespace PhysConst { - static constexpr amrex::Real c = 299792458.; + static constexpr amrex::Real c = 299'792'458.; static constexpr amrex::Real ep0 = 8.8541878128e-12; static constexpr amrex::Real mu0 = 1.25663706212e-06; static constexpr amrex::Real q_e = 1.602176634e-19; diff --git a/Source/Utils/WarpXUtil.H b/Source/Utils/WarpXUtil.H index ab28c5446..195e309cf 100644 --- a/Source/Utils/WarpXUtil.H +++ b/Source/Utils/WarpXUtil.H @@ -1,3 +1,6 @@ +#ifndef WARPX_UTILS_H_ +#define WARPX_UTILS_H_ + #include <AMReX_REAL.H> #include <AMReX_Vector.H> #include <AMReX_MultiFab.H> @@ -21,3 +24,5 @@ namespace WarpXUtilIO{ */ bool WriteBinaryDataOnFile(std::string filename, const amrex::Vector<char>& data); } + +#endif //WARPX_UTILS_H_ diff --git a/Source/WarpX.H b/Source/WarpX.H index caffa5e25..fe8b2840d 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -171,7 +171,7 @@ public: void ComputeDt (); // Compute max_step automatically for simulations in a boosted frame. - void computeMaxStepBoostAccelerator(amrex::Geometry geom); + void computeMaxStepBoostAccelerator(const amrex::Geometry& geom); int MoveWindow (bool move_j); void UpdatePlasmaInjectionPosition (amrex::Real dt); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c1eda30ba..377d103d1 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1116,9 +1116,9 @@ WarpX::Evolve (int numsteps) { } void -WarpX::ComputeDivB (MultiFab& divB, int dcomp, - const std::array<const MultiFab*, 3>& B, - const std::array<Real,3>& dx) +WarpX::ComputeDivB (amrex::MultiFab& divB, int dcomp, + const std::array<const amrex::MultiFab*, 3>& B, + const std::array<amrex::Real,3>& dx) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1150,9 +1150,9 @@ WarpX::ComputeDivB (MultiFab& divB, int dcomp, } void -WarpX::ComputeDivB (MultiFab& divB, int dcomp, - const std::array<const MultiFab*, 3>& B, - const std::array<Real,3>& dx, int ngrow) +WarpX::ComputeDivB (amrex::MultiFab& divB, int dcomp, + const std::array<const amrex::MultiFab*, 3>& B, + const std::array<amrex::Real,3>& dx, int ngrow) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1184,9 +1184,9 @@ WarpX::ComputeDivB (MultiFab& divB, int dcomp, } void -WarpX::ComputeDivE (MultiFab& divE, int dcomp, - const std::array<const MultiFab*, 3>& E, - const std::array<Real,3>& dx) +WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp, + const std::array<const amrex::MultiFab*, 3>& E, + const std::array<amrex::Real,3>& dx) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; @@ -1218,9 +1218,9 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, } void -WarpX::ComputeDivE (MultiFab& divE, int dcomp, - const std::array<const MultiFab*, 3>& E, - const std::array<Real,3>& dx, int ngrow) +WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp, + const std::array<const amrex::MultiFab*, 3>& E, + const std::array<amrex::Real,3>& dx, int ngrow) { Real dxinv = 1./dx[0], dyinv = 1./dx[1], dzinv = 1./dx[2]; diff --git a/Tools/read_raw_data.py b/Tools/read_raw_data.py index 756a5a354..c907b1fe1 100644 --- a/Tools/read_raw_data.py +++ b/Tools/read_raw_data.py @@ -297,4 +297,4 @@ def _read_buffer(snapshot, header_fn, _component_names): if __name__ == "__main__": - data = read_lab_snapshot("lab_frame_data/snapshot00012"); + data = read_lab_snapshot("lab_frame_data/snapshot00012", "lab_frame_data/Header"); |