diff options
34 files changed, 1232 insertions, 215 deletions
diff --git a/Docs/requirements.txt b/Docs/requirements.txt index 478b0ba60..68c49affd 100644 --- a/Docs/requirements.txt +++ b/Docs/requirements.txt @@ -1,4 +1,6 @@ sphinx_rtd_theme>=0.3.1 recommonmark -sphinx +sphinx>=2.0 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 2e6acfdef..597d36996 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -282,12 +282,13 @@ Particle initialization temperature parameter ``<species_name>.theta`` as an input, where theta is kb*T/(m*c^2), kb is the Boltzmann constant, c is the speed of light, and m is the mass of the species. It also includes the optional parameter ``<species_name>.beta`` where beta is equal to v/c. - The plasma will be initialized to move at drift velocity beta*c in the positive - ``<species_name>.direction = 'x', 'y', 'z'`` direction. The MB distribution is initialized - in the drifting frame by sampling three Gaussian distributions in each dimension using, - the Box Mueller method, and then the distribution is transformed to the simulation frame - using the flipping method. The flipping method can be found in Zenitani 2015 - section III. B. (Phys. Plasmas 22, 042116). + The plasma will be initialized to move at drift velocity beta*c in the + ``<species_name>.drift_vel_dir = (+/-) 'x', 'y', 'z'`` direction. Please leave no whitespace + between the sign and the character on input. A direction without a sign will be treated as + positive. The MB distribution is initialized in the drifting frame by sampling three Gaussian + distributions in each dimension using, the Box Mueller method, and then the distribution is + transformed to the simulation frame using the flipping method. The flipping method can be + found in Zenitani 2015 section III. B. (Phys. Plasmas 22, 042116). Note that though the particles may move at relativistic speeds in the simulation frame, they are not relativistic in the drift frame. This is as opposed to the Maxwell Juttner @@ -297,11 +298,12 @@ Particle initialization requires a dimensionless temperature parameter ``<species_name>.theta``, where theta is equal to kb*T/(m*c^2), where kb is the Boltzmann constant, and m is the mass of the species. It also includes the optional parameter ``<species_name>.beta`` where beta is equal to v/c. The plasma - will be initialized to move at velocity beta*c in the positive - ``<species_name>.direction = 'x', 'y', 'z'`` direction. The MJ distribution will be initialized - in the moving frame using the Sobol method, and then the distribution will be transformed to the - simulation frame using the flipping method. Both the Sobol and the flipping method can be found - in Zenitani 2015 (Phys. Plasmas 22, 042116). + will be initialized to move at velocity beta*c in the + ``<species_name>.drift_vel_dir = (+/-) 'x', 'y', 'z'`` direction. Please leave no whitespace + between the sign and the character on input. A direction without a sign will be treated as + positive. The MJ distribution will be initialized in the moving frame using the Sobol method, + and then the distribution will be transformed to the simulation frame using the flipping method. + Both the Sobol and the flipping method can be found in Zenitani 2015 (Phys. Plasmas 22, 042116). Please take notice that particles initialized with this setting can be relativistic in two ways. In the simulation frame, they can drift with a relativistic speed beta. Then, in the drifting diff --git a/Python/pywarpx/_libwarpx.py b/Python/pywarpx/_libwarpx.py index 6c3e54619..0444140e4 100755 --- a/Python/pywarpx/_libwarpx.py +++ b/Python/pywarpx/_libwarpx.py @@ -130,18 +130,30 @@ libwarpx.warpx_getEfieldCP.restype = _LP_LP_c_real libwarpx.warpx_getEfieldCPLoVects.restype = _LP_c_int libwarpx.warpx_getEfieldFP.restype = _LP_LP_c_real libwarpx.warpx_getEfieldFPLoVects.restype = _LP_c_int +libwarpx.warpx_getEfieldCP_PML.restype = _LP_LP_c_real +libwarpx.warpx_getEfieldCPLoVects_PML.restype = _LP_c_int +libwarpx.warpx_getEfieldFP_PML.restype = _LP_LP_c_real +libwarpx.warpx_getEfieldFPLoVects_PML.restype = _LP_c_int libwarpx.warpx_getBfield.restype = _LP_LP_c_real libwarpx.warpx_getBfieldLoVects.restype = _LP_c_int libwarpx.warpx_getBfieldCP.restype = _LP_LP_c_real libwarpx.warpx_getBfieldCPLoVects.restype = _LP_c_int libwarpx.warpx_getBfieldFP.restype = _LP_LP_c_real libwarpx.warpx_getBfieldFPLoVects.restype = _LP_c_int +libwarpx.warpx_getBfieldCP_PML.restype = _LP_LP_c_real +libwarpx.warpx_getBfieldCPLoVects_PML.restype = _LP_c_int +libwarpx.warpx_getBfieldFP_PML.restype = _LP_LP_c_real +libwarpx.warpx_getBfieldFPLoVects_PML.restype = _LP_c_int libwarpx.warpx_getCurrentDensity.restype = _LP_LP_c_real libwarpx.warpx_getCurrentDensityLoVects.restype = _LP_c_int libwarpx.warpx_getCurrentDensityCP.restype = _LP_LP_c_real libwarpx.warpx_getCurrentDensityCPLoVects.restype = _LP_c_int libwarpx.warpx_getCurrentDensityFP.restype = _LP_LP_c_real libwarpx.warpx_getCurrentDensityFPLoVects.restype = _LP_c_int +libwarpx.warpx_getCurrentDensityCP_PML.restype = _LP_LP_c_real +libwarpx.warpx_getCurrentDensityCPLoVects_PML.restype = _LP_c_int +libwarpx.warpx_getCurrentDensityFP_PML.restype = _LP_LP_c_real +libwarpx.warpx_getCurrentDensityFPLoVects_PML.restype = _LP_c_int #libwarpx.warpx_getPMLSigma.restype = _LP_c_real #libwarpx.warpx_getPMLSigmaStar.restype = _LP_c_real @@ -747,6 +759,66 @@ def get_mesh_electric_field_fp(level, direction, include_ghosts=True): return _get_mesh_field_list(libwarpx.warpx_getEfieldFP, level, direction, include_ghosts) +def get_mesh_electric_field_cp_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of numpy arrays containing the mesh electric field + data on each grid for this process. This version returns the field on + the coarse patch for the PML for the given level. + + The data for the numpy arrays are not copied, but share the underlying + memory buffer with WarpX. The numpy arrays are fully writeable. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A List of numpy arrays. + + ''' + + try: + return _get_mesh_field_list(libwarpx.warpx_getEfieldCP_PML, level, direction, include_ghosts) + except ValueError: + raise Exception('PML not initialized') + + +def get_mesh_electric_field_fp_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of numpy arrays containing the mesh electric field + data on each grid for this process. This version returns the field on + the fine patch for the PML for the given level. + + The data for the numpy arrays are not copied, but share the underlying + memory buffer with WarpX. The numpy arrays are fully writeable. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A List of numpy arrays. + + ''' + + try: + return _get_mesh_field_list(libwarpx.warpx_getEfieldFP_PML, level, direction, include_ghosts) + except ValueError: + raise Exception('PML not initialized') + + def get_mesh_magnetic_field(level, direction, include_ghosts=True): ''' @@ -829,6 +901,66 @@ def get_mesh_magnetic_field_fp(level, direction, include_ghosts=True): return _get_mesh_field_list(libwarpx.warpx_getBfieldFP, level, direction, include_ghosts) +def get_mesh_magnetic_field_cp_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of numpy arrays containing the mesh magnetic field + data on each grid for this process. This version returns the field on + the coarse patch for the PML for the given level. + + The data for the numpy arrays are not copied, but share the underlying + memory buffer with WarpX. The numpy arrays are fully writeable. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A List of numpy arrays. + + ''' + + try: + return _get_mesh_field_list(libwarpx.warpx_getBfieldCP_PML, level, direction, include_ghosts) + except ValueError: + raise Exception('PML not initialized') + + +def get_mesh_magnetic_field_fp_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of numpy arrays containing the mesh magnetic field + data on each grid for this process. This version returns the field on + the fine patch for the PML for the given level. + + The data for the numpy arrays are not copied, but share the underlying + memory buffer with WarpX. The numpy arrays are fully writeable. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A List of numpy arrays. + + ''' + + try: + return _get_mesh_field_list(libwarpx.warpx_getBfieldFP_PML, level, direction, include_ghosts) + except ValueError: + raise Exception('PML not initialized') + + def get_mesh_current_density(level, direction, include_ghosts=True): ''' @@ -909,6 +1041,66 @@ def get_mesh_current_density_fp(level, direction, include_ghosts=True): return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityFP, level, direction, include_ghosts) +def get_mesh_current_density_cp_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of numpy arrays containing the mesh current density + data on each grid for this process. This version returns the density for + the coarse patch for the PML for the given level. + + The data for the numpy arrays are not copied, but share the underlying + memory buffer with WarpX. The numpy arrays are fully writeable. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A List of numpy arrays. + + ''' + + try: + return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityCP_PML, level, direction, include_ghosts) + except ValueError: + raise Exception('PML not initialized') + + +def get_mesh_current_density_fp_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of numpy arrays containing the mesh current density + data on each grid for this process. This version returns the density on + the fine patch for the PML for the given level. + + The data for the numpy arrays are not copied, but share the underlying + memory buffer with WarpX. The numpy arrays are fully writeable. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A List of numpy arrays. + + ''' + + try: + return _get_mesh_field_list(libwarpx.warpx_getCurrentDensityFP_PML, level, direction, include_ghosts) + except ValueError: + raise Exception('PML not initialized') + + def _get_mesh_array_lovects(level, direction, include_ghosts=True, getarrayfunc=None): assert(0 <= level and level <= libwarpx.warpx_finestLevel()) @@ -998,6 +1190,56 @@ def get_mesh_electric_field_fp_lovects(level, direction, include_ghosts=True): return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getEfieldFPLoVects) +def get_mesh_electric_field_cp_lovects_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of the lo vectors of the arrays containing the mesh electric field + data on each PML grid for this process. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids) + + ''' + try: + return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getEfieldCPLoVects_PML) + except ValueError: + raise Exception('PML not initialized') + + +def get_mesh_electric_field_fp_lovects_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of the lo vectors of the arrays containing the mesh electric field + data on each PML grid for this process. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids) + + ''' + try: + return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getEfieldFPLoVects_PML) + except ValueError: + raise Exception('PML not initialized') + + def get_mesh_magnetic_field_lovects(level, direction, include_ghosts=True): ''' @@ -1066,6 +1308,56 @@ def get_mesh_magnetic_field_fp_lovects(level, direction, include_ghosts=True): return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getBfieldFPLoVects) +def get_mesh_magnetic_field_cp_lovects_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of the lo vectors of the arrays containing the mesh electric field + data on each PML grid for this process. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids) + + ''' + try: + return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getBfieldCPLoVects_PML) + except ValueError: + raise Exception('PML not initialized') + + +def get_mesh_magnetic_field_fp_lovects_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of the lo vectors of the arrays containing the mesh electric field + data on each PML grid for this process. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids) + + ''' + try: + return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getBfieldFPLoVects_PML) + except ValueError: + raise Exception('PML not initialized') + + def get_mesh_current_density_lovects(level, direction, include_ghosts=True): ''' @@ -1129,3 +1421,52 @@ def get_mesh_current_density_fp_lovects(level, direction, include_ghosts=True): ''' return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getCurrentDensityFPLoVects) + + +def get_mesh_current_density_cp_lovects_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of the lo vectors of the arrays containing the mesh electric field + data on each PML grid for this process. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids) + + ''' + try: + return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getCurrentDensityCPLoVects_PML) + except ValueError: + raise Exception('PML not initialized') + +def get_mesh_current_density_fp_lovects_pml(level, direction, include_ghosts=True): + ''' + + This returns a list of the lo vectors of the arrays containing the mesh electric field + data on each PML grid for this process. + + Parameters + ---------- + + level : the AMR level to get the data for + direction : the component of the data you want + include_ghosts : whether to include ghost zones or not + + Returns + ------- + + A 2d numpy array of the lo vector for each grid with the shape (dims, number of grids) + + ''' + try: + return _get_mesh_array_lovects(level, direction, include_ghosts, libwarpx.warpx_getCurrentDensityFPLoVects_PML) + except ValueError: + raise Exception('PML not initialized') diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py index fa247db07..920c109d6 100644 --- a/Python/pywarpx/fields.py +++ b/Python/pywarpx/fields.py @@ -87,14 +87,28 @@ class _MultiFABWrapper(object): def _getitem3d(self, index): """Returns slices of a 3D decomposed array, """ - ix = index[0] - iy = index[1] - iz = index[2] lovects = self._getlovects() hivects = self._gethivects() fields = self._getfields() + ix = index[0] + iy = index[1] + iz = index[2] + + if len(fields[0].shape) > self.dim: + ncomps = fields[0].shape[-1] + else: + ncomps = 1 + + if len(index) > self.dim: + if ncomps > 1: + ic = index[-1] + else: + raise Exception('Too many indices given') + else: + ic = None + nx = hivects[0,:].max() - self.nghosts ny = hivects[1,:].max() - self.nghosts nz = hivects[2,:].max() - self.nghosts @@ -124,9 +138,12 @@ class _MultiFABWrapper(object): izstop = iz + 1 # --- Setup the size of the array to be returned and create it. + # --- Space is added for multiple components if needed. sss = (max(0, ixstop - ixstart), max(0, iystop - iystart), max(0, izstop - izstart)) + if ncomps > 1 and ic is None: + sss = tuple(list(sss) + [ncomps]) resultglobal = np.zeros(sss, dtype=_libwarpx._numpy_real_dtype) datalist = [] @@ -145,6 +162,8 @@ class _MultiFABWrapper(object): sss = (slice(ix1 - lovects[0,i], ix2 - lovects[0,i]), slice(iy1 - lovects[1,i], iy2 - lovects[1,i]), slice(iz1 - lovects[2,i], iz2 - lovects[2,i])) + if ic is not None: + sss = tuple(list(sss) + [ic]) vslice = (slice(ix1 - ixstart, ix2 - ixstart), slice(iy1 - iystart, iy2 - iystart), @@ -295,6 +314,14 @@ class _MultiFABWrapper(object): hivects = self._gethivects() fields = self._getfields() + if len(index) > self.dim: + if ncomps > 1: + ic = index[-1] + else: + raise Exception('Too many indices given') + else: + ic = None + nx = hivects[0,:].max() - self.nghosts ny = hivects[1,:].max() - self.nghosts nz = hivects[2,:].max() - self.nghosts @@ -343,6 +370,8 @@ class _MultiFABWrapper(object): sss = (slice(ix1 - lovects[0,i], ix2 - lovects[0,i]), slice(iy1 - lovects[1,i], iy2 - lovects[1,i]), slice(iz1 - lovects[2,i], iz2 - lovects[2,i])) + if ic is not None: + sss = tuple(list(sss) + [ic]) if isinstance(value, np.ndarray): vslice = (slice(ix1 - ixstart, ix2 - ixstart), @@ -593,3 +622,119 @@ def JzFPWrapper(level=0, include_ghosts=False): get_lovects=_libwarpx.get_mesh_current_density_fp_lovects, get_fabs=_libwarpx.get_mesh_current_density_fp, level=level, include_ghosts=include_ghosts) +def ExCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=0, overlaps=[0,1,1], + get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_electric_field_cp_pml, + level=level, include_ghosts=include_ghosts) + +def EyCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=1, overlaps=[1,0,1], + get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_electric_field_cp_pml, + level=level, include_ghosts=include_ghosts) + +def EzCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=2, overlaps=[1,1,0], + get_lovects=_libwarpx.get_mesh_electric_field_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_electric_field_cp_pml, + level=level, include_ghosts=include_ghosts) + +def BxCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=0, overlaps=[1,0,0], + get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml, + level=level, include_ghosts=include_ghosts) + +def ByCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=1, overlaps=[0,1,0], + get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml, + level=level, include_ghosts=include_ghosts) + +def BzCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=2, overlaps=[0,0,1], + get_lovects=_libwarpx.get_mesh_magnetic_field_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_magnetic_field_cp_pml, + level=level, include_ghosts=include_ghosts) + +def JxCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=0, overlaps=[0,1,1], + get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_current_density_cp_pml, + level=level, include_ghosts=include_ghosts) + +def JyCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=1, overlaps=[1,0,1], + get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_current_density_cp_pml, + level=level, include_ghosts=include_ghosts) + +def JzCPPMLWrapper(level=1, include_ghosts=False): + assert level>0, Exception('Coarse patch only available on levels > 0') + return _MultiFABWrapper(direction=2, overlaps=[1,1,0], + get_lovects=_libwarpx.get_mesh_current_density_cp_lovects_pml, + get_fabs=_libwarpx.get_mesh_current_density_cp_pml, + level=level, include_ghosts=include_ghosts) + +def ExFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=0, overlaps=[0,1,1], + get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_electric_field_fp_pml, + level=level, include_ghosts=include_ghosts) + +def EyFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=1, overlaps=[1,0,1], + get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_electric_field_fp_pml, + level=level, include_ghosts=include_ghosts) + +def EzFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=2, overlaps=[1,1,0], + get_lovects=_libwarpx.get_mesh_electric_field_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_electric_field_fp_pml, + level=level, include_ghosts=include_ghosts) + +def BxFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=0, overlaps=[1,0,0], + get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml, + level=level, include_ghosts=include_ghosts) + +def ByFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=1, overlaps=[0,1,0], + get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml, + level=level, include_ghosts=include_ghosts) + +def BzFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=2, overlaps=[0,0,1], + get_lovects=_libwarpx.get_mesh_magnetic_field_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_magnetic_field_fp_pml, + level=level, include_ghosts=include_ghosts) + +def JxFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=0, overlaps=[0,1,1], + get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_current_density_fp_pml, + level=level, include_ghosts=include_ghosts) + +def JyFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=1, overlaps=[1,0,1], + get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_current_density_fp_pml, + level=level, include_ghosts=include_ghosts) + +def JzFPPMLWrapper(level=0, include_ghosts=False): + return _MultiFABWrapper(direction=2, overlaps=[1,1,0], + get_lovects=_libwarpx.get_mesh_current_density_fp_lovects_pml, + get_fabs=_libwarpx.get_mesh_current_density_fp_pml, + level=level, include_ghosts=include_ghosts) diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index 6e9c9153b..224a53737 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -75,6 +75,13 @@ class Species(picmistandard.PICMI_Species): if self.initial_distribution is not None: self.initial_distribution.initialize_inputs(self.species_number, layout, self.species, self.density_scale) + for interaction in self.interactions: + assert interaction[0] == 'ionization' + assert interaction[1] == 'ADK', 'WarpX only has ADK ionization model implemented' + self.species.do_field_ionization=1 + self.species.physical_element=self.particle_type + self.species.ionization_product_species = interaction[2].name + self.species.ionization_initial_level = self.charge_state picmistandard.PICMI_MultiSpecies.Species_class = Species class MultiSpecies(picmistandard.PICMI_MultiSpecies): @@ -326,13 +333,14 @@ class CylindricalGrid(picmistandard.PICMI_CylindricalGrid): pywarpx.geometry.prob_hi = self.upper_bound pywarpx.warpx.n_rz_azimuthal_modes = self.n_azimuthal_modes - if self.moving_window_velocity is not None and np.any(np.not_equal(self.moving_window_velocity, 0.)): - pywarpx.warpx.do_moving_window = 1 - if self.moving_window_velocity[0] != 0.: - raise Exception('In cylindrical coordinates, a moving window in r can not be done') - if self.moving_window_velocity[1] != 0.: - pywarpx.warpx.moving_window_dir = 'z' - pywarpx.warpx.moving_window_v = self.moving_window_velocity[1]/constants.c # in units of the speed of light + if self.moving_window_zvelocity is not None: + if np.isscalar(self.moving_window_zvelocity): + if self.moving_window_zvelocity !=0: + pywarpx.warpx.do_moving_window = 1 + pywarpx.warpx.moving_window_dir = 'z' + pywarpx.warpx.moving_window_v = self.moving_window_zvelocity/constants.c # in units of the speed of light + else: + raise Exception('RZ PICMI moving_window_velocity (only available in z direction) should be a scalar') if self.refined_regions: assert len(self.refined_regions) == 1, Exception('WarpX only supports one refined region.') diff --git a/Source/Diagnostics/ElectrostaticIO.cpp b/Source/Diagnostics/ElectrostaticIO.cpp index a023da0b7..332638cff 100644 --- a/Source/Diagnostics/ElectrostaticIO.cpp +++ b/Source/Diagnostics/ElectrostaticIO.cpp @@ -1,6 +1,3 @@ -#include <AMReX_MGT_Solver.H> -#include <AMReX_stencil_types.H> - #include <WarpX.H> #include <WarpX_f.H> @@ -98,24 +95,7 @@ WritePlotFileES (const amrex::Vector<std::unique_ptr<amrex::MultiFab> >& rho, } } - Vector<std::string> particle_varnames; - particle_varnames.push_back("weight"); - - particle_varnames.push_back("momentum_x"); - particle_varnames.push_back("momentum_y"); - particle_varnames.push_back("momentum_z"); - - particle_varnames.push_back("Ex"); - particle_varnames.push_back("Ey"); - particle_varnames.push_back("Ez"); - - particle_varnames.push_back("Bx"); - particle_varnames.push_back("By"); - particle_varnames.push_back("Bz"); - - Vector<std::string> int_names; - - mypc->Checkpoint(plotfilename, particle_varnames, int_names); + mypc->Checkpoint(plotfilename); WriteJobInfo(plotfilename); diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package index 7fd5e5f76..710e4399b 100644 --- a/Source/Diagnostics/Make.package +++ b/Source/Diagnostics/Make.package @@ -2,14 +2,17 @@ CEXE_sources += WarpXIO.cpp CEXE_sources += BackTransformedDiagnostic.cpp CEXE_sources += ParticleIO.cpp CEXE_sources += FieldIO.cpp +CEXE_sources += SliceDiagnostic.cpp +ifeq ($(DO_ELECTROSTATIC),TRUE) + CEXE_sources += ElectrostaticIO.cpp +endif + CEXE_headers += FieldIO.H CEXE_headers += BackTransformedDiagnostic.H -CEXE_headers += ElectrostaticIO.cpp CEXE_headers += SliceDiagnostic.H -CEXE_sources += SliceDiagnostic.cpp ifeq ($(USE_OPENPMD), TRUE) -CEXE_sources += WarpXOpenPMD.cpp + CEXE_sources += WarpXOpenPMD.cpp endif INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics diff --git a/Source/Evolve/WarpXEvolveES.cpp b/Source/Evolve/WarpXEvolveES.cpp index 61a8a5f32..555ab37ad 100644 --- a/Source/Evolve/WarpXEvolveES.cpp +++ b/Source/Evolve/WarpXEvolveES.cpp @@ -1,6 +1,3 @@ -#include <AMReX_MGT_Solver.H> -#include <AMReX_stencil_types.H> - #include <WarpX.H> #include <WarpX_f.H> @@ -11,16 +8,6 @@ namespace using namespace amrex; -class NoOpPhysBC - : public amrex::PhysBCFunctBase -{ -public: - NoOpPhysBC () {} - virtual ~NoOpPhysBC () {} - virtual void FillBoundary (amrex::MultiFab& mf, int, int, amrex::Real time) override { } - using amrex::PhysBCFunctBase::FillBoundary; -}; - void WarpX::EvolveES (int numsteps) { diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 5464430cf..f7c7e498f 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -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") { @@ -275,17 +276,26 @@ void PlasmaInjector::parseMomentum (ParmParse& pp) int dir = 0; std::string direction = "x"; pp.query("beta", beta); + if(beta < 0){ + amrex::Abort("Please enter a positive beta value. Drift direction is set with <s_name>.bulk_vel_dir = 'x' or '+x', '-x', 'y' or '+y', etc."); + } pp.query("theta", theta); - pp.query("direction", direction); - if(direction == "x" || direction == "X"){ + pp.query("bulk_vel_dir", direction); + if(direction[0] == '-'){ + beta = -beta; + } + if((direction == "x" || direction[1] == 'x') || + (direction == "X" || direction[1] == 'X')){ dir = 0; - } else if (direction == "y" || direction == "Y"){ + } else if ((direction == "y" || direction[1] == 'y') || + (direction == "Y" || direction[1] == 'Y')){ dir = 1; - } else if (direction == "z" || direction == "Z"){ + } else if ((direction == "z" || direction[1] == 'z') || + (direction == "Z" || direction[1] == 'Z')){ dir = 2; } else{ std::stringstream stringstream; - stringstream << "Direction " << direction << " is not recognzied. Please enter x, y, or z."; + stringstream << "Cannot interpret <s_name>.bulk_vel_dir input '" << direction << "'. Please enter +/- x, y, or z with no whitespace between the sign and other character."; direction = stringstream.str(); amrex::Abort(direction.c_str()); } @@ -297,17 +307,26 @@ void PlasmaInjector::parseMomentum (ParmParse& pp) int dir = 0; std::string direction = "x"; pp.query("beta", beta); + if(beta < 0){ + amrex::Abort("Please enter a positive beta value. Drift direction is set with <s_name>.bulk_vel_dir = 'x' or '+x', '-x', 'y' or '+y', etc."); + } pp.query("theta", theta); - pp.query("direction", direction); - if(direction == "x" || direction == "X"){ + pp.query("bulk_vel_dir", direction); + if(direction[0] == '-'){ + beta = -beta; + } + if((direction == "x" || direction[1] == 'x') || + (direction == "X" || direction[1] == 'X')){ dir = 0; - } else if (direction == "y" || direction == "Y"){ + } else if ((direction == "y" || direction[1] == 'y') || + (direction == "Y" || direction[1] == 'Y')){ dir = 1; - } else if (direction == "z" || direction == "Z"){ + } else if ((direction == "z" || direction[1] == 'z') || + (direction == "Z" || direction[1] == 'Z')){ dir = 2; } else{ std::stringstream stringstream; - stringstream << "Direction " << direction << " is not recognzied. Please enter x, y, or z."; + stringstream << "Cannot interpret <s_name>.bulk_vel_dir input '" << direction << "'. Please enter +/- x, y, or z with no whitespace between the sign and other character."; direction = stringstream.str(); amrex::Abort(direction.c_str()); } 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/Make.WarpX b/Source/Make.WarpX index 462299db9..fa708a62a 100644 --- a/Source/Make.WarpX +++ b/Source/Make.WarpX @@ -172,9 +172,6 @@ ifeq ($(USE_RZ),TRUE) endif ifeq ($(DO_ELECTROSTATIC),TRUE) - include $(AMREX_HOME)/Src/LinearSolvers/C_to_F_MG/Make.package - include $(AMREX_HOME)/Src/LinearSolvers/F_MG/FParallelMG.mak - include $(AMREX_HOME)/Src/F_BaseLib/FParallelMG.mak DEFINES += -DWARPX_DO_ELECTROSTATIC endif diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index d7ddc7a72..9db129b05 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -23,14 +23,20 @@ #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 { 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 1cd4ba621..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, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d5c08074a..94d9bc363 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1431,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 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 567bb402d..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, diff --git a/Source/Python/WarpXWrappers.cpp b/Source/Python/WarpXWrappers.cpp index 3074b1990..54dfe1955 100644 --- a/Source/Python/WarpXWrappers.cpp +++ b/Source/Python/WarpXWrappers.cpp @@ -208,113 +208,80 @@ extern "C" return myspc.TotalNumberOfParticles(); } - amrex::Real** warpx_getEfield(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getEfield(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getEfieldLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getEfield(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getEfieldCP(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getEfieldCPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getEfield_cp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getEfieldFP(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getEfieldFPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getEfield_fp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getBfield(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getBfield(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getBfieldLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getBfield(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getBfieldCP(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getBfieldCPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getBfield_cp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getBfieldFP(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getBfieldFPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getBfield_fp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getCurrentDensity(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getcurrent(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getCurrentDensityLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getcurrent(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getCurrentDensityCP(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getCurrentDensityCPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getcurrent_cp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } - - amrex::Real** warpx_getCurrentDensityFP(int lev, int direction, - int *return_size, int *ncomps, int *ngrow, int **shapes) { - auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); - return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); - } - - int* warpx_getCurrentDensityFPLoVects(int lev, int direction, - int *return_size, int *ngrow) { - auto & mf = WarpX::GetInstance().getcurrent_fp(lev, direction); - return getMultiFabLoVects(mf, return_size, ngrow); - } +#define WARPX_GET_FIELD(FIELD, GETTER) \ + amrex::Real** FIELD(int lev, int direction, \ + int *return_size, int *ncomps, int *ngrow, int **shapes) { \ + auto & mf = GETTER(lev, direction); \ + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); \ + } + +#define WARPX_GET_LOVECTS(FIELD, GETTER) \ + int* FIELD(int lev, int direction, \ + int *return_size, int *ngrow) { \ + auto & mf = GETTER(lev, direction); \ + return getMultiFabLoVects(mf, return_size, ngrow); \ + } + + WARPX_GET_FIELD(warpx_getEfield, WarpX::GetInstance().getEfield); + WARPX_GET_FIELD(warpx_getEfieldCP, WarpX::GetInstance().getEfield_cp); + WARPX_GET_FIELD(warpx_getEfieldFP, WarpX::GetInstance().getEfield_fp); + + WARPX_GET_FIELD(warpx_getBfield, WarpX::GetInstance().getBfield); + WARPX_GET_FIELD(warpx_getBfieldCP, WarpX::GetInstance().getBfield_cp); + WARPX_GET_FIELD(warpx_getBfieldFP, WarpX::GetInstance().getBfield_fp); + + WARPX_GET_FIELD(warpx_getCurrentDensity, WarpX::GetInstance().getcurrent); + WARPX_GET_FIELD(warpx_getCurrentDensityCP, WarpX::GetInstance().getcurrent_cp); + WARPX_GET_FIELD(warpx_getCurrentDensityFP, WarpX::GetInstance().getcurrent_fp); + + WARPX_GET_LOVECTS(warpx_getEfieldLoVects, WarpX::GetInstance().getEfield); + WARPX_GET_LOVECTS(warpx_getEfieldCPLoVects, WarpX::GetInstance().getEfield_cp); + WARPX_GET_LOVECTS(warpx_getEfieldFPLoVects, WarpX::GetInstance().getEfield_fp); + + WARPX_GET_LOVECTS(warpx_getBfieldLoVects, WarpX::GetInstance().getBfield); + WARPX_GET_LOVECTS(warpx_getBfieldCPLoVects, WarpX::GetInstance().getBfield_cp); + WARPX_GET_LOVECTS(warpx_getBfieldFPLoVects, WarpX::GetInstance().getBfield_fp); + + WARPX_GET_LOVECTS(warpx_getCurrentDensityLoVects, WarpX::GetInstance().getcurrent); + WARPX_GET_LOVECTS(warpx_getCurrentDensityCPLoVects, WarpX::GetInstance().getcurrent_cp); + WARPX_GET_LOVECTS(warpx_getCurrentDensityFPLoVects, WarpX::GetInstance().getcurrent_fp); + +#define WARPX_GET_FIELD_PML(FIELD, GETTER) \ + amrex::Real** FIELD(int lev, int direction, \ + int *return_size, int *ncomps, int *ngrow, int **shapes) { \ + auto * pml = WarpX::GetInstance().GetPML(lev); \ + if (pml) { \ + auto & mf = *(pml->GETTER()[direction]); \ + return getMultiFabPointers(mf, return_size, ncomps, ngrow, shapes); \ + } else { \ + return nullptr; \ + } \ + } + +#define WARPX_GET_LOVECTS_PML(FIELD, GETTER) \ + int* FIELD(int lev, int direction, \ + int *return_size, int *ngrow) { \ + auto * pml = WarpX::GetInstance().GetPML(lev); \ + if (pml) { \ + auto & mf = *(pml->GETTER()[direction]); \ + return getMultiFabLoVects(mf, return_size, ngrow); \ + } else { \ + return nullptr; \ + } \ + } + + WARPX_GET_FIELD_PML(warpx_getEfieldCP_PML, GetE_cp); + WARPX_GET_FIELD_PML(warpx_getEfieldFP_PML, GetE_fp); + WARPX_GET_FIELD_PML(warpx_getBfieldCP_PML, GetB_cp); + WARPX_GET_FIELD_PML(warpx_getBfieldFP_PML, GetB_fp); + WARPX_GET_FIELD_PML(warpx_getCurrentDensityCP_PML, Getj_cp); + WARPX_GET_FIELD_PML(warpx_getCurrentDensityFP_PML, Getj_fp); + WARPX_GET_LOVECTS_PML(warpx_getEfieldCPLoVects_PML, GetE_cp); + WARPX_GET_LOVECTS_PML(warpx_getEfieldFPLoVects_PML, GetE_fp); + WARPX_GET_LOVECTS_PML(warpx_getBfieldCPLoVects_PML, GetB_cp); + WARPX_GET_LOVECTS_PML(warpx_getBfieldFPLoVects_PML, GetB_fp); + WARPX_GET_LOVECTS_PML(warpx_getCurrentDensityCPLoVects_PML, Getj_cp); + WARPX_GET_LOVECTS_PML(warpx_getCurrentDensityFPLoVects_PML, Getj_fp); amrex::ParticleReal** warpx_getParticleStructs(int speciesnumber, int lev, int* num_tiles, int** particles_per_tile) { diff --git a/Source/Python/WarpXWrappers.h b/Source/Python/WarpXWrappers.h index a272b9250..4de885b88 100644 --- a/Source/Python/WarpXWrappers.h +++ b/Source/Python/WarpXWrappers.h @@ -67,24 +67,6 @@ extern "C" { long warpx_getNumParticles(int speciesnumber); - amrex::Real** warpx_getEfield(int lev, int direction, - int *return_size, int* ncomps, int* ngrow, int **shapes); - - int* warpx_getEfieldLoVects(int lev, int direction, - int *return_size, int* ngrow); - - amrex::Real** warpx_getBfield(int lev, int direction, - int *return_size, int* ncomps, int* ngrow, int **shapes); - - int* warpx_getBfieldLoVects(int lev, int direction, - int *return_size, int* ngrow); - - amrex::Real** warpx_getCurrentDensity(int lev, int direction, - int *return_size, int* ncomps, int* ngrow, int **shapes); - - int* warpx_getCurrentDensityLoVects(int lev, int direction, - int *return_size, int* ngrow); - amrex::ParticleReal** warpx_getParticleStructs(int speciesnumber, int lev, int* num_tiles, int** particles_per_tile); diff --git a/Source/WarpX.H b/Source/WarpX.H index 7fff82445..decde3dc2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -207,6 +207,8 @@ public: void CopyJPML (); + PML* GetPML (int lev); + void PushParticlesandDepose (int lev, amrex::Real cur_time, DtType a_dt_type=DtType::Full); void PushParticlesandDepose ( amrex::Real cur_time); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 53728b54f..377d103d1 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1251,6 +1251,17 @@ WarpX::ComputeDivE (amrex::MultiFab& divE, int dcomp, } } +PML* +WarpX::GetPML (int lev) +{ + if (do_pml) { + // This should check if pml was initialized + return pml[lev].get(); + } else { + return nullptr; + } +} + void WarpX::BuildBufferMasks () { |