diff options
74 files changed, 1870 insertions, 584 deletions
diff --git a/Docs/source/building/building.rst b/Docs/source/building/building.rst index a9ea4e755..e0bdabb84 100644 --- a/Docs/source/building/building.rst +++ b/Docs/source/building/building.rst @@ -73,8 +73,8 @@ Advanced building instructions python spack -Building for specific plateforms --------------------------------- +Building for specific platforms +------------------------------- .. toctree:: :maxdepth: 1 diff --git a/Docs/source/building/cori.rst b/Docs/source/building/cori.rst index 43c5fc31c..3429ed576 100644 --- a/Docs/source/building/cori.rst +++ b/Docs/source/building/cori.rst @@ -51,6 +51,40 @@ In order to compile for the **Knight's Landing (KNL) architecture**: module swap PrgEnv-intel PrgEnv-gnu make -j 16 COMP=gnu +GPU Build +--------- + +To compile on the experimental GPU nodes on Cori, you first need to purge +your modules, most of which won't work on the GPU nodes. + + :: + + module purge + +Then, you need to load the following modules: + + :: + + module load esslurm cuda pgi openmpi/3.1.0-ucx + +Currently, you need to use OpenMPI; mvapich2 seems not to work. + +Then, you need to use slurm to request access to a GPU node: + + :: + + salloc -C gpu -N 1 -t 30 -c 10 --gres=gpu:1 --mem=30GB -A m1759 + +This reserves 10 logical cores (5 physical), 1 GPU, and 30 GB of RAM for your job. +Note that you can't cross-compile for the GPU nodes - you have to log on to one +and then build your software. + +Finally, navigate to the base of the WarpX repository and compile in GPU mode: + + :: + + make -j 16 COMP=pgi USE_GPU=TRUE + Building WarpX with openPMD support ----------------------------------- diff --git a/Docs/source/building/summit.rst b/Docs/source/building/summit.rst index 44cc5c0ee..c1061f29e 100644 --- a/Docs/source/building/summit.rst +++ b/Docs/source/building/summit.rst @@ -15,7 +15,7 @@ correct branch: git clone --branch master https://bitbucket.org/berkeleylab/picsar.git git clone --branch development https://github.com/AMReX-Codes/amrex.git -Then, use the following set of commands to compile: +Then, ``cd`` into the directory ``WarpX`` and use the following set of commands to compile: :: diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index 664c1c2ac..ba64c1fcd 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -29,12 +29,12 @@ Overall simulation parameters (The direction ``y`` cannot be used in 2D simulations.) * ``warpx.zmax_plasma_to_compute_max_step`` (`float`) optional - Can be useful when running in a boosted frame. If specified, automatically - calculates the number of iterations required in the boosted frame for the - lower `z` end of the simulation domain to reach - ``warpx.zmax_plasma_to_compute_max_step`` (typically the plasma end, - given in the lab frame). The value of ``max_step`` is overwritten, and - printed to standard output. Currently only works if the Lorentz boost and + Can be useful when running in a boosted frame. If specified, automatically + calculates the number of iterations required in the boosted frame for the + lower `z` end of the simulation domain to reach + ``warpx.zmax_plasma_to_compute_max_step`` (typically the plasma end, + given in the lab frame). The value of ``max_step`` is overwritten, and + printed to standard output. Currently only works if the Lorentz boost and the moving window are along the z direction. * ``warpx.verbose`` (`0` or `1`) @@ -188,8 +188,8 @@ Particle initialization This requires the additional parameter ``<species_name>.num_particles_per_cell``. * ``<species_name>.do_continuous_injection`` (`0` or `1`) - Whether to inject particles during the simulation, and not only at - initialization. This can be required whith a moving window and/or when + Whether to inject particles during the simulation, and not only at + initialization. This can be required whith a moving window and/or when running in a boosted frame. * ``<species_name>.profile`` (`string`) @@ -295,20 +295,20 @@ Particle initialization * ``<species>.plot_species`` (`0` or `1` optional; default `1`) Whether to plot particle quantities for this species. -* ``<species>.plot_vars`` (list of `strings` separated by spaces, optional) - List of particle quantities to write to `plotfiles`. By defaults, all - quantities are written to file. Choices are +* ``<species>.plot_vars`` (list of `strings` separated by spaces, optional) + List of particle quantities to write to `plotfiles`. By defaults, all + quantities are written to file. Choices are * ``w`` for the particle weight, - * ``ux`` ``uy`` ``uz`` for the particle momentum, + * ``ux`` ``uy`` ``uz`` for the particle momentum, * ``Ex`` ``Ey`` ``Ez`` for the electric field on particles, * ``Bx`` ``By`` ``Bz`` for the magnetic field on particles. - The particle positions are always included. Use - ``<species>.plot_vars = none`` to plot no particle data, except + The particle positions are always included. Use + ``<species>.plot_vars = none`` to plot no particle data, except particle position. * ``<species>.do_boosted_frame_diags`` (`0` or `1` optional, default `1`) - Only used when ``warpx.do_boosted_frame_diagnostic=1``. When running in a - boosted frame, whether or not to plot back-transformed diagnostics for + Only used when ``warpx.do_boosted_frame_diagnostic=1``. When running in a + boosted frame, whether or not to plot back-transformed diagnostics for this species. * ``warpx.serialize_ics`` (`0 or 1`) @@ -467,13 +467,13 @@ Laser initialization * ``<laser_name>.do_continuous_injection`` (`0` or `1`) optional (default `0`). Whether or not to use continuous injection (`0` or not `0`). - If the antenna starts outside of the simulation domain but enters it - at some point (due to moving window or moving antenna in the boosted - frame), use this so that the laser antenna is injected when it reaches - the box boundary. If running in a boosted frame, this requires the - boost direction, moving window direction and laser propagation direction - to be along `z`. If not running in a boosted frame, this requires the - moving window and laser propagation directions to be the same (`x`, `y` + If the antenna starts outside of the simulation domain but enters it + at some point (due to moving window or moving antenna in the boosted + frame), use this so that the laser antenna is injected when it reaches + the box boundary. If running in a boosted frame, this requires the + boost direction, moving window direction and laser propagation direction + to be along `z`. If not running in a boosted frame, this requires the + moving window and laser propagation directions to be the same (`x`, `y` or `z`) * ``warpx.num_mirrors`` (`int`) optional (default `0`) @@ -512,48 +512,60 @@ Numerics and algorithms Number of passes along each direction for the bilinear filter. In 2D simulations, only the first two values are read. -* ``algo.current_deposition`` (`integer`) - The algorithm for current deposition: +* ``algo.current_deposition`` (`string`, optional) + The algorithm for current deposition. Available options are: - - ``0``: Esirkepov deposition, vectorized - - ``1``: Esirkepov deposition, non-optimized - - ``2``: Direct deposition, vectorized - - ``3``: Direct deposition, non-optimized + - ``esirkepov``: the charge-conserving Esirkepov algorithm + (see `Esirkepov, Comp. Phys. Comm. (2001) <https://www.sciencedirect.com/science/article/pii/S0010465500002289>`__) + - ``direct``: simpler current deposition algorithm, described in + the section :doc:`../theory/picsar_theory`. Note that this algorithm is not strictly charge-conserving. + - ``direct-vectorized`` (only available in 3D, and when running on CPU/KNL - as opposed to GPU): + mathematically equivalent to ``direct``, but uses an optimized algorithm + for vectorization on CPU/KNL (see `Vincenti, Comp. Phys. Comm. (2017) + <https://www.sciencedirect.com/science/article/pii/S0010465516302764>`__) - .. warning:: + If ``algo.current_deposition`` is not specified, the default is ``esirkepov``. - On GPU, use ``algo.current_deposition=0`` for Esirkepov - or ``3`` for direct deposition. +* ``algo.charge_deposition`` (`string`, optional) + The algorithm for the charge density deposition. Available options are: -* ``algo.charge_deposition`` (`integer`) - The algorithm for the charge density deposition: + - ``standard``: standard charge deposition algorithm, described in + the section :doc:`../theory/picsar_theory`. + - ``vectorized`` (only available in 3D, and when running on CPU/KNL - as opposed to GPU): + mathematically equivalent to ``standard``, but uses an optimized algorithm + for vectorization on CPU/KNL (see `Vincenti, Comp. Phys. Comm. (2017) + <https://www.sciencedirect.com/science/article/pii/S0010465516302764>`__) - - ``0``: Vectorized version - - ``1``: Non-optimized version + If ``algo.charge_deposition`` is not specified, ``vectorized`` is the default + whenever it is available ; ``standard`` is the default otherwise. -* ``algo.field_gathering`` (`integer`) - The algorithm for field gathering: +* ``algo.field_gathering`` (`string`, optional) + The algorithm for field gathering. Available options are: - - ``0``: Vectorized version - - ``1``: Non-optimized version + - ``standard``: gathers directly from the grid points (either staggered + or nodal gridpoints depending on ``warpx.do_nodal``). + - ``vectorized`` (not available when running on GPU): mathematically + equivalent to ``standard``, but uses optimized vector instructions for CPU/KNL. - .. warning:: + If ``algo.field_gathering`` is not specified, ``vectorized`` is the default + on CPU/KNL ; ``standard`` is the default on GPU. + +* ``algo.particle_pusher`` (`string`, optional) + The algorithm for the particle pusher. Available options are: - The vectorized version does not run on GPU. Use - ``algo.field_gather=1`` when running on GPU. + - ``boris``: Boris pusher. + - ``vay``: Vay pusher (see `Vay, Phys. Plasmas (2008) <https://aip.scitation.org/doi/10.1063/1.2837054>`__) -* ``algo.particle_pusher`` (`integer`) - The algorithm for the particle pusher: + If ``algo.particle_pusher`` is not specified, ``boris`` is the default. - - ``0``: Boris pusher - - ``1``: Vay pusher +* ``algo.maxwell_fdtd_solver`` (`string`, optional) + The algorithm for the FDTD Maxwell field solver. Available options are: -* ``algo.maxwell_fdtd_solver`` (`string`) - The algorithm for the FDTD Maxwell field solver: + - ``yee``: Yee FDTD solver. + - ``ckc``: (not available in ``RZ`` geometry) Cole-Karkkainen solver with Cowan + coefficients (see `Cowan, PRSTAB 16 (2013) <https://journals.aps.org/prab/abstract/10.1103/PhysRevSTAB.16.041303>`__) - - ``yee``: Yee FDTD solver - - ``ckc``: Cole-Karkkainen solver with Cowan - coefficients (see Cowan - PRST-AB 16, 041303 (2013)) + If ``algo.maxwell_fdtd_solver`` is not specified, ``yee`` is the default. * ``interpolation.nox``, ``interpolation.noy``, ``interpolation.noz`` (`integer`) The order of the shape factors for the macroparticles, for the 3 dimensions of space. @@ -576,17 +588,17 @@ Numerics and algorithms fields are defined at different points in space) * ``warpx.do_subcycling`` (`0` or `1`; default: 0) - Whether or not to use sub-cycling. Different refinement levels have a - different cell size, which results in different Courant–Friedrichs–Lewy - (CFL) limits for the time step. By default, when using mesh refinement, - the same time step is used for all levels. This time step is - taken as the CFL limit of the finest level. Hence, for coarser - levels, the timestep is only a fraction of the CFL limit for this - level, which may lead to numerical artifacts. With sub-cycling, each level - evolves with its own time step, set to its own CFL limit. In practice, it - means that when level 0 performs one iteration, level 1 performs two - iterations. Currently, this option is only supported when - ``amr.max_level = 1``. More information can be found at + Whether or not to use sub-cycling. Different refinement levels have a + different cell size, which results in different Courant–Friedrichs–Lewy + (CFL) limits for the time step. By default, when using mesh refinement, + the same time step is used for all levels. This time step is + taken as the CFL limit of the finest level. Hence, for coarser + levels, the timestep is only a fraction of the CFL limit for this + level, which may lead to numerical artifacts. With sub-cycling, each level + evolves with its own time step, set to its own CFL limit. In practice, it + means that when level 0 performs one iteration, level 1 performs two + iterations. Currently, this option is only supported when + ``amr.max_level = 1``. More information can be found at https://ieeexplore.ieee.org/document/8659392. * ``psatd.nox``, ``psatd.noy``, ``pstad.noz`` (`integer`) optional (default `16` for all) @@ -653,7 +665,7 @@ Diagnostics and output The directory in which to save the lab frame data when using the **back-transformed diagnostics**. If not specified, the default is is `lab_frame_data`. - + * ``warpx.num_snapshots_lab`` (`integer`) Only used when ``warpx.do_boosted_frame_diagnostic`` is ``1``. The number of lab-frame snapshots that will be written. @@ -667,9 +679,9 @@ Diagnostics and output Whether to use the **back-transformed diagnostics** for the fields. * ``warpx.boosted_frame_diag_fields`` (space-separated list of `string`) - Which fields to dumped in back-transformed diagnostics. Choices are - 'Ex', 'Ey', Ez', 'Bx', 'By', Bz', 'jx', 'jy', jz' and 'rho'. Example: - ``warpx.boosted_frame_diag_fields = Ex Ez By``. By default, all fields + Which fields to dumped in back-transformed diagnostics. Choices are + 'Ex', 'Ey', Ez', 'Bx', 'By', Bz', 'jx', 'jy', jz' and 'rho'. Example: + ``warpx.boosted_frame_diag_fields = Ex Ez By``. By default, all fields are dumped. * ``warpx.plot_raw_fields`` (`0` or `1`) optional (default `0`) @@ -706,6 +718,19 @@ Diagnostics and output * ``warpx.plot_B_field`` (`0` or `1` optional; default `1`) Whether to plot the magnetic field. +* ``slice.dom_lo`` and ``slice.dom_hi`` (`2 floats in 2D`, `3 floats in 3D`; in meters similar to the units of the simulation box.) + The extent of the slice are defined by the co-ordinates of the lower corner (``slice.dom_lo``) and upper corner (``slice.dom_hi``). The slice could be 1D, 2D, or 3D, aligned with the co-ordinate axes and the first axis of the coordinates is x. For example: if for a 3D simulation, an x-z slice is to be extracted at y = 0.0, then the y-value of slice.dom_lo and slice.dom_hi must be equal to 0.0 + +* ``slice.coarsening_ratio`` (`2 integers in 2D`, `3 integers in 3D`; default `1`) + The coarsening ratio input must be greater than 0. Default is 1 in all directions. + In the directions that is reduced, i.e., for an x-z slice in 3D, the reduced y-dimension has a default coarsening ratio equal to 1. + +* ``slice.plot_int`` (`integer`) + The number of PIC cycles inbetween two consecutive data dumps for the slice. Use a + negative number to disable slice generation and slice data dumping. + + + Checkpoints and restart ----------------------- WarpX supports checkpoints/restart via AMReX. diff --git a/Docs/source/theory/picsar_theory.rst b/Docs/source/theory/picsar_theory.rst index 7338d5c36..135d78dea 100644 --- a/Docs/source/theory/picsar_theory.rst +++ b/Docs/source/theory/picsar_theory.rst @@ -328,6 +328,8 @@ a collocated and a staggered formulation is application-dependent. Spectral solvers used to be very popular in the years 1970s to early 1990s, before being replaced by finite-difference methods with the advent of parallel supercomputers that favored local methods. However, it was shown recently that standard domain decomposition with Fast Fourier Transforms that are local to each subdomain could be used effectively with PIC spectral methods (Jean-Luc Vay, Haber, and Godfrey 2013), at the cost of truncation errors in the guard cells that could be neglected. A detailed analysis of the effectiveness of the method with exact evaluation of the magnitude of the effect of the truncation error is given in (Vincenti and Vay 2016) for stencils of arbitrary order (up-to the infinite “spectral” order). +.. _current_deposition: + Current deposition ------------------ diff --git a/Examples/Modules/RigidInjection/inputs b/Examples/Modules/RigidInjection/inputs index d7e3b77e8..3c0d52ed0 100644 --- a/Examples/Modules/RigidInjection/inputs +++ b/Examples/Modules/RigidInjection/inputs @@ -23,10 +23,6 @@ geometry.prob_hi = 2. 2. 4. warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 # interpolation interpolation.nox = 3 diff --git a/Examples/Modules/boosted_diags/inputs.2d b/Examples/Modules/boosted_diags/inputs.2d index 33eee5392..528eb6cd9 100644 --- a/Examples/Modules/boosted_diags/inputs.2d +++ b/Examples/Modules/boosted_diags/inputs.2d @@ -24,10 +24,7 @@ geometry.prob_hi = 150.e-6 150.e-6 0. warpx.verbose = 1 # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = direct # Numerics interpolation.nox = 3 diff --git a/Examples/Modules/boosted_diags/inputs.3d b/Examples/Modules/boosted_diags/inputs.3d index 33eee5392..528eb6cd9 100644 --- a/Examples/Modules/boosted_diags/inputs.3d +++ b/Examples/Modules/boosted_diags/inputs.3d @@ -24,10 +24,7 @@ geometry.prob_hi = 150.e-6 150.e-6 0. warpx.verbose = 1 # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = direct # Numerics interpolation.nox = 3 diff --git a/Examples/Modules/charged_beam/inputs b/Examples/Modules/charged_beam/inputs index 36b991c3e..18b645281 100644 --- a/Examples/Modules/charged_beam/inputs +++ b/Examples/Modules/charged_beam/inputs @@ -23,10 +23,7 @@ geometry.prob_hi = 20.e-6 20.e-6 20.e-6 warpx.verbose = 1 # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = direct # CFL warpx.cfl = 1.0 diff --git a/Examples/Modules/gaussian_beam/gaussian_beam_PICMI.py b/Examples/Modules/gaussian_beam/gaussian_beam_PICMI.py index b9bfc2864..a3a2330d7 100644 --- a/Examples/Modules/gaussian_beam/gaussian_beam_PICMI.py +++ b/Examples/Modules/gaussian_beam/gaussian_beam_PICMI.py @@ -17,7 +17,7 @@ number_sim_particles = 32768 total_charge = 8.010883097437485e-07 beam_rms_size = 0.25 -electron_beam_divergence = -0.04 +electron_beam_divergence = -0.04*picmi.c em_order = 3 diff --git a/Examples/Modules/gaussian_beam/inputs b/Examples/Modules/gaussian_beam/inputs index 0cab8d140..46cd785f2 100644 --- a/Examples/Modules/gaussian_beam/inputs +++ b/Examples/Modules/gaussian_beam/inputs @@ -23,10 +23,6 @@ geometry.prob_hi = 2. 2. 2. warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 # interpolation interpolation.nox = 3 diff --git a/Examples/Modules/laser_injection/inputs b/Examples/Modules/laser_injection/inputs index 4186f2682..91ea135d1 100644 --- a/Examples/Modules/laser_injection/inputs +++ b/Examples/Modules/laser_injection/inputs @@ -21,10 +21,6 @@ geometry.prob_hi = 20.e-6 20.e-6 12.e-6 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 diff --git a/Examples/Modules/laser_injection/inputs.2d.rt b/Examples/Modules/laser_injection/inputs.2d.rt index 1c40525f8..db726a98c 100644 --- a/Examples/Modules/laser_injection/inputs.2d.rt +++ b/Examples/Modules/laser_injection/inputs.2d.rt @@ -23,10 +23,7 @@ warpx.serialize_ics = 1 warpx.verbose = 1 # Algorithms -algo.current_deposition = 1 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = esirkepov # CFL warpx.cfl = 1.0 diff --git a/Examples/Modules/laser_injection/inputs.rt b/Examples/Modules/laser_injection/inputs.rt index b6f49d661..cec710403 100644 --- a/Examples/Modules/laser_injection/inputs.rt +++ b/Examples/Modules/laser_injection/inputs.rt @@ -23,10 +23,6 @@ warpx.serialize_ics = 1 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 diff --git a/Examples/Modules/nci_corrector/inputs2d b/Examples/Modules/nci_corrector/inputs2d index 2f022f4cf..8db312fc6 100644 --- a/Examples/Modules/nci_corrector/inputs2d +++ b/Examples/Modules/nci_corrector/inputs2d @@ -25,10 +25,7 @@ warpx.verbose = 1 warpx.use_filter = 1 # Algorithms -algo.current_deposition = 1 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = esirkepov interpolation.nox = 3 interpolation.noy = 3 interpolation.noz = 3 diff --git a/Examples/Physics_applications/laser_acceleration/inputs.2d b/Examples/Physics_applications/laser_acceleration/inputs.2d index 11ddc6b0b..1a76feacf 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs.2d +++ b/Examples/Physics_applications/laser_acceleration/inputs.2d @@ -17,10 +17,6 @@ warpx.fine_tag_hi = 5.e-6 -25.e-6 ################################# ############ NUMERICS ########### ################################# -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z interpolation.noy = 3 interpolation.noz = 3 diff --git a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost index 18c354634..fb03ae3f6 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost +++ b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost @@ -1,8 +1,8 @@ ################################# ######### BOX PARAMETERS ######## ################################# -# max_step = 2700 -stop_time = 1.9e-12 +max_step = 1000 +# stop_time = 1.9e-12 amr.n_cell = 128 1024 amr.max_grid_size = 64 amr.blocking_factor = 32 @@ -19,10 +19,8 @@ geometry.prob_hi = 128.e-6 0.96e-6 ################################# warpx.verbose = 1 amrex.v = 1 -algo.current_deposition = 2 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 1 +algo.current_deposition = direct-vectorized +algo.particle_pusher = vay algo.maxwell_fdtd_solver = ckc interpolation.nox = 3 interpolation.noy = 3 @@ -39,7 +37,7 @@ warpx.serialize_ics = 1 ################################# ####### BOOST PARAMETERS ######## ################################# -warpx.gamma_boost = 10. +warpx.gamma_boost = 30. warpx.boost_direction = z warpx.do_boosted_frame_diagnostic = 1 warpx.num_snapshots_lab = 7 @@ -62,11 +60,11 @@ electrons.momentum_distribution_type = "gaussian" electrons.xmin = -120.e-6 electrons.xmax = 120.e-6 electrons.zmin = 0.5e-3 -electrons.zmax = .0035 +electrons.zmax = 1. electrons.profile = "predefined" electrons.predefined_profile_name = "parabolic_channel" # predefined_profile_params = z_start ramp_up plateau ramp_down rc n0 -electrons.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 +electrons.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e25 electrons.do_continuous_injection = 1 ions.charge = q_e @@ -77,11 +75,11 @@ ions.momentum_distribution_type = "gaussian" ions.xmin = -120.e-6 ions.xmax = 120.e-6 ions.zmin = 0.5e-3 -ions.zmax = .0035 +ions.zmax = 1. ions.profile = "predefined" ions.predefined_profile_name = "parabolic_channel" # predefined_profile_params = z_start ramp_up plateau ramp_down rc n0 -ions.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 +ions.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e25 ions.do_continuous_injection = 1 beam.charge = -q_e diff --git a/Examples/Physics_applications/laser_acceleration/inputs.3d b/Examples/Physics_applications/laser_acceleration/inputs.3d index 8e854cbf4..f055708d1 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs.3d +++ b/Examples/Physics_applications/laser_acceleration/inputs.3d @@ -17,10 +17,6 @@ warpx.fine_tag_hi = 5.e-6 5.e-6 -30.e-6 ################################# ############ NUMERICS ########### ################################# -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z interpolation.noy = 3 interpolation.noz = 3 diff --git a/Examples/Physics_applications/plasma_acceleration/inputs.2d b/Examples/Physics_applications/plasma_acceleration/inputs.2d index e092cf932..9d0edbc7e 100644 --- a/Examples/Physics_applications/plasma_acceleration/inputs.2d +++ b/Examples/Physics_applications/plasma_acceleration/inputs.2d @@ -18,10 +18,6 @@ warpx.fine_tag_hi = 12.e-6 -100.e-6 ################################# ############ NUMERICS ########### ################################# -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 algo.maxwell_fdtd_solver = ckc warpx.use_filter = 1 warpx.do_pml = 1 diff --git a/Examples/Physics_applications/plasma_acceleration/inputs.2d.boost b/Examples/Physics_applications/plasma_acceleration/inputs.2d.boost index 59e7c432b..ba9166dee 100644 --- a/Examples/Physics_applications/plasma_acceleration/inputs.2d.boost +++ b/Examples/Physics_applications/plasma_acceleration/inputs.2d.boost @@ -16,10 +16,6 @@ geometry.prob_hi = 125.e-6 1.e-6 ################################# ############ NUMERICS ########### ################################# -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 algo.maxwell_fdtd_solver = ckc warpx.verbose = 1 warpx.plot_raw_fields = 0 diff --git a/Examples/Physics_applications/plasma_acceleration/inputs.3d.boost b/Examples/Physics_applications/plasma_acceleration/inputs.3d.boost index 1d58ef940..1af4fac86 100644 --- a/Examples/Physics_applications/plasma_acceleration/inputs.3d.boost +++ b/Examples/Physics_applications/plasma_acceleration/inputs.3d.boost @@ -16,10 +16,6 @@ geometry.prob_hi = 0.00015 0.00015 1.e-06 ################################# ############ NUMERICS ########### ################################# -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 algo.maxwell_fdtd_solver = ckc warpx.verbose = 1 warpx.plot_raw_fields = 0 diff --git a/Examples/Physics_applications/plasma_mirror/inputs.2d b/Examples/Physics_applications/plasma_mirror/inputs.2d index e16313b42..dc0a48ebf 100644 --- a/Examples/Physics_applications/plasma_mirror/inputs.2d +++ b/Examples/Physics_applications/plasma_mirror/inputs.2d @@ -18,10 +18,6 @@ warpx.serialize_ics = 1 ################################# ############ NUMERICS ########### ################################# -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 interpolation.nox = 3 interpolation.noy = 3 interpolation.noz = 3 diff --git a/Examples/Physics_applications/uniform_plasma/inputs.2d b/Examples/Physics_applications/uniform_plasma/inputs.2d index f3eb2decb..4304638dd 100644 --- a/Examples/Physics_applications/uniform_plasma/inputs.2d +++ b/Examples/Physics_applications/uniform_plasma/inputs.2d @@ -16,10 +16,6 @@ geometry.prob_hi = 20.e-6 20.e-6 ################################# warpx.serialize_ics = 1 warpx.verbose = 1 -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 warpx.cfl = 1.0 amr.plot_int = 10 diff --git a/Examples/Physics_applications/uniform_plasma/inputs.3d b/Examples/Physics_applications/uniform_plasma/inputs.3d index 245f3204c..6792dac2c 100644 --- a/Examples/Physics_applications/uniform_plasma/inputs.3d +++ b/Examples/Physics_applications/uniform_plasma/inputs.3d @@ -16,10 +16,6 @@ geometry.prob_hi = 20.e-6 20.e-6 20.e-6 ################################# warpx.serialize_ics = 1 warpx.verbose = 1 -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 warpx.cfl = 1.0 amr.plot_int = 10 diff --git a/Examples/Tests/Langmuir/inputs b/Examples/Tests/Langmuir/inputs index cd13aebbd..a989f2249 100644 --- a/Examples/Tests/Langmuir/inputs +++ b/Examples/Tests/Langmuir/inputs @@ -26,10 +26,7 @@ warpx.moving_window_dir = x warpx.moving_window_v = 0.0 # in units of the speed of light # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = direct # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/Langmuir/inputs.2d.rz b/Examples/Tests/Langmuir/inputs.2d.rz index c4132d493..a4464fe1d 100644 --- a/Examples/Tests/Langmuir/inputs.2d.rz +++ b/Examples/Tests/Langmuir/inputs.2d.rz @@ -27,10 +27,7 @@ warpx.moving_window_dir = z warpx.moving_window_v = 0.0 # in units of the speed of light # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = direct interpolation.nox = 1 interpolation.noy = 1 interpolation.noz = 1 diff --git a/Examples/Tests/Langmuir/inputs.lb b/Examples/Tests/Langmuir/inputs.lb index 5900cb1bb..73ededb23 100644 --- a/Examples/Tests/Langmuir/inputs.lb +++ b/Examples/Tests/Langmuir/inputs.lb @@ -32,10 +32,7 @@ warpx.moving_window_v = 0.0 # in units of the speed of light warpx.regrid_int = 1 # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = direct # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/Langmuir/inputs.multi.2d.rt b/Examples/Tests/Langmuir/inputs.multi.2d.rt index 79bc2b383..54d4e2c16 100644 --- a/Examples/Tests/Langmuir/inputs.multi.2d.rt +++ b/Examples/Tests/Langmuir/inputs.multi.2d.rt @@ -25,10 +25,7 @@ warpx.serialize_ics = 1 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.field_gathering = standard # Interpolation interpolation.nox = 1 diff --git a/Examples/Tests/Langmuir/inputs.multi.rt b/Examples/Tests/Langmuir/inputs.multi.rt index c9969d39a..8e6d6d8e1 100644 --- a/Examples/Tests/Langmuir/inputs.multi.rt +++ b/Examples/Tests/Langmuir/inputs.multi.rt @@ -26,10 +26,8 @@ warpx.serialize_ics = 1 warpx.verbose = 1 # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.current_deposition = direct +algo.field_gathering = standard # Interpolation interpolation.nox = 1 diff --git a/Examples/Tests/Langmuir/inputs.multi.rz.rt b/Examples/Tests/Langmuir/inputs.multi.rz.rt index 6bc84a434..e4099f9c5 100644 --- a/Examples/Tests/Langmuir/inputs.multi.rz.rt +++ b/Examples/Tests/Langmuir/inputs.multi.rz.rt @@ -25,10 +25,7 @@ warpx.serialize_ics = 1 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.field_gathering = standard # Interpolation interpolation.nox = 1 diff --git a/Examples/Tests/Langmuir/inputs.nolb b/Examples/Tests/Langmuir/inputs.nolb index cd2900b29..6ac0629cb 100644 --- a/Examples/Tests/Langmuir/inputs.nolb +++ b/Examples/Tests/Langmuir/inputs.nolb @@ -32,10 +32,7 @@ warpx.moving_window_v = 0.0 # in units of the speed of light warpx.regrid_int = -1 # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 +algo.current_deposition = direct # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/Langmuir/inputs.rt b/Examples/Tests/Langmuir/inputs.rt index 3ec8fc47d..be5a730c2 100644 --- a/Examples/Tests/Langmuir/inputs.rt +++ b/Examples/Tests/Langmuir/inputs.rt @@ -25,10 +25,7 @@ warpx.serialize_ics = 1 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.field_gathering = standard # Interpolation interpolation.nox = 1 diff --git a/Examples/Tests/Langmuir/langmuir2d_PICMI.py b/Examples/Tests/Langmuir/langmuir2d_PICMI.py index c8c5a9559..8108a39ba 100644 --- a/Examples/Tests/Langmuir/langmuir2d_PICMI.py +++ b/Examples/Tests/Langmuir/langmuir2d_PICMI.py @@ -9,7 +9,9 @@ ymin = -20.e-6 xmax = +20.e-6 ymax = +20.e-6 -uniform_plasma = picmi.UniformDistribution(density=1.e25, upper_bound=[0., None, None], directed_velocity=[0.1, 0., 0.]) +uniform_plasma = picmi.UniformDistribution(density = 1.e25, + upper_bound = [0., None, None], + directed_velocity = [0.1*picmi.c, 0., 0.]) electrons = picmi.Species(particle_type='electron', name='electrons', initial_distribution=uniform_plasma) diff --git a/Examples/Tests/Langmuir/langmuir_PICMI.py b/Examples/Tests/Langmuir/langmuir_PICMI.py index 67f37ec65..d09c90714 100644 --- a/Examples/Tests/Langmuir/langmuir_PICMI.py +++ b/Examples/Tests/Langmuir/langmuir_PICMI.py @@ -14,7 +14,9 @@ xmax = +20.e-6 ymax = +20.e-6 zmax = +20.e-6 -uniform_plasma = picmi.UniformDistribution(density=1.e25, upper_bound=[0., None, None], directed_velocity=[0.1, 0., 0.]) +uniform_plasma = picmi.UniformDistribution(density = 1.e25, + upper_bound = [0., None, None], + directed_velocity = [0.1*picmi.c, 0., 0.]) electrons = picmi.Species(particle_type='electron', name='electrons', initial_distribution=uniform_plasma) diff --git a/Examples/Tests/Langmuir/langmuir_PICMI_rt.py b/Examples/Tests/Langmuir/langmuir_PICMI_rt.py index 30540313e..e629e268d 100644 --- a/Examples/Tests/Langmuir/langmuir_PICMI_rt.py +++ b/Examples/Tests/Langmuir/langmuir_PICMI_rt.py @@ -14,7 +14,9 @@ xmax = +20.e-6 ymax = +20.e-6 zmax = +20.e-6 -uniform_plasma = picmi.UniformDistribution(density=1.e25, upper_bound=[0., None, None], directed_velocity=[0.1, 0., 0.]) +uniform_plasma = picmi.UniformDistribution(density = 1.e25, + upper_bound = [0., None, None], + directed_velocity = [0.1*picmi.c, 0., 0.]) electrons = picmi.Species(particle_type='electron', name='electrons', initial_distribution=uniform_plasma) diff --git a/Examples/Tests/Langmuir/langmuir_PICMI_rz.py b/Examples/Tests/Langmuir/langmuir_PICMI_rz.py index 4eb998177..88b3947a6 100644 --- a/Examples/Tests/Langmuir/langmuir_PICMI_rz.py +++ b/Examples/Tests/Langmuir/langmuir_PICMI_rz.py @@ -9,7 +9,9 @@ zmin = -20.e-6 rmax = +20.e-6 zmax = +20.e-6 -uniform_plasma = picmi.UniformDistribution(density=1.e25, upper_bound=[None, None, 0.], directed_velocity=[0., 0., 0.1]) +uniform_plasma = picmi.UniformDistribution(density = 1.e25, + upper_bound = [None, None, 0.], + directed_velocity = [0., 0., 0.1*picmi.c]) electrons = picmi.Species(particle_type='electron', name='electrons', initial_distribution=uniform_plasma) diff --git a/Examples/Tests/Larmor/inputs b/Examples/Tests/Larmor/inputs index 729d315e9..da8fb5bae 100644 --- a/Examples/Tests/Larmor/inputs +++ b/Examples/Tests/Larmor/inputs @@ -41,10 +41,6 @@ warpx.B_external = 0.0 0.00078110417851950768 0.0 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/Larmor/inputs.ml b/Examples/Tests/Larmor/inputs.ml index 9a0c9eb3d..3cfb978d7 100644 --- a/Examples/Tests/Larmor/inputs.ml +++ b/Examples/Tests/Larmor/inputs.ml @@ -41,10 +41,6 @@ warpx.B_external = 0.0 0.00078110417851950768 0.0 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/PML/inputs2d b/Examples/Tests/PML/inputs2d index c84cf0555..5b936a333 100644 --- a/Examples/Tests/PML/inputs2d +++ b/Examples/Tests/PML/inputs2d @@ -23,10 +23,6 @@ geometry.prob_hi = 30.e-6 120.e-6 warpx.verbose = 1 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 warpx.cfl = 1.0 warpx.do_pml = 1 diff --git a/Examples/Tests/SingleParticle/bilinear_filter_analysis.py b/Examples/Tests/SingleParticle/bilinear_filter_analysis.py new file mode 100755 index 000000000..1a47796fb --- /dev/null +++ b/Examples/Tests/SingleParticle/bilinear_filter_analysis.py @@ -0,0 +1,38 @@ +#! /usr/bin/env python + +import sys +import yt ; yt.funcs.mylog.setLevel(0) +import numpy as np +from scipy import signal + +# Build Jx without filter (from other simulation) +my_F_nofilter = np.zeros([16,16]) +my_F_nofilter[8,8] = -1.601068065642412e-11 +my_F_nofilter[8,7] = -1.601068065642412e-11 + +# Build 2D filter +filter0 = np.array([.25,.5,.25]) +my_order = [1,5] +my_filterx = filter0 +my_filtery = filter0 +while my_order[0]>1: + my_filterx = np.convolve(my_filterx,filter0) + my_order[0] -= 1 +while my_order[1]>1: + my_filtery = np.convolve(my_filtery,filter0) + my_order[1] -= 1 +my_filter = my_filterx[:,None]*my_filtery + +# Apply filter. my_F_filetered is the theoretical value for filtered field +my_F_filtered = signal.convolve2d(my_F_nofilter, my_filter, boundary='symm', mode='same') + +# Get simulation result for F_filtered +filename = sys.argv[1] +ds = yt.load( filename ) +sl = yt.SlicePlot(ds, 2, 'jx', aspect=1) +all_data_level_0 = ds.covering_grid(level=0,left_edge=ds.domain_left_edge, dims=ds.domain_dimensions) +F_filtered = all_data_level_0['boxlib', 'jx'].v.squeeze() + +# Compare theory and PIC for filtered value +error = np.sum( np.abs(F_filtered - my_F_filtered) ) / np.sum( np.abs(my_F_filtered) ) +assert( error < 1.e-14 ) diff --git a/Examples/Tests/SingleParticle/inputs b/Examples/Tests/SingleParticle/inputs new file mode 100644 index 000000000..548848d79 --- /dev/null +++ b/Examples/Tests/SingleParticle/inputs @@ -0,0 +1,23 @@ +max_step = 1 +amr.n_cell = 16 16 +amr.max_level = 0 +amr.blocking_factor = 8 +amr.max_grid_size = 8 +amr.plot_int = 1 +geometry.coord_sys = 0 +geometry.is_periodic = 0 0 +geometry.prob_lo = -8 -12 +geometry.prob_hi = 8 12 +warpx.do_pml = 0 +algo.charge_deposition = standard +algo.field_gathering = standard +warpx.cfl = 1.0 + +particles.nspecies = 1 +particles.species_names = electron +electron.charge = -q_e +electron.mass = m_e +electron.injection_style = "SingleParticle" +electron.single_particle_pos = 0.0 0.0 0.0 +electron.single_particle_vel = 1.e20 0.0 0.0 # gamma*beta +electron.single_particle_weight = 1.0 diff --git a/Examples/Tests/gpu_test/inputs b/Examples/Tests/gpu_test/inputs index 2b719655d..e4ae27469 100644 --- a/Examples/Tests/gpu_test/inputs +++ b/Examples/Tests/gpu_test/inputs @@ -24,10 +24,8 @@ warpx.do_pml = 0 warpx.verbose = 1 # Algorithms -algo.current_deposition = 3 -algo.charge_deposition = 0 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.current_deposition = direct +algo.field_gathering = standard interpolation.nox = 1 interpolation.noy = 1 diff --git a/Examples/Tests/laser_on_fine/inputs b/Examples/Tests/laser_on_fine/inputs index 8fc8876cb..9d3133dd7 100644 --- a/Examples/Tests/laser_on_fine/inputs +++ b/Examples/Tests/laser_on_fine/inputs @@ -38,10 +38,9 @@ warpx.pml_ncell = 10 warpx.verbose = 1 # Algorithms -algo.current_deposition = 1 -algo.charge_deposition = 1 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = standard # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/laser_on_fine/inputs.2d b/Examples/Tests/laser_on_fine/inputs.2d index c4af2df56..d011ff4ee 100644 --- a/Examples/Tests/laser_on_fine/inputs.2d +++ b/Examples/Tests/laser_on_fine/inputs.2d @@ -38,10 +38,9 @@ warpx.pml_ncell = 10 warpx.verbose = 1 # Algorithms -algo.current_deposition = 1 -algo.charge_deposition = 1 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = standard # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/self_force_test/inputs b/Examples/Tests/self_force_test/inputs index 1e55c1174..5e7bcc4f8 100644 --- a/Examples/Tests/self_force_test/inputs +++ b/Examples/Tests/self_force_test/inputs @@ -42,10 +42,9 @@ warpx.pml_ncell = 10 warpx.verbose = 1 # Algorithms -algo.current_deposition = 1 -algo.charge_deposition = 1 -algo.field_gathering = 1 -algo.particle_pusher = 0 +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = standard # CFL warpx.cfl = 1.0 diff --git a/Examples/Tests/subcycling/inputs.2d b/Examples/Tests/subcycling/inputs.2d index 9790ee973..15034be52 100644 --- a/Examples/Tests/subcycling/inputs.2d +++ b/Examples/Tests/subcycling/inputs.2d @@ -37,10 +37,6 @@ warpx.n_current_deposition_buffer = 0 warpx.n_field_gather_buffer = 0 # Algorithms -algo.current_deposition = 0 -algo.charge_deposition = 0 -algo.field_gathering = 0 -algo.particle_pusher = 0 algo.maxwell_fdtd_solver = "ckc" # CFL diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py index f242f8589..f6dc47fec 100644 --- a/Python/pywarpx/picmi.py +++ b/Python/pywarpx/picmi.py @@ -94,24 +94,25 @@ class GaussianBunchDistribution(picmistandard.PICMI_GaussianBunchDistribution): #species.zmin #species.zmax + # --- Note that WarpX takes gamma*beta as input if np.any(np.not_equal(self.velocity_divergence, 0.)): species.momentum_distribution_type = "radial_expansion" - species.u_over_r = self.velocity_divergence[0] - #species.u_over_y = self.velocity_divergence[1] - #species.u_over_z = self.velocity_divergence[2] + species.u_over_r = self.velocity_divergence[0]/c + #species.u_over_y = self.velocity_divergence[1]/c + #species.u_over_z = self.velocity_divergence[2]/c elif np.any(np.not_equal(self.rms_velocity, 0.)): species.momentum_distribution_type = "gaussian" - species.ux_m = self.centroid_velocity[0] - species.uy_m = self.centroid_velocity[1] - species.uz_m = self.centroid_velocity[2] - species.ux_th = self.rms_velocity[0] - species.uy_th = self.rms_velocity[1] - species.uz_th = self.rms_velocity[2] + species.ux_m = self.centroid_velocity[0]/c + species.uy_m = self.centroid_velocity[1]/c + species.uz_m = self.centroid_velocity[2]/c + species.ux_th = self.rms_velocity[0]/c + species.uy_th = self.rms_velocity[1]/c + species.uz_th = self.rms_velocity[2]/c else: species.momentum_distribution_type = "constant" - species.ux = self.centroid_velocity[0] - species.uy = self.centroid_velocity[1] - species.uz = self.centroid_velocity[2] + species.ux = self.centroid_velocity[0]/c + species.uy = self.centroid_velocity[1]/c + species.uz = self.centroid_velocity[2]/c class UniformDistribution(picmistandard.PICMI_UniformDistribution): @@ -139,19 +140,20 @@ class UniformDistribution(picmistandard.PICMI_UniformDistribution): species.profile = "constant" species.density = self.density + # --- Note that WarpX takes gamma*beta as input if np.any(np.not_equal(self.rms_velocity, 0.)): species.momentum_distribution_type = "gaussian" - species.ux_m = self.directed_velocity[0] - species.uy_m = self.directed_velocity[1] - species.uz_m = self.directed_velocity[2] - species.ux_th = self.rms_velocity[0] - species.uy_th = self.rms_velocity[1] - species.uz_th = self.rms_velocity[2] + species.ux_m = self.directed_velocity[0]/c + species.uy_m = self.directed_velocity[1]/c + species.uz_m = self.directed_velocity[2]/c + species.ux_th = self.rms_velocity[0]/c + species.uy_th = self.rms_velocity[1]/c + species.uz_th = self.rms_velocity[2]/c else: species.momentum_distribution_type = "constant" - species.ux = self.directed_velocity[0] - species.uy = self.directed_velocity[1] - species.uz = self.directed_velocity[2] + species.ux = self.directed_velocity[0]/c + species.uy = self.directed_velocity[1]/c + species.uz = self.directed_velocity[2]/c if self.fill_in: species.do_continuous_injection = 1 @@ -185,19 +187,20 @@ class AnalyticDistribution(picmistandard.PICMI_AnalyticDistribution): for k,v in self.user_defined_kw.items(): setattr(pywarpx.my_constants, k, v) + # --- Note that WarpX takes gamma*beta as input if np.any(np.not_equal(self.rms_velocity, 0.)): species.momentum_distribution_type = "gaussian" - species.ux_m = self.directed_velocity[0] - species.uy_m = self.directed_velocity[1] - species.uz_m = self.directed_velocity[2] - species.ux_th = self.rms_velocity[0] - species.uy_th = self.rms_velocity[1] - species.uz_th = self.rms_velocity[2] + species.ux_m = self.directed_velocity[0]/c + species.uy_m = self.directed_velocity[1]/c + species.uz_m = self.directed_velocity[2]/c + species.ux_th = self.rms_velocity[0]/c + species.uy_th = self.rms_velocity[1]/c + species.uz_th = self.rms_velocity[2]/c else: species.momentum_distribution_type = "constant" - species.ux = self.directed_velocity[0] - species.uy = self.directed_velocity[1] - species.uz = self.directed_velocity[2] + species.ux = self.directed_velocity[0]/c + species.uy = self.directed_velocity[1]/c + species.uz = self.directed_velocity[2]/c if self.fill_in: species.do_continuous_injection = 1 diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index 34fbc2769..103bec201 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -102,6 +102,21 @@ compileTest = 0 doVis = 0 analysisRoutine = Examples/Modules/nci_corrector/ncicorr_analysis.py +[bilinear_filter] +buildDir = . +inputFile = Examples/Tests/SingleParticle/inputs +runtime_params = warpx.use_filter=1 warpx.filter_npass_each_dir=1 5 +dim = 2 +addToCompileString = +restartTest = 0 +useMPI = 1 +numprocs = 2 +useOMP = 1 +numthreads = 2 +compileTest = 0 +doVis = 0 +analysisRoutine = Examples/Tests/SingleParticle/bilinear_filter_analysis.py + [Langmuir_2d] buildDir = . inputFile = Examples/Tests/Langmuir/inputs.rt diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index e3d44d1fc..b1181f22f 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -677,7 +677,7 @@ getInterpolatedScalar( if ( F_fp.is_nodal() ){ IntVect refinement_vector{AMREX_D_DECL(r_ratio, r_ratio, r_ratio)}; node_bilinear_interp.interp(cfab, 0, ffab, 0, 1, - finebx, refinement_vector, {}, {}, {}, 0, 0); + finebx, refinement_vector, {}, {}, {}, 0, 0, RunOn::Cpu); } else { amrex::Abort("Unknown field staggering."); } diff --git a/Source/Diagnostics/Make.package b/Source/Diagnostics/Make.package index 3b3b3d202..dfd947d53 100644 --- a/Source/Diagnostics/Make.package +++ b/Source/Diagnostics/Make.package @@ -5,6 +5,8 @@ CEXE_sources += FieldIO.cpp CEXE_headers += FieldIO.H CEXE_headers += BoostedFrameDiagnostic.H CEXE_headers += ElectrostaticIO.cpp +CEXE_headers += SliceDiagnostic.H +CEXE_sources += SliceDiagnostic.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics VPATH_LOCATIONS += $(WARPX_HOME)/Source/Diagnostics diff --git a/Source/Diagnostics/SliceDiagnostic.H b/Source/Diagnostics/SliceDiagnostic.H new file mode 100644 index 000000000..31eea83be --- /dev/null +++ b/Source/Diagnostics/SliceDiagnostic.H @@ -0,0 +1,42 @@ +#ifndef WARPX_SliceDiagnostic_H_ +#define WARPX_SliceDiagnostic_H_ + +#include <AMReX_VisMF.H> +#include <AMReX_PlotFileUtil.H> +#include <AMReX_ParallelDescriptor.H> +#include <AMReX_Geometry.H> + +#include <WarpX.H> +#include <AMReX_FArrayBox.H> +#include <AMReX_IArrayBox.H> +#include <AMReX_Vector.H> +#include <AMReX_BLassert.H> +#include <AMReX_MultiFabUtil.H> +#include <AMReX_MultiFabUtil_C.H> + + + +std::unique_ptr<amrex::MultiFab> CreateSlice( const amrex::MultiFab& mf, + const amrex::Vector<amrex::Geometry> &dom_geom, + amrex::RealBox &slice_realbox, + amrex::IntVect &slice_cr_ratio ); + +void CheckSliceInput( const amrex::RealBox real_box, + amrex::RealBox &slice_cc_nd_box, amrex::RealBox &slice_realbox, + amrex::IntVect &slice_cr_ratio, amrex::Vector<amrex::Geometry> dom_geom, + amrex::IntVect const SliceType, amrex::IntVect &slice_lo, + amrex::IntVect &slice_hi, amrex::IntVect &interp_lo); + +void InterpolateSliceValues( amrex::MultiFab& smf, + amrex::IntVect interp_lo, amrex::RealBox slice_realbox, + amrex::Vector<amrex::Geometry> geom, int ncomp, int nghost, + amrex::IntVect slice_lo, amrex::IntVect slice_hi, + amrex::IntVect SliceType, const amrex::RealBox real_box); + +void InterpolateLo( const amrex::Box& bx, amrex::FArrayBox &fabox, + amrex::IntVect slice_lo, amrex::Vector<amrex::Geometry> geom, + int idir, amrex::IntVect IndType, amrex::RealBox slice_realbox, + int srccomp, int ncomp, int nghost, const amrex::RealBox real_box); + +#endif + diff --git a/Source/Diagnostics/SliceDiagnostic.cpp b/Source/Diagnostics/SliceDiagnostic.cpp new file mode 100644 index 000000000..994f990c6 --- /dev/null +++ b/Source/Diagnostics/SliceDiagnostic.cpp @@ -0,0 +1,475 @@ +#include "SliceDiagnostic.H" +#include <AMReX_MultiFabUtil.H> +#include <AMReX_PlotFileUtil.H> +#include <AMReX_FillPatchUtil_F.H> + +#include <WarpX.H> + +using namespace amrex; + + +/* \brief + * The functions creates the slice for diagnostics based on the user-input. + * The slice can be 1D, 2D, or 3D and it inherts the index type of the underlying data. + * The implementation assumes that the slice is aligned with the coordinate axes. + * The input parameters are modified if the user-input does not comply with requirements of coarsenability or if the slice extent is not contained within the simulation domain. + * First a slice multifab (smf) with cell size equal to that of the simulation grid is created such that it extends from slice.dim_lo to slice.dim_hi and shares the same index space as the source multifab (mf) + * The values are copied from src mf to dst smf using amrex::ParallelCopy + * If interpolation is required, then on the smf, using data points stored in the ghost cells, the data in interpolated. + * If coarsening is required, then a coarse slice multifab is generated (cs_mf) and the + * values of the refined slice (smf) is averaged down to obtain the coarse slice. + * \param mf is the source multifab containing the field data + * \param dom_geom is the geometry of the domain and used in the function to obtain the + * CellSize of the underlying grid. + * \param slice_realbox defines the extent of the slice + * \param slice_cr_ratio provides the coarsening ratio for diagnostics + */ + +std::unique_ptr<MultiFab> +CreateSlice( const MultiFab& mf, const Vector<Geometry> &dom_geom, + RealBox &slice_realbox, IntVect &slice_cr_ratio ) +{ + std::unique_ptr<MultiFab> smf; + std::unique_ptr<MultiFab> cs_mf; + + Vector<int> slice_ncells(AMREX_SPACEDIM); + int nghost = 1; + int nlevels = dom_geom.size(); + int ncomp = (mf).nComp(); + + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nlevels==1, + "Slice diagnostics does not work with mesh refinement yet (TO DO)."); + + const auto conversionType = (mf).ixType(); + IntVect SliceType(AMREX_D_DECL(0,0,0)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim ) + { + SliceType[idim] = conversionType.nodeCentered(idim); + } + + const RealBox& real_box = dom_geom[0].ProbDomain(); + RealBox slice_cc_nd_box; + int slice_grid_size = 32; + + bool interpolate = false; + bool coarsen = false; + + // same index space as domain // + IntVect slice_lo(AMREX_D_DECL(0,0,0)); + IntVect slice_hi(AMREX_D_DECL(1,1,1)); + IntVect interp_lo(AMREX_D_DECL(0,0,0)); + + CheckSliceInput(real_box, slice_cc_nd_box, slice_realbox, slice_cr_ratio, + dom_geom, SliceType, slice_lo, + slice_hi, interp_lo); + // Determine if interpolation is required and number of cells in slice // + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + + // Flag for interpolation if required // + if ( interp_lo[idim] == 1) { + interpolate = 1; + } + + // For the case when a dimension is reduced // + if ( ( slice_hi[idim] - slice_lo[idim]) == 1) { + slice_ncells[idim] = 1; + } + else { + slice_ncells[idim] = ( slice_hi[idim] - slice_lo[idim] + 1 ) + / slice_cr_ratio[idim]; + + int refined_ncells = slice_hi[idim] - slice_lo[idim] + 1 ; + if ( slice_cr_ratio[idim] > 1) { + coarsen = true; + + // modify slice_grid_size if >= refines_cells // + if ( slice_grid_size >= refined_ncells ) { + slice_grid_size = refined_ncells - 1; + } + + } + } + } + + // Slice generation with index type inheritance // + Box slice(slice_lo, slice_hi); + + Vector<BoxArray> sba(1); + sba[0].define(slice); + sba[0].maxSize(slice_grid_size); + + // Distribution mapping for slice can be different from that of domain // + Vector<DistributionMapping> sdmap(1); + sdmap[0] = DistributionMapping{sba[0]}; + + smf.reset(new MultiFab(amrex::convert(sba[0],SliceType), sdmap[0], + ncomp, nghost)); + + // Copy data from domain to slice that has same cell size as that of // + // the domain mf. src and dst have the same number of ghost cells // + smf->ParallelCopy(mf, 0, 0, ncomp,nghost,nghost); + + // inteprolate if required on refined slice // + if (interpolate == 1 ) { + InterpolateSliceValues( *smf, interp_lo, slice_cc_nd_box, dom_geom, + ncomp, nghost, slice_lo, slice_hi, SliceType, real_box); + } + + + if (coarsen == false) { + return smf; + } + else if ( coarsen == true ) { + Vector<BoxArray> crse_ba(1); + crse_ba[0] = sba[0]; + crse_ba[0].coarsen(slice_cr_ratio); + + AMREX_ALWAYS_ASSERT(crse_ba[0].size() == sba[0].size()); + + cs_mf.reset( new MultiFab(amrex::convert(crse_ba[0],SliceType), + sdmap[0], ncomp,nghost)); + + MultiFab& mfSrc = *smf; + MultiFab& mfDst = *cs_mf; + + MFIter mfi_dst(mfDst); + for (MFIter mfi(mfSrc); mfi.isValid(); ++mfi) { + + FArrayBox& Src_fabox = mfSrc[mfi]; + + const Box& Dst_bx = mfi_dst.validbox(); + FArrayBox& Dst_fabox = mfDst[mfi_dst]; + + int scomp = 0; + int dcomp = 0; + + IntVect cctype(AMREX_D_DECL(0,0,0)); + if( SliceType==cctype ) { + amrex::amrex_avgdown(Dst_bx, Dst_fabox, Src_fabox, dcomp, scomp, + ncomp, slice_cr_ratio); + } + IntVect ndtype(AMREX_D_DECL(1,1,1)); + if( SliceType == ndtype ) { + amrex::amrex_avgdown_nodes(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio); + } + if( SliceType == WarpX::Ex_nodal_flag ) { + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 0); + } + if( SliceType == WarpX::Ey_nodal_flag) { + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 1); + } + if( SliceType == WarpX::Ez_nodal_flag ) { + amrex::amrex_avgdown_edges(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 2); + } + if( SliceType == WarpX::Bx_nodal_flag) { + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 0); + } + if( SliceType == WarpX::By_nodal_flag ) { + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 1); + } + if( SliceType == WarpX::Bz_nodal_flag ) { + amrex::amrex_avgdown_faces(Dst_bx, Dst_fabox, Src_fabox, dcomp, + scomp, ncomp, slice_cr_ratio, 2); + } + + if ( mfi_dst.isValid() ) { + ++mfi_dst; + } + + } + return cs_mf; + + } + amrex::Abort("Should not hit this return statement."); + return smf; +} + + +/* \brief + * This function modifies the slice input parameters under certain conditions. + * The coarsening ratio, slice_cr_ratio is modified if the input is not an exponent of 2. + * for example, if the coarsening ratio is 3, 5 or 6, which is not an exponent of 2, + * then the value of coarsening ratio is modified to the nearest exponent of 2. + * The default value for coarsening ratio is 1. + * slice_realbox.lo and slice_realbox.hi are set equal to the simulation domain lo and hi + * if for the user-input for the slice lo and hi coordinates are outside the domain. + * If the slice_realbox.lo and slice_realbox.hi coordinates do not align with the data + * points and the number of cells in that dimension is greater than 1, and if the extent of + * the slice in that dimension is not coarsenable, then the value lo and hi coordinates are + * shifted to the nearest coarsenable point to include some extra data points in the slice. + * If slice_realbox.lo==slice_realbox.hi, then that dimension has only one cell and no + * modifications are made to the value. If the lo and hi do not align with a data point, + * then it is flagged for interpolation. + * \param real_box a Real box defined for the underlying domain. + * \param slice_realbox a Real box for defining the slice dimension. + * \param slice_cc_nd_box a Real box for defining the modified lo and hi of the slice + * such that the coordinates align with the underlying data points. + * If the dimension is reduced to have only one cell, the slice_realbox is not modified and * instead the values are interpolated to the coordinate from the nearest data points. + * \param slice_cr_ratio contains values of the coarsening ratio which may be modified + * if the input values do not satisfy coarsenability conditions. + * \param slice_lo and slice_hi are the index values of the slice + * \param interp_lo are set to 0 or 1 if they are flagged for interpolation. + * The slice shares the same index space as that of the simulation domain. + */ + + +void +CheckSliceInput( const RealBox real_box, RealBox &slice_cc_nd_box, + RealBox &slice_realbox, IntVect &slice_cr_ratio, + Vector<Geometry> dom_geom, IntVect const SliceType, + IntVect &slice_lo, IntVect &slice_hi, IntVect &interp_lo) +{ + + IntVect slice_lo2(AMREX_D_DECL(0,0,0)); + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + // Modify coarsening ratio if the input value is not an exponent of 2 for AMR // + if ( slice_cr_ratio[idim] > 0 ) { + int log_cr_ratio = floor ( log2( double(slice_cr_ratio[idim]))); + slice_cr_ratio[idim] = exp2( double(log_cr_ratio) ); + } + + //// Default coarsening ratio is 1 // + // Modify lo if input is out of bounds // + if ( slice_realbox.lo(idim) < real_box.lo(idim) ) { + slice_realbox.setLo( idim, real_box.lo(idim)); + amrex::Print() << " slice lo is out of bounds. " << + " Modified it in dimension " << idim << + " to be aligned with the domain box\n"; + } + + // Modify hi if input in out od bounds // + if ( slice_realbox.hi(idim) > real_box.hi(idim) ) { + slice_realbox.setHi( idim, real_box.hi(idim)); + amrex::Print() << " slice hi is out of bounds." << + " Modified it in dimension " << idim << + " to be aligned with the domain box\n"; + } + + // Factor to ensure index values computation depending on index type // + double fac = ( 1.0 - SliceType[idim] )*dom_geom[0].CellSize(idim)*0.5; + // if dimension is reduced to one cell length // + if ( slice_realbox.hi(idim) - slice_realbox.lo(idim) <= 0) + { + slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) ); + slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) ); + + if ( slice_cr_ratio[idim] > 1) slice_cr_ratio[idim] = 1; + + // check for interpolation -- compute index lo with floor and ceil + if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) >= fac ) { + slice_lo[idim] = floor( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac ) ) + / dom_geom[0].CellSize(idim)) + fac * 1E-10); + slice_lo2[idim] = ceil( ( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) + fac) ) + / dom_geom[0].CellSize(idim)) - fac * 1E-10 ); + } + else { + slice_lo[idim] = round( (slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) + / dom_geom[0].CellSize(idim)); + slice_lo2[idim] = ceil((slice_cc_nd_box.lo(idim) + - (real_box.lo(idim) ) ) + / dom_geom[0].CellSize(idim) ); + } + + // flag for interpolation -- if reduced dimension location // + // does not align with data point // + if ( slice_lo[idim] == slice_lo2[idim]) { + if ( slice_cc_nd_box.lo(idim) - real_box.lo(idim) < fac ) { + interp_lo[idim] = 1; + } + } + else { + interp_lo[idim] = 1; + } + + // ncells = 1 if dimension is reduced // + slice_hi[idim] = slice_lo[idim] + 1; + } + else + { + // moving realbox.lo and reabox.hi to nearest coarsenable grid point // + int index_lo = floor(((slice_realbox.lo(idim) + 1E-10 + - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); + int index_hi = ceil(((slice_realbox.hi(idim) - 1E-10 + - (real_box.lo(idim))) / dom_geom[0].CellSize(idim))); + + bool modify_cr = true; + + while ( modify_cr == true) { + int lo_new = index_lo; + int hi_new = index_hi; + int mod_lo = index_lo % slice_cr_ratio[idim]; + int mod_hi = index_hi % slice_cr_ratio[idim]; + modify_cr = false; + + // To ensure that the index.lo is coarsenable // + if ( mod_lo > 0) { + lo_new = index_lo - mod_lo; + } + // To ensure that the index.hi is coarsenable // + if ( mod_hi > 0) { + hi_new = index_hi + (slice_cr_ratio[idim] - mod_hi); + } + + // If modified index.hi is > baselinebox.hi, move the point // + // to the previous coarsenable point // + if ( (hi_new * dom_geom[0].CellSize(idim)) + > real_box.hi(idim) - real_box.lo(idim) + dom_geom[0].CellSize(idim)*0.01 ) + { + hi_new = index_hi - mod_hi; + } + + if ( (hi_new - lo_new) == 0 ){ + amrex::Print() << " Diagnostic Warning :: "; + amrex::Print() << " Coarsening ratio "; + amrex::Print() << slice_cr_ratio[idim] << " in dim "<< idim; + amrex::Print() << "is leading to zero cells for slice."; + amrex::Print() << " Thus reducing cr_ratio by half.\n"; + + slice_cr_ratio[idim] = slice_cr_ratio[idim]/2; + modify_cr = true; + } + + if ( modify_cr == false ) { + index_lo = lo_new; + index_hi = hi_new; + } + slice_lo[idim] = index_lo; + slice_hi[idim] = index_hi - 1; // since default is cell-centered + } + slice_realbox.setLo( idim, index_lo * dom_geom[0].CellSize(idim) + + real_box.lo(idim) ); + slice_realbox.setHi( idim, index_hi * dom_geom[0].CellSize(idim) + + real_box.lo(idim) ); + slice_cc_nd_box.setLo( idim, slice_realbox.lo(idim) + fac ); + slice_cc_nd_box.setHi( idim, slice_realbox.hi(idim) - fac ); + } + } +} + + +/* \brief + * This function is called if the coordinates of the slice do not align with data points + * \param interp_lo is an IntVect which is flagged as 1, if interpolation + is required in the dimension. + */ +void +InterpolateSliceValues(MultiFab& smf, IntVect interp_lo, RealBox slice_realbox, + Vector<Geometry> geom, int ncomp, int nghost, + IntVect slice_lo, IntVect slice_hi, IntVect SliceType, + const RealBox real_box) +{ + for (MFIter mfi(smf); mfi.isValid(); ++mfi) + { + const Box& bx = mfi.tilebox(); + const auto IndType = smf.ixType(); + const auto lo = amrex::lbound(bx); + const auto hi = amrex::ubound(bx); + FArrayBox& fabox = smf[mfi]; + + for ( int idim = 0; idim < AMREX_SPACEDIM; ++idim) { + if ( interp_lo[idim] == 1 ) { + InterpolateLo( bx, fabox, slice_lo, geom, idim, SliceType, + slice_realbox, 0, ncomp, nghost, real_box); + } + } + } + +} + +void +InterpolateLo(const Box& bx, FArrayBox &fabox, IntVect slice_lo, + Vector<Geometry> geom, int idir, IntVect IndType, + RealBox slice_realbox, int srccomp, int ncomp, + int nghost, const RealBox real_box ) +{ + auto fabarr = fabox.array(); + const auto lo = amrex::lbound(bx); + const auto hi = amrex::ubound(bx); + double fac = ( 1.0-IndType[idir] )*geom[0].CellSize(idir) * 0.5; + int imin = slice_lo[idir]; + double minpos = imin*geom[0].CellSize(idir) + fac + real_box.lo(idir); + double maxpos = (imin+1)*geom[0].CellSize(idir) + fac + real_box.lo(idir); + double slice_minpos = slice_realbox.lo(idir) ; + + switch (idir) { + case 0: + { + if ( imin >= lo.x && imin <= lo.x) { + for (int n = srccomp; n < srccomp + ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i+1,j,k,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + case 1: + { + if ( imin >= lo.y && imin <= lo.y) { + for (int n = srccomp; n < srccomp+ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i,j+1,k,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + case 2: + { + if ( imin >= lo.z && imin <= lo.z) { + for (int n = srccomp; n < srccomp+ncomp; ++n) { + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + double minval = fabarr(i,j,k,n); + double maxval = fabarr(i,j,k+1,n); + double ratio = (maxval - minval) / (maxpos - minpos); + double xdiff = slice_minpos - minpos; + double newval = minval + xdiff * ratio; + fabarr(i,j,k,n) = newval; + } + } + } + } + } + break; + } + + } + +} + + + + + + + diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 6d646dd42..7e7ed278d 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -11,6 +11,8 @@ #include <AMReX_AmrMeshInSituBridge.H> #endif +#include "SliceDiagnostic.H" + #ifdef AMREX_USE_ASCENT #include <ascent.hpp> #include <AMReX_Conduit_Blueprint.H> @@ -771,3 +773,94 @@ WarpX::WriteJobInfo (const std::string& dir) const jobInfoFile.close(); } } + + +/* \brief + * The slice is ouput using visMF and can be visualized used amrvis. + */ +void +WarpX::WriteSlicePlotFile () const +{ + if (F_fp[0] ) { + VisMF::Write( (*F_slice[0]), "vismf_F_slice"); + } + + if (rho_fp[0]) { + VisMF::Write( (*rho_slice[0]), "vismf_rho_slice"); + } + + VisMF::Write( (*Efield_slice[0][0]), amrex::Concatenate("vismf_Ex_slice_",istep[0])); + VisMF::Write( (*Efield_slice[0][1]), amrex::Concatenate("vismf_Ey_slice_",istep[0])); + VisMF::Write( (*Efield_slice[0][2]), amrex::Concatenate("vismf_Ez_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][0]), amrex::Concatenate("vismf_Bx_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][1]), amrex::Concatenate("vismf_By_slice_",istep[0])); + VisMF::Write( (*Bfield_slice[0][2]), amrex::Concatenate("vismf_Bz_slice_",istep[0])); + VisMF::Write( (*current_slice[0][0]), amrex::Concatenate("vismf_jx_slice_",istep[0])); + VisMF::Write( (*current_slice[0][1]), amrex::Concatenate("vismf_jy_slice_",istep[0])); + VisMF::Write( (*current_slice[0][2]), amrex::Concatenate("vismf_jz_slice_",istep[0])); + +} + + +void +WarpX::InitializeSliceMultiFabs () +{ + + int nlevels = Geom().size(); + + F_slice.resize(nlevels); + rho_slice.resize(nlevels); + current_slice.resize(nlevels); + Efield_slice.resize(nlevels); + Bfield_slice.resize(nlevels); + +} + + +// To generate slice that inherits index type of underlying data // +void +WarpX::SliceGenerationForDiagnostics () +{ + + Vector<Geometry> dom_geom; + dom_geom = Geom(); + + if (F_fp[0] ) { + F_slice[0] = CreateSlice( *F_fp[0].get(), dom_geom, slice_realbox, + slice_cr_ratio ); + } + if (rho_fp[0]) { + rho_slice[0] = CreateSlice( *rho_fp[0].get(), dom_geom, slice_realbox, + slice_cr_ratio ); + } + + for (int idim = 0; idim < 3; ++idim) { + Efield_slice[0][idim] = CreateSlice( *Efield_fp[0][idim].get(), + dom_geom, slice_realbox, slice_cr_ratio ); + Bfield_slice[0][idim] = CreateSlice( *Bfield_fp[0][idim].get(), + dom_geom, slice_realbox, slice_cr_ratio ); + current_slice[0][idim] = CreateSlice( *current_fp[0][idim].get(), + dom_geom, slice_realbox, slice_cr_ratio ); + } + + +} + + +void +WarpX::ClearSliceMultiFabs () +{ + + F_slice.clear(); + rho_slice.clear(); + current_slice.clear(); + Efield_slice.clear(); + Bfield_slice.clear(); + F_slice.shrink_to_fit(); + rho_slice.shrink_to_fit(); + current_slice.shrink_to_fit(); + Efield_slice.shrink_to_fit(); + Bfield_slice.shrink_to_fit(); + +} + diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index ad7c7d840..a5d68e4f9 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -132,6 +132,8 @@ WarpX::EvolveEM (int numsteps) bool to_make_plot = (plot_int > 0) && ((step+1) % plot_int == 0); + // slice generation // + bool to_make_slice_plot = (slice_plot_int > 0) && ( (step+1)% slice_plot_int == 0); bool do_insitu = ((step+1) >= insitu_start) && (insitu_int > 0) && ((step+1) % insitu_int == 0); @@ -176,7 +178,8 @@ WarpX::EvolveEM (int numsteps) myBFD->writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); } - if (to_make_plot || do_insitu) + // slice gen // + if (to_make_plot || do_insitu || to_make_slice_plot) { FillBoundaryE(); FillBoundaryB(); @@ -192,11 +195,19 @@ WarpX::EvolveEM (int numsteps) last_insitu_step = step+1; if (to_make_plot) - WritePlotFile(); + WritePlotFile(); + + if (to_make_slice_plot) + { + InitializeSliceMultiFabs (); + SliceGenerationForDiagnostics(); + WriteSlicePlotFile(); + ClearSliceMultiFabs (); + } if (do_insitu) UpdateInSitu(); - } + } if (check_int > 0 && (step+1) % check_int == 0) { last_check_file_step = step+1; diff --git a/Source/Filter/BilinearFilter.H b/Source/Filter/BilinearFilter.H index 9bade3c77..4fcfec14b 100644 --- a/Source/Filter/BilinearFilter.H +++ b/Source/Filter/BilinearFilter.H @@ -1,30 +1,17 @@ -#include <AMReX_MultiFab.H> -#include <AMReX_CudaContainers.H> +#include <Filter.H> -#ifndef WARPX_FILTER_H_ -#define WARPX_FILTER_H_ +#ifndef WARPX_BILIN_FILTER_H_ +#define WARPX_BILIN_FILTER_H_ -class BilinearFilter +class BilinearFilter : public Filter { public: - BilinearFilter () = default; + + BilinearFilter() = default; void ComputeStencils(); - void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); amrex::IntVect npass_each_dir; - amrex::IntVect stencil_length_each_dir; - - // public for cuda - void Filter(const amrex::Box& tbx, - amrex::Array4<amrex::Real const> const& tmp, - amrex::Array4<amrex::Real > const& dst, - int scomp, int dcomp, int ncomp); - -private: - - amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z; - amrex::Dim3 slen; }; -#endif +#endif // #ifndef WARPX_BILIN_FILTER_H_ diff --git a/Source/Filter/BilinearFilter.cpp b/Source/Filter/BilinearFilter.cpp index f6acaa5df..68ebde687 100644 --- a/Source/Filter/BilinearFilter.cpp +++ b/Source/Filter/BilinearFilter.cpp @@ -68,173 +68,3 @@ void BilinearFilter::ComputeStencils(){ slen.z = 1; #endif } - - -#ifdef AMREX_USE_CUDA - -void -BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) -{ - BL_PROFILE("BilinearFilter::ApplyStencil()"); - ncomp = std::min(ncomp, srcmf.nComp()); - - for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) - { - const auto& src = srcmf.array(mfi); - const auto& dst = dstmf.array(mfi); - const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - - // tmpfab has enough ghost cells for the stencil - FArrayBox tmp_fab(gbx,ncomp); - Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early - auto const& tmp = tmp_fab.array(); - - // Copy values in srcfab into tmpfab - const Box& ibx = gbx & srcmf[mfi].box(); - AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, - { - if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { - tmp(i,j,k,n) = src(i,j,k,n+scomp); - } else { - tmp(i,j,k,n) = 0.0; - } - }); - - // Apply filter - Filter(tbx, tmp, dst, 0, dcomp, ncomp); - } -} - -void BilinearFilter::Filter (const Box& tbx, - Array4<Real const> const& tmp, - Array4<Real > const& dst, - int scomp, int dcomp, int ncomp) -{ - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); - Dim3 slen_local = slen; - AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, - { - Real d = 0.0; - - for (int iz=0; iz < slen_local.z; ++iz){ - for (int iy=0; iy < slen_local.y; ++iy){ - for (int ix=0; ix < slen_local.x; ++ix){ -#if (AMREX_SPACEDIM == 3) - Real sss = sx[ix]*sy[iy]*sz[iz]; -#else - Real sss = sx[ix]*sz[iy]; -#endif -#if (AMREX_SPACEDIM == 3) - d += sss*( tmp(i-ix,j-iy,k-iz,scomp+n) - +tmp(i+ix,j-iy,k-iz,scomp+n) - +tmp(i-ix,j+iy,k-iz,scomp+n) - +tmp(i+ix,j+iy,k-iz,scomp+n) - +tmp(i-ix,j-iy,k+iz,scomp+n) - +tmp(i+ix,j-iy,k+iz,scomp+n) - +tmp(i-ix,j+iy,k+iz,scomp+n) - +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else - d += sss*( tmp(i-ix,j-iy,k,scomp+n) - +tmp(i+ix,j-iy,k,scomp+n) - +tmp(i-ix,j+iy,k,scomp+n) - +tmp(i+ix,j+iy,k,scomp+n)); -#endif - } - } - } - - dst(i,j,k,dcomp+n) = d; - }); -} - -#else - -void -BilinearFilter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) -{ - BL_PROFILE("BilinearFilter::ApplyStencil()"); - ncomp = std::min(ncomp, srcmf.nComp()); -#ifdef _OPENMP -#pragma omp parallel -#endif - { - FArrayBox tmpfab; - for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ - const auto& srcfab = srcmf[mfi]; - auto& dstfab = dstmf[mfi]; - const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - // tmpfab has enough ghost cells for the stencil - tmpfab.resize(gbx,ncomp); - tmpfab.setVal(0.0, gbx, 0, ncomp); - // Copy values in srcfab into tmpfab - const Box& ibx = gbx & srcfab.box(); - tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - // Apply filter - Filter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); - } - } -} - -void BilinearFilter::Filter (const Box& tbx, - Array4<Real const> const& tmp, - Array4<Real > const& dst, - int scomp, int dcomp, int ncomp) -{ - const auto lo = amrex::lbound(tbx); - const auto hi = amrex::ubound(tbx); - // tmp and dst are of type Array4 (Fortran ordering) - amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); - amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); - amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); - for (int n = 0; n < ncomp; ++n) { - // Set dst value to 0. - for (int k = lo.z; k <= hi.z; ++k) { - for (int j = lo.y; j <= hi.y; ++j) { - for (int i = lo.x; i <= hi.x; ++i) { - dst(i,j,k,dcomp+n) = 0.0; - } - } - } - // 3 nested loop on 3D stencil - for (int iz=0; iz < slen.z; ++iz){ - for (int iy=0; iy < slen.y; ++iy){ - for (int ix=0; ix < slen.x; ++ix){ -#if (AMREX_SPACEDIM == 3) - Real sss = sx[ix]*sy[iy]*sz[iz]; -#else - Real sss = sx[ix]*sz[iy]; -#endif - // 3 nested loop on 3D array - for (int k = lo.z; k <= hi.z; ++k) { - for (int j = lo.y; j <= hi.y; ++j) { - AMREX_PRAGMA_SIMD - for (int i = lo.x; i <= hi.x; ++i) { -#if (AMREX_SPACEDIM == 3) - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n) - +tmp(i+ix,j-iy,k-iz,scomp+n) - +tmp(i-ix,j+iy,k-iz,scomp+n) - +tmp(i+ix,j+iy,k-iz,scomp+n) - +tmp(i-ix,j-iy,k+iz,scomp+n) - +tmp(i+ix,j-iy,k+iz,scomp+n) - +tmp(i-ix,j+iy,k+iz,scomp+n) - +tmp(i+ix,j+iy,k+iz,scomp+n)); -#else - dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) - +tmp(i+ix,j-iy,k,scomp+n) - +tmp(i-ix,j+iy,k,scomp+n) - +tmp(i+ix,j+iy,k,scomp+n)); -#endif - } - } - } - } - } - } - } -} - -#endif diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H new file mode 100644 index 000000000..eaa0498c9 --- /dev/null +++ b/Source/Filter/Filter.H @@ -0,0 +1,43 @@ +#include <AMReX_MultiFab.H> +#include <AMReX_CudaContainers.H> + +#ifndef WARPX_FILTER_H_ +#define WARPX_FILTER_H_ + +class Filter +{ +public: + Filter () = default; + + // Apply stencil on MultiFab. + // Guard cells are handled inside this function + void ApplyStencil(amrex::MultiFab& dstmf, + const amrex::MultiFab& srcmf, int scomp=0, + int dcomp=0, int ncomp=10000); + + // Apply stencil on a FabArray. + void ApplyStencil (amrex::FArrayBox& dstfab, + const amrex::FArrayBox& srcfab, const amrex::Box& tbx, + int scomp=0, int dcomp=0, int ncomp=10000); + + // public for cuda + void DoFilter(const amrex::Box& tbx, + amrex::Array4<amrex::Real const> const& tmp, + amrex::Array4<amrex::Real > const& dst, + int scomp, int dcomp, int ncomp); + + // In 2D, stencil_length_each_dir = {length(stencil_x), length(stencil_z)} + amrex::IntVect stencil_length_each_dir; + +protected: + // Stencil along each direction. + // in 2D, stencil_y is not initialized. + amrex::Gpu::ManagedVector<amrex::Real> stencil_x, stencil_y, stencil_z; + // Length of each stencil. + // In 2D, slen = {length(stencil_x), length(stencil_z), 1} + amrex::Dim3 slen; + +private: + +}; +#endif // #ifndef WARPX_FILTER_H_ diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp new file mode 100644 index 000000000..f8a2dd6c2 --- /dev/null +++ b/Source/Filter/Filter.cpp @@ -0,0 +1,257 @@ +#include <WarpX.H> +#include <Filter.H> + +#ifdef _OPENMP +#include <omp.h> +#endif + +using namespace amrex; + +#ifdef AMREX_USE_CUDA + +/* \brief Apply stencil on MultiFab (GPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \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) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(MultiFab)"); + ncomp = std::min(ncomp, srcmf.nComp()); + + for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) + { + const auto& src = srcmf.array(mfi); + const auto& dst = dstmf.array(mfi); + const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); + + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcmf[mfi].box(); + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + { + if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + tmp(i,j,k,n) = src(i,j,k,n+scomp); + } else { + tmp(i,j,k,n) = 0.0; + } + }); + + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); + } +} + +/* \brief Apply stencil on FArrayBox (GPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \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) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + const auto& src = srcfab.array(); + const auto& dst = dstfab.array(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); + + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + { + if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + tmp(i,j,k,n) = src(i,j,k,n+scomp); + } else { + tmp(i,j,k,n) = 0.0; + } + }); + + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); +} + +/* \brief Apply stencil (2D/3D, CPU/GPU) + */ +void Filter::DoFilter (const Box& tbx, + Array4<Real const> const& tmp, + Array4<Real > const& dst, + int scomp, int dcomp, int ncomp) +{ + amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); + amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); + amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + Dim3 slen_local = slen; + AMREX_PARALLEL_FOR_4D ( tbx, ncomp, i, j, k, n, + { + Real d = 0.0; + + for (int iz=0; iz < slen_local.z; ++iz){ + for (int iy=0; iy < slen_local.y; ++iy){ + for (int ix=0; ix < slen_local.x; ++ix){ +#if (AMREX_SPACEDIM == 3) + Real sss = sx[ix]*sy[iy]*sz[iz]; +#else + Real sss = sx[ix]*sz[iy]; +#endif +#if (AMREX_SPACEDIM == 3) + d += sss*( tmp(i-ix,j-iy,k-iz,scomp+n) + +tmp(i+ix,j-iy,k-iz,scomp+n) + +tmp(i-ix,j+iy,k-iz,scomp+n) + +tmp(i+ix,j+iy,k-iz,scomp+n) + +tmp(i-ix,j-iy,k+iz,scomp+n) + +tmp(i+ix,j-iy,k+iz,scomp+n) + +tmp(i-ix,j+iy,k+iz,scomp+n) + +tmp(i+ix,j+iy,k+iz,scomp+n)); +#else + d += sss*( tmp(i-ix,j-iy,k,scomp+n) + +tmp(i+ix,j-iy,k,scomp+n) + +tmp(i-ix,j+iy,k,scomp+n) + +tmp(i+ix,j+iy,k,scomp+n)); +#endif + } + } + } + + dst(i,j,k,dcomp+n) = d; + }); +} + +#else + +/* \brief Apply stencil on MultiFab (CPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \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) +{ + BL_PROFILE("BilinearFilter::ApplyStencil()"); + ncomp = std::min(ncomp, srcmf.nComp()); +#ifdef _OPENMP +#pragma omp parallel +#endif + { + FArrayBox tmpfab; + for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ + const auto& srcfab = srcmf[mfi]; + auto& dstfab = dstmf[mfi]; + const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); + } + } +} + +/* \brief Apply stencil on FArrayBox (CPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \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) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + FArrayBox tmpfab; + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + const Box& ibx = gbx & srcfab.box(); + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); +} + +void Filter::DoFilter (const Box& tbx, + Array4<Real const> const& tmp, + Array4<Real > const& dst, + int scomp, int dcomp, int ncomp) +{ + const auto lo = amrex::lbound(tbx); + const auto hi = amrex::ubound(tbx); + // tmp and dst are of type Array4 (Fortran ordering) + amrex::Real const* AMREX_RESTRICT sx = stencil_x.data(); + amrex::Real const* AMREX_RESTRICT sy = stencil_y.data(); + amrex::Real const* AMREX_RESTRICT sz = stencil_z.data(); + for (int n = 0; n < ncomp; ++n) { + // Set dst value to 0. + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + for (int i = lo.x; i <= hi.x; ++i) { + dst(i,j,k,dcomp+n) = 0.0; + } + } + } + // 3 nested loop on 3D stencil + for (int iz=0; iz < slen.z; ++iz){ + for (int iy=0; iy < slen.y; ++iy){ + for (int ix=0; ix < slen.x; ++ix){ +#if (AMREX_SPACEDIM == 3) + Real sss = sx[ix]*sy[iy]*sz[iz]; +#else + Real sss = sx[ix]*sz[iy]; +#endif + // 3 nested loop on 3D array + for (int k = lo.z; k <= hi.z; ++k) { + for (int j = lo.y; j <= hi.y; ++j) { + AMREX_PRAGMA_SIMD + for (int i = lo.x; i <= hi.x; ++i) { +#if (AMREX_SPACEDIM == 3) + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k-iz,scomp+n) + +tmp(i+ix,j-iy,k-iz,scomp+n) + +tmp(i-ix,j+iy,k-iz,scomp+n) + +tmp(i+ix,j+iy,k-iz,scomp+n) + +tmp(i-ix,j-iy,k+iz,scomp+n) + +tmp(i+ix,j-iy,k+iz,scomp+n) + +tmp(i-ix,j+iy,k+iz,scomp+n) + +tmp(i+ix,j+iy,k+iz,scomp+n)); +#else + dst(i,j,k,dcomp+n) += sss*(tmp(i-ix,j-iy,k,scomp+n) + +tmp(i+ix,j-iy,k,scomp+n) + +tmp(i-ix,j+iy,k,scomp+n) + +tmp(i+ix,j+iy,k,scomp+n)); +#endif + } + } + } + } + } + } + } +} + +#endif // #ifdef AMREX_USE_CUDA diff --git a/Source/Filter/Make.package b/Source/Filter/Make.package index 36e74bfb0..8b8e9b50b 100644 --- a/Source/Filter/Make.package +++ b/Source/Filter/Make.package @@ -1,5 +1,9 @@ +CEXE_headers += Filter.H +CEXE_sources += Filter.cpp CEXE_sources += BilinearFilter.cpp CEXE_headers += BilinearFilter.H +CEXE_sources += NCIGodfreyFilter.cpp +CEXE_headers += NCIGodfreyFilter.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Filter VPATH_LOCATIONS += $(WARPX_HOME)/Source/Filter diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H new file mode 100644 index 000000000..a53039dfa --- /dev/null +++ b/Source/Filter/NCIGodfreyFilter.H @@ -0,0 +1,30 @@ +#include <Filter.H> + +#ifndef WARPX_GODFREY_FILTER_H_ +#define WARPX_GODFREY_FILTER_H_ + +enum class godfrey_coeff_set { Ex_Ey_Bz=0, Bx_By_Ez=1 }; + +class NCIGodfreyFilter : public Filter +{ +public: + + NCIGodfreyFilter () = default; + + NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_); + + void ComputeStencils(); + + void getGodfreyCoeffs(godfrey_coeff_set coeff_set_in); + + static constexpr int stencil_width = 4; + +private: + + godfrey_coeff_set coeff_set; + amrex::Real cdtodz; + int l_lower_order_in_v; + +}; + +#endif // #ifndef WARPX_GODFREY_FILTER_H_ diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp new file mode 100644 index 000000000..3f589260a --- /dev/null +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -0,0 +1,78 @@ +#include <WarpX.H> +#include <NCIGodfreyFilter.H> +#include <NCIGodfreyTables.H> + +#ifdef _OPENMP +#include <omp.h> +#endif + +using namespace amrex; + +NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_){ + // Store parameters into class data members + coeff_set = coeff_set_; + cdtodz = cdtodz_; + l_lower_order_in_v = l_lower_order_in_v_; + // NCI Godfrey filter has fixed size, and is applied along z only. +#if (AMREX_SPACEDIM == 3) + stencil_length_each_dir = {1,1,5}; + slen = {1,1,5}; +#else + stencil_length_each_dir = {1,5}; + slen = {1,5,1}; +#endif +} + +void NCIGodfreyFilter::ComputeStencils(){ + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( +#if ( AMREX_SPACEDIM == 3 ) + slen.z==5, +#else + slen.y==5, +#endif + "ERROR: NCI filter requires 5 points in z"); + AMREX_ALWAYS_ASSERT_WITH_MESSAGE(l_lower_order_in_v==1, + "ERROR: NCI corrector requires l_lower_order_in_v=1, i.e., Galerkin scheme"); + + // Interpolate coefficients from the table, and store into prestencil. + int index = tab_length*cdtodz; + index = min(index, tab_length-2); + index = max(index, 0); + Real weight_right = cdtodz - index/tab_length; + Real prestencil[4]; + for(int i=0; i<tab_width; i++){ + if (coeff_set == godfrey_coeff_set::Ex_Ey_Bz) { + prestencil[i] = (1-weight_right)*table_nci_godfrey_Ex_Ey_Bz[index ][i] + + weight_right*table_nci_godfrey_Ex_Ey_Bz[index+1][i]; + } else if (coeff_set == godfrey_coeff_set::Bx_By_Ez) { + prestencil[i] = (1-weight_right)*table_nci_godfrey_Bx_By_Ez[index ][i] + + weight_right*table_nci_godfrey_Bx_By_Ez[index+1][i]; + } + } + + // Compute stencil_z + stencil_z.resize( 5 ); + stencil_z[0] = (256 + 128*prestencil[0] + 96*prestencil[1] + 80*prestencil[2] + 70*prestencil[3]) / 256; + stencil_z[1] = -( 64*prestencil[0] + 64*prestencil[1] + 60*prestencil[2] + 56*prestencil[3]) / 256; + stencil_z[2] = ( 16*prestencil[1] + 24*prestencil[2] + 28*prestencil[3]) / 256; + stencil_z[3] = -( 4*prestencil[2] + 8*prestencil[3]) / 256; + stencil_z[4] = ( 1*prestencil[3]) / 256; + + // Compute stencil_x and stencil_y (no filter in these directions, + // so only 1 coeff, equal to 1) + stencil_x.resize(1); + stencil_x[0] = 1.; +#if (AMREX_SPACEDIM == 3) + stencil_y.resize(1); + stencil_y[0] = 1.; +#endif + + // Due to the way Filter::DoFilter() is written, + // coefficient 0 has to be /2 + stencil_x[0] /= 2.; +#if (AMREX_SPACEDIM == 3) + stencil_y[0] /= 2.; +#endif + stencil_z[0] /= 2.; + +} diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 52996a60a..1053ace89 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -39,10 +39,6 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_3d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_3d -#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs -#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_3d - - #elif (AMREX_SPACEDIM == 2) #define WRPX_COMPUTE_DIVB warpx_compute_divb_2d @@ -66,9 +62,6 @@ #define WRPX_PUSH_LEAPFROG warpx_push_leapfrog_2d #define WRPX_PUSH_LEAPFROG_POSITIONS warpx_push_leapfrog_positions_2d -#define WRPX_PXR_NCI_CORR_INIT init_godfrey_filter_coeffs -#define WRPX_PXR_GODFREY_FILTER apply_godfrey_filter_z_2d - #ifdef WARPX_RZ #define WRPX_COMPUTE_DIVE warpx_compute_dive_rz #else @@ -455,14 +448,6 @@ extern "C" const BL_FORT_FAB_ARG_ANYD(fine), const int* ncomp); - void WRPX_PXR_NCI_CORR_INIT(amrex::Real*, amrex::Real*, const int, - const amrex::Real, const int); - - void WRPX_PXR_GODFREY_FILTER (const int* lo, const int* hi, - amrex_real* fou, const int* olo, const int* ohi, - const amrex_real* fin, const int* ilo, const int* ihi, - const amrex_real* stencil, const int* nsten); - #ifdef WARPX_USE_PSATD void warpx_fft_mpi_init (int fcomm); void warpx_fft_domain_decomp (int* warpx_local_nz, int* warpx_local_z0, diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index c17e8861b..12d541b08 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -117,7 +117,7 @@ contains pxr_ll4symtry = ll4symtry .eq. 1 pxr_l_lower_order_in_v = l_lower_order_in_v .eq. 1 pxr_l_nodal = l_nodal .eq. 1 - + exg_nguards = ixyzmin - exg_lo eyg_nguards = ixyzmin - eyg_lo ezg_nguards = ixyzmin - ezg_lo @@ -218,6 +218,11 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n IF ((nox.eq.1).and.(noy.eq.1).and.(noz.eq.1)) THEN CALL depose_rho_vecHVv2_1_1_1(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& nxguard,nyguard,nzguard,lvect) + + ELSE IF ((nox.eq.2).and.(noy.eq.2).and.(noz.eq.2)) THEN + CALL depose_rho_vecHVv2_2_2_2(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& + nxguard,nyguard,nzguard,lvect) + ELSE CALL pxr_depose_rho_n(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,nx,ny,nz,& nxguard,nyguard,nzguard,nox,noy,noz, & @@ -268,7 +273,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n #ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 #endif - + ! Compute the number of valid cells and guard cells integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) rho_nvalid = rho_ntot - 2*rho_ng @@ -387,7 +392,7 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) real(amrex_real), intent(IN) :: rmin, dr -#ifdef WARPX_RZ +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 #endif ! Compute the number of valid cells and guard cells diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index 23637ec97..0f33d1a0f 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -7,6 +7,7 @@ #include <WarpX.H> #include <WarpX_f.H> #include <BilinearFilter.H> +#include <NCIGodfreyFilter.H> #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> @@ -161,8 +162,6 @@ WarpX::InitNCICorrector () { if (WarpX::use_fdtd_nci_corr) { - mypc->fdtd_nci_stencilz_ex.resize(max_level+1); - mypc->fdtd_nci_stencilz_by.resize(max_level+1); for (int lev = 0; lev <= max_level; ++lev) { const Geometry& gm = Geom(lev); @@ -174,10 +173,15 @@ WarpX::InitNCICorrector () dz = dx[1]; } cdtodz = PhysConst::c * dt[lev] / dz; - WRPX_PXR_NCI_CORR_INIT( (mypc->fdtd_nci_stencilz_ex)[lev].data(), - (mypc->fdtd_nci_stencilz_by)[lev].data(), - mypc->nstencilz_fdtd_nci_corr, cdtodz, - WarpX::l_lower_order_in_v); + + // Initialize Godfrey filters + // Same filter for fields Ex, Ey and Bz + nci_godfrey_filter_exeybz[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Ex_Ey_Bz, cdtodz, WarpX::l_lower_order_in_v) ); + // Same filter for fields Bx, By and Ez + nci_godfrey_filter_bxbyez[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Bx_By_Ez, cdtodz, WarpX::l_lower_order_in_v) ); + // Compute Godfrey filters stencils + nci_godfrey_filter_exeybz[lev]->ComputeStencils(); + nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); } } } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 513f30b99..aaf7ead0e 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -181,17 +181,6 @@ public: void UpdateContinuousInjectionPosition(amrex::Real dt) const; int doContinuousInjection() const; - // - // Parameters for the Cherenkov corrector in the FDTD solver. - // Both stencils are calculated ar runtime. - // - // Number of coefficients for the stencil of the NCI corrector. - // The stencil is applied in the z direction only. - static constexpr int nstencilz_fdtd_nci_corr=5; - - amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_ex; - amrex::Vector<amrex::Array<amrex::Real, nstencilz_fdtd_nci_corr> > fdtd_nci_stencilz_by; - std::vector<std::string> GetSpeciesNames() const { return species_names; } PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 6d618c096..9d39ec2f9 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -8,8 +8,6 @@ using namespace amrex; -constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr; - MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) { ReadParameters(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 212084e64..ab37bf6ff 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1132,8 +1132,9 @@ PhysicalParticleContainer::Evolve (int lev, const std::array<Real,3>& dx = WarpX::CellSize(lev); const std::array<Real,3>& cdx = WarpX::CellSize(std::max(lev-1,0)); - const auto& mypc = WarpX::GetInstance().GetPartContainer(); - const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; + // Get instances of NCI Godfrey filters + const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; + const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; BL_ASSERT(OnSameGrids(lev,jx)); @@ -1165,7 +1166,7 @@ PhysicalParticleContainer::Evolve (int lev, { Real wt = amrex::second(); - const Box& box = pti.validbox(); + const Box& box = pti.validbox(); auto& attribs = pti.GetAttribs(); @@ -1182,7 +1183,7 @@ PhysicalParticleContainer::Evolve (int lev, const long np = pti.numParticles(); - // Data on the grid + // Data on the grid FArrayBox const* exfab = &(Ex[pti]); FArrayBox const* eyfab = &(Ey[pti]); FArrayBox const* ezfab = &(Ez[pti]); @@ -1190,6 +1191,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* byfab = &(By[pti]); FArrayBox const* bzfab = &(Bz[pti]); + Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1201,54 +1203,43 @@ PhysicalParticleContainer::Evolve (int lev, static_cast<int>(WarpX::noz)}); #endif - // both 2d and 3d + // Filter Ex (Both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD(Ex[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box()); + // Update exfab reference exfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD(Ez[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box()); ezfab = &filtered_Ez; + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD(By[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box()); byfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD(Ey[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box()); eyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD(Bx[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box()); bxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD(Bz[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box()); bzfab = &filtered_Bz; #endif } @@ -1424,53 +1415,43 @@ PhysicalParticleContainer::Evolve (int lev, static_cast<int>(WarpX::noz)}); #endif - // both 2d and 3d + // Filter Ex (both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD((*cEx)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + // Safeguard for GPU + exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box()); + // Update exfab reference cexfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD((*cEz)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box()); cezfab = &filtered_Ez; + + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD((*cBy)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box()); cbyfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD((*cEy)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], filtered_Ey.box()); ceyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD((*cBx)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], filtered_Bx.box()); cbxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD((*cBz)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box()); cbzfab = &filtered_Bz; #endif } diff --git a/Source/Utils/Make.package b/Source/Utils/Make.package index d8e2d2dab..cd335dbcf 100644 --- a/Source/Utils/Make.package +++ b/Source/Utils/Make.package @@ -4,6 +4,9 @@ CEXE_sources += WarpXTagging.cpp CEXE_sources += WarpXUtil.cpp CEXE_headers += WarpXConst.H CEXE_headers += WarpXUtil.H +CEXE_headers += WarpXAlgorithmSelection.H +CEXE_sources += WarpXAlgorithmSelection.cpp +CEXE_headers += NCIGodfreyTables.H INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/Utils VPATH_LOCATIONS += $(WARPX_HOME)/Source/Utils diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H new file mode 100644 index 000000000..8cb105aa0 --- /dev/null +++ b/Source/Utils/NCIGodfreyTables.H @@ -0,0 +1,224 @@ +#include <AMReX_AmrCore.H> + +#ifndef WARPX_GODFREY_COEFF_TABLE_H_ +#define WARPX_GODFREY_COEFF_TABLE_H_ + +// Table width. This is related to the stencil length +const int tab_width = 4; +// table length. Each line correspond to 1 value of cdtodz +// (here 101 values). +const int tab_length = 101; + +// Table of coefficient for Ex, Ey abd Bz +// We typically interpolate between two lines +const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ + -2.47536,2.04288,-0.598163,0.0314711, + -2.47536,2.04288,-0.598163,0.0314711, + -2.47545,2.04309,-0.598307,0.0315029, + -2.4756,2.04342,-0.598549,0.0315558, + -2.47581,2.0439,-0.598886,0.0316298, + -2.47608,2.0445,-0.59932,0.031725, + -2.47641,2.04525,-0.59985,0.0318412, + -2.4768,2.04612,-0.600477,0.0319785, + -2.47725,2.04714,-0.6012,0.0321367, + -2.47776,2.04829,-0.602019,0.0323158, + -2.47833,2.04957,-0.602934,0.0325158, + -2.47896,2.05099,-0.603944,0.0327364, + -2.47965,2.05254,-0.605051,0.0329777, + -2.4804,2.05423,-0.606253,0.0332396, + -2.48121,2.05606,-0.60755,0.0335218, + -2.48208,2.05802,-0.608942,0.0338243, + -2.48301,2.06012,-0.610429,0.0341469, + -2.48401,2.06235,-0.61201,0.0344895, + -2.48506,2.06471,-0.613685,0.0348519, + -2.48618,2.06721,-0.615453,0.0352339, + -2.48735,2.06984,-0.617314,0.0356353, + -2.48859,2.07261,-0.619268,0.0360559, + -2.48988,2.0755,-0.621312,0.0364954, + -2.49123,2.07853,-0.623447,0.0369536, + -2.49265,2.08169,-0.625672,0.0374302, + -2.49412,2.08498,-0.627986,0.0379248, + -2.49565,2.0884,-0.630386,0.0384372, + -2.49724,2.09194,-0.632873,0.0389669, + -2.49888,2.09561,-0.635443,0.0395135, + -2.50058,2.09939,-0.638096,0.0400766, + -2.50234,2.1033,-0.640829,0.0406557, + -2.50415,2.10732,-0.64364,0.0412502, + -2.50601,2.11145,-0.646526,0.0418594, + -2.50791,2.1157,-0.649485,0.0424828, + -2.50987,2.12004,-0.652512,0.0431196, + -2.51187,2.12448,-0.655604,0.0437688, + -2.51392,2.12901,-0.658756,0.0444297, + -2.516,2.13363,-0.661964,0.0451011, + -2.51812,2.13832,-0.665221,0.0457818, + -2.52027,2.14308,-0.668521,0.0464705, + -2.52244,2.14789,-0.671856,0.0471658, + -2.52464,2.15274,-0.675218,0.0478658, + -2.52684,2.15762,-0.678596,0.0485687, + -2.52906,2.16251,-0.68198,0.0492723, + -2.53126,2.16738,-0.685355,0.049974, + -2.53345,2.17222,-0.688706,0.0506708, + -2.53561,2.177,-0.692015,0.0513594, + -2.53773,2.18168,-0.69526,0.0520359, + -2.53978,2.18623,-0.698416,0.0526955, + -2.54175,2.19059,-0.701452,0.053333, + -2.5436,2.19471,-0.704331,0.0539417, + -2.54531,2.19852,-0.70701,0.0545141, + -2.54683,2.20193,-0.709433,0.0550409, + -2.5481,2.20483,-0.711533,0.0555106, + -2.54906,2.20709,-0.713224,0.0559094, + -2.54963,2.20852,-0.714397,0.0562198, + -2.54968,2.20888,-0.714907,0.0564196, + -2.54905,2.20785,-0.714562,0.0564797, + -2.54751,2.20496,-0.713094,0.0563618, + -2.54472,2.19955,-0.710118,0.0560124, + -2.54014,2.19058,-0.705048,0.0553544, + -2.53286,2.1763,-0.69693,0.0542684, + -2.52115,2.15344,-0.684027,0.05255, + -2.50098,2.11466,-0.66255,0.0497817, + -2.45797,2.03459,-0.620099,0.0446889, + -2.28371,1.72254,-0.465905,0.0283268, + -2.4885,2.04899,-0.599292,0.0390466, + -2.1433,1.36735,-0.220924,-0.00215633, + -2.4943,2.07019,-0.610552,0.035166, + -2.84529,2.77303,-1.00018,0.0724884, + -2.72242,2.51888,-0.847226,0.0509964, + -2.65633,2.3744,-0.750392,0.0326366, + -2.59601,2.23412,-0.646421,0.00868027, + -2.51477,2.0369,-0.491066,-0.0306397, + -2.35935,1.65155,-0.178971,-0.112713, + -1.84315,0.361693,0.876104,-0.393844, + -2.65422,2.39262,-0.789663,0.0516265, + -3.46529,4.42354,-2.45543,0.497097, + -3.15747,3.65311,-1.824,0.328432, + -3.04694,3.37613,-1.59668,0.267631, + -2.99205,3.23814,-1.48302,0.237103, + -2.96075,3.15894,-1.41733,0.219317, + -2.94172,3.11028,-1.37649,0.20811, + -2.92994,3.07962,-1.35025,0.200755, + -2.92283,3.06054,-1.33338,0.195859, + -2.91894,3.04938,-1.3229,0.192637, + -2.91736,3.04394,-1.31702,0.190612, + -2.91753,3.04278,-1.31456,0.189477, + -2.91905,3.04494,-1.31475,0.189026, + -2.92165,3.04973,-1.31705,0.189117, + -2.92512,3.05667,-1.32105,0.189646, + -2.92933,3.06539,-1.32646,0.190538, + -2.93416,3.07562,-1.33308,0.191735, + -2.93952,3.08715,-1.34072,0.193194, + -2.94535,3.09982,-1.34925,0.194881, + -2.95159,3.11349,-1.35858,0.196769, + -2.9582,3.12805,-1.36861,0.198838, + -2.96514,3.14342,-1.37929,0.201068, + -2.97239,3.15953,-1.39055,0.203448, + -2.97991,3.17632,-1.40234,0.205964, + -2.98769,3.19374,-1.41463,0.208607 + }; + +// Table of coefficient for Bx, By and Ez +// We typically interpolate between two lines +const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ + -2.80862,2.80104,-1.14615,0.154077, + -2.80862,2.80104,-1.14615,0.154077, + -2.80851,2.80078,-1.14595,0.154027, + -2.80832,2.80034,-1.14561,0.153945, + -2.80807,2.79973,-1.14514,0.153829, + -2.80774,2.79894,-1.14454,0.15368, + -2.80733,2.79798,-1.1438,0.153498, + -2.80685,2.79685,-1.14292,0.153284, + -2.8063,2.79554,-1.14192,0.153036, + -2.80568,2.79405,-1.14077,0.152756, + -2.80498,2.79239,-1.1395,0.152443, + -2.80421,2.79056,-1.13809,0.152098, + -2.80337,2.78856,-1.13656,0.151721, + -2.80246,2.78638,-1.13488,0.151312, + -2.80147,2.78404,-1.13308,0.150871, + -2.80041,2.78152,-1.13115,0.150397, + -2.79927,2.77882,-1.12908,0.149893, + -2.79807,2.77596,-1.12689,0.149356, + -2.79679,2.77292,-1.12456,0.148789, + -2.79543,2.76972,-1.12211,0.14819, + -2.79401,2.76634,-1.11953,0.14756, + -2.79251,2.76279,-1.11681,0.1469, + -2.79094,2.75907,-1.11397,0.146208, + -2.78929,2.75517,-1.111,0.145486, + -2.78757,2.7511,-1.10789,0.144733, + -2.78578,2.74686,-1.10466,0.14395, + -2.78391,2.74245,-1.1013,0.143137, + -2.78196,2.73786,-1.09781,0.142293, + -2.77994,2.73309,-1.09419,0.141419, + -2.77784,2.72814,-1.09043,0.140514, + -2.77566,2.72301,-1.08654,0.139578, + -2.7734,2.7177,-1.08252,0.138612, + -2.77106,2.7122,-1.07836,0.137614, + -2.76864,2.70651,-1.07406,0.136586, + -2.76613,2.70062,-1.06962,0.135525, + -2.76353,2.69453,-1.06503,0.134432, + -2.76084,2.68824,-1.0603,0.133307, + -2.75806,2.68173,-1.05541,0.132148, + -2.75518,2.675,-1.05037,0.130954, + -2.75219,2.66804,-1.04516,0.129725, + -2.7491,2.66084,-1.03978,0.12846, + -2.7459,2.65339,-1.03423,0.127156, + -2.74257,2.64566,-1.02848,0.125813, + -2.73912,2.63765,-1.02254,0.124428, + -2.73552,2.62934,-1.01638,0.122999, + -2.73178,2.62069,-1.01,0.121523, + -2.72787,2.61169,-1.00337,0.119996, + -2.72379,2.6023,-0.996479,0.118417, + -2.71951,2.59248,-0.989294,0.116778, + -2.71501,2.58218,-0.981786,0.115076, + -2.71026,2.57135,-0.97392,0.113303, + -2.70524,2.55991,-0.965651,0.111453, + -2.69989,2.54778,-0.956922,0.109514, + -2.69416,2.53484,-0.947666,0.107476, + -2.68799,2.52096,-0.937795,0.105324, + -2.68129,2.50596,-0.927197,0.103039, + -2.67394,2.48959,-0.915724,0.100597, + -2.66578,2.47153,-0.903179,0.097968, + -2.65657,2.4513,-0.889283,0.0951084, + -2.64598,2.42824,-0.873638,0.0919592, + -2.63347,2.40127,-0.855632,0.0884325, + -2.61813,2.36864,-0.834261,0.0843898, + -2.59821,2.32701,-0.807691,0.0795876, + -2.56971,2.26887,-0.77188,0.0735132, + -2.51823,2.16823,-0.713448,0.0645399, + -2.33537,1.8294,-0.533852,0.0409941, + -2.53143,2.14818,-0.670502,0.053982, + -2.17737,1.43641,-0.259095,0.00101255, + -2.51929,2.12931,-0.654743,0.0452381, + -2.86122,2.82221,-1.05039,0.0894636, + -2.72908,2.54506,-0.87834,0.0626188, + -2.6536,2.37954,-0.7665,0.0409117, + -2.58374,2.21923,-0.649738,0.0146791, + -2.49284,2.00346,-0.48457,-0.0255348, + -2.32762,1.60337,-0.1698,-0.105287, + -1.80149,0.316787,0.855414,-0.369652, + -2.60242,2.28418,-0.721378,0.040091, + -3.40335,4.25157,-2.29817,0.449834, + -3.0852,3.47341,-1.67791,0.28982, + -2.9642,3.17856,-1.44399,0.229852, + -2.89872,3.01966,-1.31861,0.197945, + -2.85668,2.91811,-1.23894,0.17783, + -2.82679,2.84621,-1.18287,0.163785, + -2.80401,2.79167,-1.14058,0.153278, + -2.78577,2.74819,-1.10706,0.145015, + -2.77061,2.7122,-1.07946,0.138267, + -2.75764,2.68152,-1.05606,0.132589, + -2.74627,2.65475,-1.03575,0.127695, + -2.73612,2.63093,-1.01777,0.123395, + -2.72692,2.6094,-1.00159,0.119553, + -2.71846,2.58968,-0.986841,0.116074, + -2.71061,2.57142,-0.973239,0.112887, + -2.70323,2.55434,-0.960573,0.109937, + -2.69626,2.53824,-0.948678,0.107185, + -2.68962,2.52294,-0.937429,0.104598, + -2.68327,2.50833,-0.926722,0.102151, + -2.67714,2.4943,-0.916477,0.0998223, + -2.67122,2.48076,-0.906627,0.0975966, + -2.66546,2.46764,-0.897118,0.0954599, + -2.65985,2.45489,-0.887903,0.0934011, + -2.65437,2.44244,-0.878945,0.0914107 + }; + +#endif // #ifndef WARPX_GODFREY_COEFF_TABLE_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H new file mode 100644 index 000000000..3fb23698a --- /dev/null +++ b/Source/Utils/WarpXAlgorithmSelection.H @@ -0,0 +1,57 @@ +#ifndef UTILS_WARPXALGORITHMSELECTION_H_ +#define UTILS_WARPXALGORITHMSELECTION_H_ + +#include <AMReX_ParmParse.H> +#include <string> + +struct MaxwellSolverAlgo { + // These numbers corresponds to the algorithm code in WarpX's + // `warpx_push_bvec` and `warpx_push_evec_f` + enum { + Yee = 0, + CKC = 1 + }; +}; + +struct ParticlePusherAlgo { + // These numbers correspond to the algorithm code in WarpX's + // `warpx_particle_pusher` + enum { + Boris = 0, + Vay = 1 + }; +}; + +struct CurrentDepositionAlgo { + enum { + // These numbers corresponds to the algorithm code in PICSAR's + // `depose_jxjyjz_generic` and `depose_jxjyjz_generic_2d` + Direct = 3, + DirectVectorized = 2, + EsirkepovNonOptimized = 1, + Esirkepov = 0 + }; +}; + +struct ChargeDepositionAlgo { + // These numbers corresponds to the algorithm code in WarpX's + // `warpx_charge_deposition` function + enum { + Vectorized = 0, + Standard = 1 + }; +}; + +struct GatheringAlgo { + // These numbers corresponds to the algorithm code in PICSAR's + // `geteb3d_energy_conserving_generic` function + enum { + Vectorized = 0, + Standard = 1 + }; +}; + +int +GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ); + +#endif // UTILS_WARPXALGORITHMSELECTION_H_ diff --git a/Source/Utils/WarpXAlgorithmSelection.cpp b/Source/Utils/WarpXAlgorithmSelection.cpp new file mode 100644 index 000000000..89802e064 --- /dev/null +++ b/Source/Utils/WarpXAlgorithmSelection.cpp @@ -0,0 +1,94 @@ +#include <WarpXAlgorithmSelection.H> +#include <map> +#include <algorithm> + +// Define dictionary with correspondance between user-input strings, +// and corresponding integer for use inside the code (e.g. in PICSAR). + +const std::map<std::string, int> maxwell_solver_algo_to_int = { + {"yee", MaxwellSolverAlgo::Yee }, +#ifndef WARPX_RZ // Not available in RZ + {"ckc", MaxwellSolverAlgo::CKC }, +#endif + {"default", MaxwellSolverAlgo::Yee } +}; + +const std::map<std::string, int> particle_pusher_algo_to_int = { + {"boris", ParticlePusherAlgo::Boris }, + {"vay", ParticlePusherAlgo::Vay }, + {"default", ParticlePusherAlgo::Boris } +}; + +const std::map<std::string, int> current_deposition_algo_to_int = { + {"esirkepov", CurrentDepositionAlgo::Esirkepov }, + {"direct", CurrentDepositionAlgo::Direct }, +#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D + {"direct-vectorized", CurrentDepositionAlgo::DirectVectorized }, +#endif + {"default", CurrentDepositionAlgo::Esirkepov } +}; + +const std::map<std::string, int> charge_deposition_algo_to_int = { + {"standard", ChargeDepositionAlgo::Standard }, +#if (!defined AMREX_USE_GPU)&&(AMREX_SPACEDIM == 3) // Only available on CPU and 3D + {"vectorized", ChargeDepositionAlgo::Vectorized }, + {"default", ChargeDepositionAlgo::Vectorized } +#else + {"default", ChargeDepositionAlgo::Standard } +#endif +}; + +const std::map<std::string, int> gathering_algo_to_int = { + {"standard", GatheringAlgo::Standard }, +#ifndef AMREX_USE_GPU // Only available on CPU + {"vectorized", GatheringAlgo::Vectorized }, + {"default", GatheringAlgo::Vectorized } +#else + {"default", GatheringAlgo::Standard } +#endif +}; + + +int +GetAlgorithmInteger( amrex::ParmParse& pp, const char* pp_search_key ){ + + // Read user input ; use "default" if it is not found + std::string algo = "default"; + pp.query( pp_search_key, algo ); + // Convert to lower case + std::transform(algo.begin(), algo.end(), algo.begin(), ::tolower); + + // Pick the right dictionary + std::map<std::string, int> algo_to_int; + if (pp_search_key == "maxwell_fdtd_solver") { + algo_to_int = maxwell_solver_algo_to_int; + } else if (pp_search_key == "particle_pusher") { + algo_to_int = particle_pusher_algo_to_int; + } else if (pp_search_key == "current_deposition") { + algo_to_int = current_deposition_algo_to_int; + } else if (pp_search_key == "charge_deposition") { + algo_to_int = charge_deposition_algo_to_int; + } else if (pp_search_key == "field_gathering") { + algo_to_int = gathering_algo_to_int; + } else { + std::string pp_search_string = pp_search_key; + amrex::Abort("Unknown algorithm type: " + pp_search_string); + } + + // Check if the user-input is a valid key for the dictionary + if (algo_to_int.count(algo) == 0){ + // Not a valid key ; print error message + std::string pp_search_string = pp_search_key; + std::string error_message = "Invalid string for algo." + pp_search_string + + ": " + algo + ".\nThe valid values are:\n"; + for ( const auto &valid_pair : algo_to_int ) { + if (valid_pair.first != "default"){ + error_message += " - " + valid_pair.first + "\n"; + } + } + amrex::Abort(error_message); + } + + // If the input is a valid key, return the value + return algo_to_int[algo]; +} diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index a0ab1f26f..ae781f9aa 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -65,6 +65,26 @@ WarpX::MoveWindow (bool move_j) ResetProbDomain(RealBox(new_lo, new_hi)); + // Moving slice coordinates - lo and hi - with moving window // + // slice box is modified only if slice diagnostics is initialized in input // + if ( slice_plot_int > 0 ) + { + Real new_slice_lo[AMREX_SPACEDIM]; + Real new_slice_hi[AMREX_SPACEDIM]; + const Real* current_slice_lo = slice_realbox.lo(); + const Real* current_slice_hi = slice_realbox.hi(); + for ( int i = 0; i < AMREX_SPACEDIM; i++) { + new_slice_lo[i] = current_slice_lo[i]; + new_slice_hi[i] = current_slice_hi[i]; + } + int num_shift_base_slice = static_cast<int> ((moving_window_x - + current_slice_lo[dir]) / cdx[dir]); + new_slice_lo[dir] = current_slice_lo[dir] + num_shift_base_slice*cdx[dir]; + new_slice_hi[dir] = current_slice_hi[dir] + num_shift_base_slice*cdx[dir]; + slice_realbox.setLo(new_slice_lo); + slice_realbox.setHi(new_slice_hi); + } + int num_shift = num_shift_base; int num_shift_crse = num_shift; diff --git a/Source/WarpX.H b/Source/WarpX.H index 35b072142..2580c5320 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -23,6 +23,7 @@ #include <PML.H> #include <BoostedFrameDiagnostic.H> #include <BilinearFilter.H> +#include <NCIGodfreyFilter.H> #ifdef WARPX_USE_PSATD #include <fftw3.h> @@ -146,7 +147,9 @@ public: static amrex::IntVect filter_npass_each_dir; BilinearFilter bilinear_filter; - + amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_exeybz; + amrex::Vector< std::unique_ptr<NCIGodfreyFilter> > nci_godfrey_filter_bxbyez; + static int num_mirrors; amrex::Vector<amrex::Real> mirror_z; amrex::Vector<amrex::Real> mirror_z_width; @@ -242,6 +245,12 @@ public: static int do_moving_window; static int moving_window_dir; + // slice generation // + void InitializeSliceMultiFabs (); + void SliceGenerationForDiagnostics (); + void WriteSlicePlotFile () const; + void ClearSliceMultiFabs (); + protected: //! Tagging cells for refinement @@ -553,6 +562,17 @@ private: bool is_synchronized = true; + //Slice Parameters + int slice_max_grid_size; + int slice_plot_int = -1; + amrex::RealBox slice_realbox; + amrex::IntVect slice_cr_ratio; + amrex::Vector< std::unique_ptr<amrex::MultiFab> > F_slice; + amrex::Vector< std::unique_ptr<amrex::MultiFab> > rho_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_slice; + amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Bfield_slice; + #ifdef WARPX_USE_PSATD // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) // This includes data for the FFT guard cells (between FFT groups) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 3d7f7dcc5..c2cf97f30 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -17,6 +17,7 @@ #include <WarpXConst.H> #include <WarpXWrappers.h> #include <WarpXUtil.H> +#include <WarpXAlgorithmSelection.H> #ifdef BL_USE_SENSEI_INSITU #include <AMReX_AmrMeshInSituBridge.H> @@ -35,11 +36,11 @@ Vector<int> WarpX::boost_direction = {0,0,0}; int WarpX::do_compute_max_step_from_zmax = 0; Real WarpX::zmax_plasma_to_compute_max_step = 0.; -long WarpX::current_deposition_algo = 3; -long WarpX::charge_deposition_algo = 0; -long WarpX::field_gathering_algo = 1; -long WarpX::particle_pusher_algo = 0; -int WarpX::maxwell_fdtd_solver_id = 0; +long WarpX::current_deposition_algo; +long WarpX::charge_deposition_algo; +long WarpX::field_gathering_algo; +long WarpX::particle_pusher_algo; +int WarpX::maxwell_fdtd_solver_id; long WarpX::nox = 1; long WarpX::noy = 1; @@ -118,7 +119,7 @@ WarpX::ResetInstance () { delete m_instance; m_instance = nullptr; -} +} WarpX::WarpX () { @@ -228,6 +229,11 @@ WarpX::WarpX () #ifdef BL_USE_SENSEI_INSITU insitu_bridge = nullptr; #endif + + // NCI Godfrey filters can have different stencils + // at different levels (the stencil depends on c*dt/dz) + nci_godfrey_filter_exeybz.resize(nlevs_max); + nci_godfrey_filter_bxbyez.resize(nlevs_max); } WarpX::~WarpX () @@ -276,10 +282,10 @@ WarpX::ReadParameters () ReadBoostedFrameParameters(gamma_boost, beta_boost, boost_direction); - // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is + // pp.query returns 1 if argument zmax_plasma_to_compute_max_step is // specified by the user, 0 otherwise. - do_compute_max_step_from_zmax = - pp.query("zmax_plasma_to_compute_max_step", + do_compute_max_step_from_zmax = + pp.query("zmax_plasma_to_compute_max_step", zmax_plasma_to_compute_max_step); pp.queryarr("B_external", B_external); @@ -318,7 +324,7 @@ WarpX::ReadParameters () "gamma_boost must be > 1 to use the boosted frame diagnostic."); pp.query("lab_data_directory", lab_data_directory); - + std::string s; pp.get("boost_direction", s); AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (s == "z" || s == "Z"), @@ -476,29 +482,12 @@ WarpX::ReadParameters () } { - ParmParse pp("algo"); - pp.query("current_deposition", current_deposition_algo); - pp.query("charge_deposition", charge_deposition_algo); - pp.query("field_gathering", field_gathering_algo); - pp.query("particle_pusher", particle_pusher_algo); - std::string s_solver = ""; - pp.query("maxwell_fdtd_solver", s_solver); - std::transform(s_solver.begin(), - s_solver.end(), - s_solver.begin(), - ::tolower); - // if maxwell_fdtd_solver is specified, set the value - // of maxwell_fdtd_solver_id accordingly. - // Otherwise keep the default value maxwell_fdtd_solver_id=0 - if (s_solver != "") { - if (s_solver == "yee") { - maxwell_fdtd_solver_id = 0; - } else if (s_solver == "ckc") { - maxwell_fdtd_solver_id = 1; - } else { - amrex::Abort("Unknown FDTD Solver type " + s_solver); - } - } + ParmParse pp("algo"); + current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); + charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition"); + field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); + particle_pusher_algo = GetAlgorithmInteger(pp, "particle_pusher"); + maxwell_fdtd_solver_id = GetAlgorithmInteger(pp, "maxwell_fdtd_solver"); } #ifdef WARPX_USE_PSATD @@ -527,6 +516,34 @@ WarpX::ReadParameters () pp.query("config", insitu_config); pp.query("pin_mesh", insitu_pin_mesh); } + + // for slice generation // + { + ParmParse pp("slice"); + amrex::Vector<Real> slice_lo(AMREX_SPACEDIM); + amrex::Vector<Real> slice_hi(AMREX_SPACEDIM); + Vector<int> slice_crse_ratio(AMREX_SPACEDIM); + // set default slice_crse_ratio // + for (int idim=0; idim < AMREX_SPACEDIM; ++idim ) + { + slice_crse_ratio[idim] = 1; + } + pp.queryarr("dom_lo",slice_lo,0,AMREX_SPACEDIM); + pp.queryarr("dom_hi",slice_hi,0,AMREX_SPACEDIM); + pp.queryarr("coarsening_ratio",slice_crse_ratio,0,AMREX_SPACEDIM); + pp.query("plot_int",slice_plot_int); + slice_realbox.setLo(slice_lo); + slice_realbox.setHi(slice_hi); + slice_cr_ratio = IntVect(AMREX_D_DECL(1,1,1)); + for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) + { + if (slice_crse_ratio[idim] > 1 ) { + slice_cr_ratio[idim] = slice_crse_ratio[idim]; + } + } + + } + } // This is a virtual function. @@ -628,7 +645,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number int ngz; if (WarpX::use_fdtd_nci_corr) { - int ng = ngz_tmp + (mypc->nstencilz_fdtd_nci_corr-1); + int ng = ngz_tmp + NCIGodfreyFilter::stencil_width; ngz = (ng % 2) ? ng+1 : ng; } else { ngz = ngz_nonci; @@ -833,8 +850,6 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm if (load_balance_int > 0) { costs[lev].reset(new MultiFab(ba, dm, 1, 0)); } - - } std::array<Real,3> |