From 1ce5957cb32a3b7ae39ad68b528e5a370854c58f Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 7 Jul 2019 11:14:31 -0700 Subject: avoid duplication in current deposition when using buffers --- Source/Particles/PhysicalParticleContainer.cpp | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7e7c9534e..1ed318c13 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1529,8 +1529,14 @@ PhysicalParticleContainer::Evolve (int lev, // // Current Deposition // - DepositCurrent(pti, wp, uxp, uyp, uzp, jx, jy, jz, - cjx, cjy, cjz, np_current, np, thread_num, lev, dt); + // Deposit inside domains + DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + 0, np_current, thread_num, + lev, lev, dt, jx.nGrow()); + // Deposit in buffers + DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + np_current, np-np_current, thread_num, + lev, lev-1, dt, jx.nGrow()); // // copy particle data back -- cgit v1.2.3 From 6d3848aa366d40a3b8f48131da172a6aef5ebae0 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 8 Jul 2019 16:01:37 -0700 Subject: bug when tiling on CPU --- Examples/Tests/SingleParticle/Untitled.ipynb | 114 +++++++++++++++++++++++++ Examples/Tests/SingleParticle/inputs | 10 +-- Source/Particles/PhysicalParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.cpp | 66 +++++++++++--- 4 files changed, 172 insertions(+), 20 deletions(-) create mode 100644 Examples/Tests/SingleParticle/Untitled.ipynb (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Tests/SingleParticle/Untitled.ipynb b/Examples/Tests/SingleParticle/Untitled.ipynb new file mode 100644 index 000000000..7a98be71a --- /dev/null +++ b/Examples/Tests/SingleParticle/Untitled.ipynb @@ -0,0 +1,114 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "# Import statements\n", + "import yt ; yt.funcs.mylog.setLevel(50)\n", + "import numpy as np\n", + "import scipy.constants as scc\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib notebook" + ] + }, + { + "cell_type": "code", + "execution_count": 82, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/mthevenet/anaconda3/lib/python3.7/site-packages/yt/units/yt_array.py:1394: RuntimeWarning: invalid value encountered in true_divide\n", + " out=out, **kwargs)\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ds = yt.load( './old_diags/plt00001' ) # Create a dataset object\n", + "sl = yt.SlicePlot(ds, 2, 'jz', aspect=.4) # Create a sliceplot object\n", + "sl.annotate_grids() # Show grids\n", + "sl.show() # Show the plot" + ] + }, + { + "cell_type": "code", + "execution_count": 83, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/mthevenet/anaconda3/lib/python3.7/site-packages/yt/units/yt_array.py:1394: RuntimeWarning: invalid value encountered in true_divide\n", + " out=out, **kwargs)\n" + ] + }, + { + "data": { + "text/html": [ + "
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ds = yt.load( './new_diags/plt00001' ) # Create a dataset object\n", + "sl = yt.SlicePlot(ds, 2, 'jx', aspect=.4) # Create a sliceplot object\n", + "sl.annotate_grids() # Show grids\n", + "sl.show() # Show the plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + " " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.1" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/Examples/Tests/SingleParticle/inputs b/Examples/Tests/SingleParticle/inputs index 89e3eb1c6..feecb996c 100644 --- a/Examples/Tests/SingleParticle/inputs +++ b/Examples/Tests/SingleParticle/inputs @@ -1,11 +1,11 @@ -max_step = 1 -amr.n_cell = 8 8 +max_step = 10 +amr.n_cell = 8 16 amr.max_level = 0 amr.blocking_factor = 8 amr.max_grid_size = 16 amr.plot_int = 1 geometry.coord_sys = 0 -geometry.is_periodic = 0 0 +geometry.is_periodic = 1 1 geometry.prob_lo = -8 -12 geometry.prob_hi = 8 12 warpx.do_pml = 0 @@ -23,6 +23,6 @@ particles.species_names = electron electron.charge = -q_e electron.mass = 1.e10 electron.injection_style = "SingleParticle" -electron.single_particle_pos = 0.0 0.0 0.0 -electron.single_particle_vel = 1.e20 0.0 1.e20 # gamma*beta +electron.single_particle_pos = 4 0 6 +electron.single_particle_vel = -2 0 -5 # -2. 1. 4. electron.single_particle_weight = 1.0 diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 43b46ec49..b78d36c87 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1530,6 +1530,7 @@ PhysicalParticleContainer::Evolve (int lev, // Current Deposition // // Deposit inside domains + Print()<<"np_current "<& xyzmin_tile = WarpX::LowerCorner(tilebox, depos_lev); - const std::array& xyzmin = xyzmin_tile; - Print()<<" xyzmin_tile "<& xyzmin_tile = WarpX::LowerCorner(tilebox, depos_lev); + // const std::array& xyzmin = xyzmin_tile; + const std::array& xyzmin = WarpX::LowerCorner(tilebox, depos_lev);; + + // Print()<<" xyzmin_tile "<array(pti); auto jy_arr = jy->array(pti); auto jz_arr = jz->array(pti); - - // auto const& jx_arr = jx->array(); - // auto const& jy_arr = jy->array(); - // auto const& jz_arr = jz->array(); + */ + + // const auto& jx_arr = local_jx[thread_num].array(); + // const auto& jy_arr = local_jy[thread_num].array(); + // const auto& jz_arr = local_jz[thread_num].array(); + + //Array4 const& jx_arr = local_jx[thread_num].array(tbx); + //Array4 const& jy_arr = local_jy[thread_num].array(tby); + //Array4 const& jz_arr = local_jz[thread_num].array(tbz); + Array4 const& jx_arr = local_jx[thread_num].array(); + Array4 const& jy_arr = local_jy[thread_num].array(); + Array4 const& jz_arr = local_jz[thread_num].array(); + + Dim3 lox = lbound(jx_arr); + Dim3 loy = lbound(jy_arr); + Dim3 loz = lbound(jz_arr); + /* + Dim3 loy = lbound(local_jy[thread_num]); + Dim3 loz = lbound(local_jz[thread_num]); + Dim3 hix = hbound(local_jx[thread_num]); + Dim3 hiy = hbound(local_jy[thread_num]); + Dim3 hiz = hbound(local_jz[thread_num]); + */ + Print()<<" lox "< Vx, Vy, Vz; //pti.GetPosition(Vx, Vy, Vz); @@ -446,7 +474,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, std::cout<<" vx "< Date: Tue, 9 Jul 2019 16:22:09 -0700 Subject: Changed warpx_copy_attribs function from F90 to CPP inside PhysicalParticleContainer class. --- .DS_Store | Bin 0 -> 6148 bytes Source/.DS_Store | Bin 0 -> 6148 bytes Source/FortranInterface/WarpX_f.F90 | 33 ------------ Source/FortranInterface/WarpX_f.H | 6 --- Source/Particles/PhysicalParticleContainer.H | 3 ++ Source/Particles/PhysicalParticleContainer.cpp | 56 ++++++++++++++------- .../Particles/RigidInjectedParticleContainer.cpp | 20 ++------ 7 files changed, 47 insertions(+), 71 deletions(-) create mode 100644 .DS_Store create mode 100644 Source/.DS_Store (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 000000000..b96111fac Binary files /dev/null and b/.DS_Store differ diff --git a/Source/.DS_Store b/Source/.DS_Store new file mode 100644 index 000000000..01640e062 Binary files /dev/null and b/Source/.DS_Store differ diff --git a/Source/FortranInterface/WarpX_f.F90 b/Source/FortranInterface/WarpX_f.F90 index 0560a9981..542cc95a1 100644 --- a/Source/FortranInterface/WarpX_f.F90 +++ b/Source/FortranInterface/WarpX_f.F90 @@ -8,39 +8,6 @@ module warpx_module contains - subroutine warpx_copy_attribs(np, xp, yp, zp, uxp, uyp, uzp, & - xpold, ypold, zpold, uxpold, uypold, uzpold) & - bind(c,name='warpx_copy_attribs') - integer(c_long), intent(in) :: np - real(amrex_real), intent(in) :: xp(np) - real(amrex_real), intent(in) :: yp(np) - real(amrex_real), intent(in) :: zp(np) - real(amrex_real), intent(in) :: uxp(np) - real(amrex_real), intent(in) :: uyp(np) - real(amrex_real), intent(in) :: uzp(np) - real(amrex_real), intent(inout) :: xpold(np) - real(amrex_real), intent(inout) :: ypold(np) - real(amrex_real), intent(inout) :: zpold(np) - real(amrex_real), intent(inout) :: uxpold(np) - real(amrex_real), intent(inout) :: uypold(np) - real(amrex_real), intent(inout) :: uzpold(np) - - integer n - - do n = 1, np - - xpold(n) = xp(n) - ypold(n) = yp(n) - zpold(n) = zp(n) - - uxpold(n) = uxp(n) - uypold(n) = uyp(n) - uzpold(n) = uzp(n) - - end do - - end subroutine warpx_copy_attribs - subroutine warpx_compute_E (lo, hi, & phi, phlo, phhi, & Ex, Exlo, Exhi, & diff --git a/Source/FortranInterface/WarpX_f.H b/Source/FortranInterface/WarpX_f.H index 98dd8685d..0440148eb 100644 --- a/Source/FortranInterface/WarpX_f.H +++ b/Source/FortranInterface/WarpX_f.H @@ -75,12 +75,6 @@ extern "C" { #endif - void warpx_copy_attribs(const long* np, - const amrex_real* xp, const amrex_real* yp, const amrex_real* zp, - const amrex_real* uxp, const amrex_real* uyp, const amrex_real* uzp, - amrex_real* xpold, amrex_real* ypold, amrex_real* zpold, - amrex_real* uxpold, amrex_real* uypold, amrex_real* uzpold); - // Charge deposition void warpx_charge_deposition(amrex::Real* rho, const long* np, const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp, const amrex::Real* w, diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index bd8cfb47e..99fc0f8da 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -77,6 +77,9 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; + + void warpx_copy_attribs(WarpXParIter& pti,const amrex::Real* xp, + const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 43b46ec49..69c27cd21 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1729,7 +1729,13 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, Real dt) { - // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + warpx_copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + } + + // The following attributes should be included in CPP version of warpx_particle_pusher + // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; @@ -1741,22 +1747,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, auto& Byp = attribs[PIdx::By]; auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& xpold = pti.GetAttribs(particle_comps["xold"]); - auto& ypold = pti.GetAttribs(particle_comps["yold"]); - auto& zpold = pti.GetAttribs(particle_comps["zold"]); - auto& uxpold = pti.GetAttribs(particle_comps["uxold"]); - auto& uypold = pti.GetAttribs(particle_comps["uyold"]); - auto& uzpold = pti.GetAttribs(particle_comps["uzold"]); - - warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), - uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); - } - + warpx_particle_pusher(&np, xp.dataPtr(), yp.dataPtr(), @@ -1870,6 +1861,37 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } +void PhysicalParticleContainer::warpx_copy_attribs(WarpXParIter& pti,const amrex::Real* xp, + const amrex::Real* yp, const amrex::Real* zp) +{ + auto& attribs = pti.GetAttribs(); + + Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); + Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); + Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); + + Real* AMREX_RESTRICT xpold = attribs[PIdx::nattribs].dataPtr(); + Real* AMREX_RESTRICT ypold = attribs[PIdx::nattribs+1].dataPtr(); + Real* AMREX_RESTRICT zpold = attribs[PIdx::nattribs+2].dataPtr(); + Real* AMREX_RESTRICT uxpold = attribs[PIdx::nattribs+3].dataPtr(); + Real* AMREX_RESTRICT uypold = attribs[PIdx::nattribs+4].dataPtr(); + Real* AMREX_RESTRICT uzpold = attribs[PIdx::nattribs+5].dataPtr(); + + const long np = pti.numParticles(); + + ParallelFor( np, + [=] AMREX_GPU_DEVICE (long i) { + xpold[i]=xp[i]; + ypold[i]=yp[i]; + zpold[i]=zp[i]; + + uxpold[i]=uxp[i]; + uypold[i]=uyp[i]; + uzpold[i]=uzp[i]; + } + ); +} + void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real z_old, const Real z_new, const Real t_boost, const Real t_lab, const Real dt, diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index 2a3e8dd0d..ca02d1458 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -211,6 +211,11 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, Real dt) { + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + warpx_copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + } + // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. auto& attribs = pti.GetAttribs(); auto& uxp = attribs[PIdx::ux]; @@ -224,21 +229,6 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& xpold = pti.GetAttribs(particle_comps["xold"]); - auto& ypold = pti.GetAttribs(particle_comps["yold"]); - auto& zpold = pti.GetAttribs(particle_comps["zold"]); - auto& uxpold = pti.GetAttribs(particle_comps["uxold"]); - auto& uypold = pti.GetAttribs(particle_comps["uyold"]); - auto& uzpold = pti.GetAttribs(particle_comps["uzold"]); - - warpx_copy_attribs(&np, xp.dataPtr(), yp.dataPtr(), zp.dataPtr(), - uxp.dataPtr(), uyp.dataPtr(), uzp.dataPtr(), - xpold.dataPtr(), ypold.dataPtr(), zpold.dataPtr(), - uxpold.dataPtr(), uypold.dataPtr(), uzpold.dataPtr()); - } - // Save the position and momenta, making copies Cuda::ManagedDeviceVector xp_save, yp_save, zp_save; RealVector uxp_save, uyp_save, uzp_save; -- cgit v1.2.3 From 9957d80bce7bf685f96ccffe15b420b2216d7f11 Mon Sep 17 00:00:00 2001 From: Diana Amorim Date: Tue, 9 Jul 2019 18:06:24 -0700 Subject: Attempt to fix segfault issue - Unsuccessful. Tested alternatives to . --- .DS_Store | Bin 6148 -> 0 bytes .gitignore | 1 + Source/Particles/PhysicalParticleContainer.cpp | 22 ++++++++++++---------- 3 files changed, 13 insertions(+), 10 deletions(-) delete mode 100644 .DS_Store (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index b96111fac..000000000 Binary files a/.DS_Store and /dev/null differ diff --git a/.gitignore b/.gitignore index 02bd9eafa..13147b8da 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +.DS_Store plt* chk* main?d.*.ex diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 69c27cd21..e8654e573 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1866,20 +1866,22 @@ void PhysicalParticleContainer::warpx_copy_attribs(WarpXParIter& pti,const amrex { auto& attribs = pti.GetAttribs(); - Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); - Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); - Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); + //Real* AMREX_RESTRICT - Real* AMREX_RESTRICT xpold = attribs[PIdx::nattribs].dataPtr(); - Real* AMREX_RESTRICT ypold = attribs[PIdx::nattribs+1].dataPtr(); - Real* AMREX_RESTRICT zpold = attribs[PIdx::nattribs+2].dataPtr(); - Real* AMREX_RESTRICT uxpold = attribs[PIdx::nattribs+3].dataPtr(); - Real* AMREX_RESTRICT uypold = attribs[PIdx::nattribs+4].dataPtr(); - Real* AMREX_RESTRICT uzpold = attribs[PIdx::nattribs+5].dataPtr(); + const amrex::Real* uxp = attribs[PIdx::ux].dataPtr(); + const amrex::Real* uyp = attribs[PIdx::uy].dataPtr(); + const amrex::Real* uzp = attribs[PIdx::uz].dataPtr(); + + amrex::Real* xpold = attribs[PIdx::nattribs].dataPtr(); + amrex::Real* ypold = attribs[PIdx::nattribs+1].dataPtr(); + amrex::Real* zpold = attribs[PIdx::nattribs+2].dataPtr(); + amrex::Real* uxpold = attribs[PIdx::nattribs+3].dataPtr(); + amrex::Real* uypold = attribs[PIdx::nattribs+4].dataPtr(); + amrex::Real* uzpold = attribs[PIdx::nattribs+5].dataPtr(); const long np = pti.numParticles(); - ParallelFor( np, + amrex::ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { xpold[i]=xp[i]; ypold[i]=yp[i]; -- cgit v1.2.3 From b69ea2ccf999d80c4f2bb5e295a5ae170bb58a34 Mon Sep 17 00:00:00 2001 From: Diana Amorim Date: Wed, 10 Jul 2019 13:33:48 -0700 Subject: Fixed segfault error. --- Source/Particles/PhysicalParticleContainer.cpp | 27 +++++++++++++------------- 1 file changed, 13 insertions(+), 14 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index e8654e573..c45f73f0f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1861,27 +1861,26 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::warpx_copy_attribs(WarpXParIter& pti,const amrex::Real* xp, - const amrex::Real* yp, const amrex::Real* zp) +void PhysicalParticleContainer::warpx_copy_attribs(WarpXParIter& pti,const Real* xp, + const Real* yp, const Real* zp) { + auto& attribs = pti.GetAttribs(); - //Real* AMREX_RESTRICT - - const amrex::Real* uxp = attribs[PIdx::ux].dataPtr(); - const amrex::Real* uyp = attribs[PIdx::uy].dataPtr(); - const amrex::Real* uzp = attribs[PIdx::uz].dataPtr(); + Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); + Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); + Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - amrex::Real* xpold = attribs[PIdx::nattribs].dataPtr(); - amrex::Real* ypold = attribs[PIdx::nattribs+1].dataPtr(); - amrex::Real* zpold = attribs[PIdx::nattribs+2].dataPtr(); - amrex::Real* uxpold = attribs[PIdx::nattribs+3].dataPtr(); - amrex::Real* uypold = attribs[PIdx::nattribs+4].dataPtr(); - amrex::Real* uzpold = attribs[PIdx::nattribs+5].dataPtr(); + Real* AMREX_RESTRICT xpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); + Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); + Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); + Real* AMREX_RESTRICT uxpold = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); + Real* AMREX_RESTRICT uypold = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); + Real* AMREX_RESTRICT uzpold = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); const long np = pti.numParticles(); - amrex::ParallelFor( np, + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { xpold[i]=xp[i]; ypold[i]=yp[i]; -- cgit v1.2.3 From fb25dde6d6770f01bdfee71e10137c5fe8a8e035 Mon Sep 17 00:00:00 2001 From: Diana Amorim Date: Wed, 10 Jul 2019 16:56:58 -0700 Subject: Changed name of function to copy_attribs() --- Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- Source/Particles/RigidInjectedParticleContainer.cpp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 99fc0f8da..d55764682 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -78,7 +78,7 @@ public: const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - void warpx_copy_attribs(WarpXParIter& pti,const amrex::Real* xp, + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c45f73f0f..addde5e37 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1731,7 +1731,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - warpx_copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); } // The following attributes should be included in CPP version of warpx_particle_pusher @@ -1861,7 +1861,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } } -void PhysicalParticleContainer::warpx_copy_attribs(WarpXParIter& pti,const Real* xp, +void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const Real* yp, const Real* zp) { diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index ca02d1458..9bd4cb4fc 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -213,7 +213,7 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - warpx_copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); + copy_attribs(pti, xp.dataPtr(), yp.dataPtr(), zp.dataPtr()); } // This wraps the call to warpx_particle_pusher so that inheritors can modify the call. -- cgit v1.2.3 From 5baf28e5115f450d4f2fd821c33df5a37054792d Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 11 Jul 2019 09:07:26 -0700 Subject: remove print statements --- Source/Particles/PhysicalParticleContainer.cpp | 1 - Source/Particles/WarpXParticleContainer.cpp | 4 ---- 2 files changed, 5 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b78d36c87..df08d77a7 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1530,7 +1530,6 @@ PhysicalParticleContainer::Evolve (int lev, // Current Deposition // // Deposit inside domains - Print()<<"np_current "< Date: Thu, 11 Jul 2019 14:26:38 -0700 Subject: option use_picsar_deposition --- Source/Particles/Deposition/CurrentDeposition.H | 3 ++- Source/Particles/PhysicalParticleContainer.cpp | 32 ++++++++++++++++++------- Source/Particles/WarpXParticleContainer.H | 15 ++++++++++++ Source/Particles/WarpXParticleContainer.cpp | 2 +- Source/WarpX.cpp | 6 ++--- 5 files changed, 44 insertions(+), 14 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/Deposition/CurrentDeposition.H b/Source/Particles/Deposition/CurrentDeposition.H index 006ff05fe..66db003a3 100644 --- a/Source/Particles/Deposition/CurrentDeposition.H +++ b/Source/Particles/Deposition/CurrentDeposition.H @@ -109,7 +109,7 @@ void doDepositionShapeN(Real* xp, Real* yp, Real* zp, const amrex::Real q) { BL_PROFILE("PICSAR::CurrentDeposition"); - + Real dxi = 1.0/dx[0]; Real dzi = 1.0/dx[2]; Real dts2dx = 0.5*dt*dxi; @@ -130,6 +130,7 @@ void doDepositionShapeN(Real* xp, Real* yp, Real* zp, compute_shape_factor_lo(shape_lo); Real clightsq = 1.0/PhysConst::c/PhysConst::c; + // Loop over particles and deposit into jx_arr, jy_arr and jz_arr ParallelFor( np_to_depose, [=] AMREX_GPU_DEVICE (long ip) { diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index df08d77a7..86cd0827f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1529,16 +1529,30 @@ PhysicalParticleContainer::Evolve (int lev, // // Current Deposition // - // Deposit inside domains - DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, - 0, np_current, thread_num, - lev, lev, dt); - if (has_buffer){ - // Deposit in buffers - DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, - np_current, np-np_current, thread_num, - lev, lev-1, dt); + if (WarpX::use_picsar_deposition) { + // Deposit inside domains + DepositCurrentFortran(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + 0, np_current, thread_num, + lev, lev, dt); + if (has_buffer){ + // Deposit in buffers + DepositCurrentFortran(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + np_current, np-np_current, thread_num, + lev, lev-1, dt); + } + } else { + // Deposit inside domains + DepositCurrent(pti, wp, uxp, uyp, uzp, &jx, &jy, &jz, + 0, np_current, thread_num, + lev, lev, dt); + if (has_buffer){ + // Deposit in buffers + DepositCurrent(pti, wp, uxp, uyp, uzp, cjx, cjy, cjz, + np_current, np-np_current, thread_num, + lev, lev-1, dt); + } } + // // copy particle data back // diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 6182ec5bd..662b2e1b8 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -189,6 +189,21 @@ public: int depos_lev, amrex::Real dt); + virtual void DepositCurrentFortran(WarpXParIter& pti, + RealVector& wp, + RealVector& uxp, + RealVector& uyp, + RealVector& uzp, + amrex::MultiFab* jx, + amrex::MultiFab* jy, + amrex::MultiFab* jz, + const long offset, + const long np_to_depose, + int thread_num, + int lev, + int depos_lev, + amrex::Real dt); + // If particles start outside of the domain, ContinuousInjection // makes sure that they are initialized when they enter the domain, and // NOT before. Virtual function, overriden by derived classes. diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 42b99efcb..ac123f9b6 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -435,7 +435,7 @@ WarpXParticleContainer::DepositCurrent(WarpXParIter& pti, const long offset, const long np_to_depose, int thread_num, int lev, int depos_lev, Real dt) -{ +{ AMREX_ALWAYS_ASSERT_WITH_MESSAGE((depos_lev==(lev-1)) || (depos_lev==(lev )), "Deposition buffers only work for lev-1"); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 24011fb2f..fd4840b8d 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -485,11 +485,11 @@ WarpX::ReadParameters () { ParmParse pp("algo"); - pp.query("use_picsar_deposition", use_picsar_deposition) + pp.query("use_picsar_deposition", use_picsar_deposition); current_deposition_algo = GetAlgorithmInteger(pp, "current_deposition"); - if (!use_picsar_algo){ + if (!use_picsar_deposition){ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( current_deposition_algo >= 2, - "if not use_picsar_deposition, cannot use Esirkepov deposition.") + "if not use_picsar_deposition, cannot use Esirkepov deposition."); } charge_deposition_algo = GetAlgorithmInteger(pp, "charge_deposition"); field_gathering_algo = GetAlgorithmInteger(pp, "field_gathering"); -- cgit v1.2.3 From cd7f25974133d0aef31a90af6d79dd36372fc3fd Mon Sep 17 00:00:00 2001 From: Diana Amorim Date: Fri, 12 Jul 2019 20:24:28 -0700 Subject: Corrected typo of xold incorrect repetition. --- .../laser_acceleration/inputs.2d.boost | 21 +++++++++++---------- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- 2 files changed, 13 insertions(+), 12 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost index 5fe4bfcf6..75a246ea6 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 = 1000 -# stop_time = 1.9e-12 +# max_step = 2700 +stop_time = 1.9e-12 amr.n_cell = 128 1024 amr.max_grid_size = 64 amr.blocking_factor = 32 @@ -19,8 +19,10 @@ geometry.prob_hi = 128.e-6 0.96e-6 ################################# warpx.verbose = 1 amrex.v = 1 -algo.current_deposition = direct -algo.particle_pusher = vay +algo.current_deposition = direct #2 +algo.charge_deposition = standard # 0 +algo.field_gathering = standard # 0 +algo.particle_pusher = vay # 1 algo.maxwell_fdtd_solver = ckc interpolation.nox = 3 interpolation.noy = 3 @@ -37,12 +39,11 @@ warpx.serialize_ics = 1 ################################# ####### BOOST PARAMETERS ######## ################################# -warpx.gamma_boost = 30. +warpx.gamma_boost = 10. warpx.boost_direction = z warpx.do_boosted_frame_diagnostic = 1 warpx.num_snapshots_lab = 7 warpx.dt_snapshots_lab = 1.6678204759907604e-12 -warpx.boosted_frame_diag_fields = Ex Ez By jz ################################# ############ PLASMA ############# @@ -60,11 +61,11 @@ electrons.momentum_distribution_type = "gaussian" electrons.xmin = -120.e-6 electrons.xmax = 120.e-6 electrons.zmin = 0.5e-3 -electrons.zmax = 1. +electrons.zmax = .0035 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.5e25 +electrons.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 electrons.do_continuous_injection = 1 ions.charge = q_e @@ -75,11 +76,11 @@ ions.momentum_distribution_type = "gaussian" ions.xmin = -120.e-6 ions.xmax = 120.e-6 ions.zmin = 0.5e-3 -ions.zmax = 1. +ions.zmax = .0035 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.5e25 +ions.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 ions.do_continuous_injection = 1 beam.charge = -q_e diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index addde5e37..d6c4973c3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1872,8 +1872,8 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); Real* AMREX_RESTRICT xpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); - Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["xold"]).dataPtr(); + Real* AMREX_RESTRICT ypold = pti.GetAttribs(particle_comps["yold"]).dataPtr(); + Real* AMREX_RESTRICT zpold = pti.GetAttribs(particle_comps["zold"]).dataPtr(); Real* AMREX_RESTRICT uxpold = pti.GetAttribs(particle_comps["uxold"]).dataPtr(); Real* AMREX_RESTRICT uypold = pti.GetAttribs(particle_comps["uyold"]).dataPtr(); Real* AMREX_RESTRICT uzpold = pti.GetAttribs(particle_comps["uzold"]).dataPtr(); -- cgit v1.2.3 From 7a450310547e94df56a7dd089b40bec47f2ba772 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 15 Jul 2019 17:59:13 -0700 Subject: fix indentation in physicalparticlecontainer.cpp, couldn't read --- Source/Particles/PhysicalParticleContainer.cpp | 588 ++++++++++++------------- 1 file changed, 294 insertions(+), 294 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 28d611020..d47a7b220 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -12,7 +12,7 @@ using namespace amrex; long PhysicalParticleContainer:: NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, - const RealBox& tile_realbox, const RealBox& particle_real_box) + const RealBox& tile_realbox, const RealBox& particle_real_box) { const int lev = 0; const Geometry& geom = Geom(lev); @@ -24,43 +24,43 @@ NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (do_continuous_injection) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) - Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; - Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; - Real z = overlap_corner[2] + (iv[2] + 0.5)*dx[2]; + Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; + Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; + Real z = overlap_corner[2] + (iv[2] + 0.5)*dx[2]; #elif ( AMREX_SPACEDIM == 2 ) - Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; - Real y = 0; - Real z = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; + Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; + Real y = 0; + Real z = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; #endif - fac = GetRefineFac(x, y, z); - } else { - fac = 1.0; - } + fac = GetRefineFac(x, y, z); + } else { + fac = 1.0; + } - int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac); - for (int i_part=0; i_part r; - plasma_injector->getPositionUnitBox(r, i_part, fac); + int ref_num_ppc = num_ppc * AMREX_D_TERM(fac, *fac, *fac); + for (int i_part=0; i_part r; + plasma_injector->getPositionUnitBox(r, i_part, fac); #if ( AMREX_SPACEDIM == 3 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; - Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1]; - Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2]; + Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; + Real y = overlap_corner[1] + (iv[1] + r[1])*dx[1]; + Real z = overlap_corner[2] + (iv[2] + r[2])*dx[2]; #elif ( AMREX_SPACEDIM == 2 ) - Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; - Real y = 0; - Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1]; + Real x = overlap_corner[0] + (iv[0] + r[0])*dx[0]; + Real y = 0; + Real z = overlap_corner[1] + (iv[1] + r[1])*dx[1]; #endif - // If the new particle is not inside the tile box, - // go to the next generated particle. + // If the new particle is not inside the tile box, + // go to the next generated particle. #if ( AMREX_SPACEDIM == 3 ) - if(!tile_realbox.contains( RealVect{x, y, z} )) continue; + if(!tile_realbox.contains( RealVect{x, y, z} )) continue; #elif ( AMREX_SPACEDIM == 2 ) - if(!tile_realbox.contains( RealVect{x, z} )) continue; + if(!tile_realbox.contains( RealVect{x, z} )) continue; #endif - ++np; - } + ++np; + } } return np; @@ -170,8 +170,8 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real // Move the particles to where they will be at t = 0 in the boosted frame if (boost_adjust_transverse_positions) { - x = xpr - tpr*vxpr; - y = ypr - tpr*vypr; + x = xpr - tpr*vxpr; + y = ypr - tpr*vypr; } z = zpr - tpr*vzpr; @@ -323,9 +323,9 @@ void PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) { #ifdef AMREX_USE_GPU - AddPlasmaGPU(lev, part_realbox); + AddPlasmaGPU(lev, part_realbox); #else - AddPlasmaCPU(lev, part_realbox); + AddPlasmaCPU(lev, part_realbox); #endif } @@ -416,7 +416,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) // Count the number of cells in this direction in overlap_realbox overlap_box.setSmall( dir, 0 ); overlap_box.setBig( dir, - int( round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )) - 1); + int( round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )) - 1); } if (no_overlap == 1) { continue; // Go to the next tile @@ -483,54 +483,54 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) Real dens; std::array u; if (WarpX::gamma_boost == 1.){ - // Lab-frame simulation - // If the particle is not within the species's - // xmin, xmax, ymin, ymax, zmin, zmax, go to - // the next generated particle. - if (!plasma_injector->insideBounds(xb, yb, z)) continue; - plasma_injector->getMomentum(u, x, y, z); - dens = plasma_injector->getDensity(x, y, z); + // Lab-frame simulation + // If the particle is not within the species's + // xmin, xmax, ymin, ymax, zmin, zmax, go to + // the next generated particle. + if (!plasma_injector->insideBounds(xb, yb, z)) continue; + plasma_injector->getMomentum(u, x, y, z); + dens = plasma_injector->getDensity(x, y, z); } else { - // Boosted-frame simulation - Real c = PhysConst::c; - Real gamma_boost = WarpX::gamma_boost; - Real beta_boost = WarpX::beta_boost; - // Since the user provides the density distribution - // at t_lab=0 and in the lab-frame coordinates, - // we need to find the lab-frame position of this - // particle at t_lab=0, from its boosted-frame coordinates - // Assuming ballistic motion, this is given by: - // z0_lab = gamma*( z_boost*(1-beta*betaz_lab) - ct_boost*(betaz_lab-beta) ) - // where betaz_lab is the speed of the particle in the lab frame - // - // In order for this equation to be solvable, betaz_lab - // is explicitly assumed to have no dependency on z0_lab - plasma_injector->getMomentum(u, x, y, 0.); // No z0_lab dependency - // At this point u is the lab-frame momentum - // => Apply the above formula for z0_lab - Real gamma_lab = std::sqrt( 1 + (u[0]*u[0] + u[1]*u[1] + u[2]*u[2])/(c*c) ); - Real betaz_lab = u[2]/gamma_lab/c; - Real t = WarpX::GetInstance().gett_new(lev); - Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) ); - // If the particle is not within the lab-frame zmin, zmax, etc. - // go to the next generated particle. - if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue; - // call `getDensity` with lab-frame parameters - dens = plasma_injector->getDensity(x, y, z0_lab); - // At this point u and dens are the lab-frame quantities - // => Perform Lorentz transform - dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab ); - u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab ); + // Boosted-frame simulation + Real c = PhysConst::c; + Real gamma_boost = WarpX::gamma_boost; + Real beta_boost = WarpX::beta_boost; + // Since the user provides the density distribution + // at t_lab=0 and in the lab-frame coordinates, + // we need to find the lab-frame position of this + // particle at t_lab=0, from its boosted-frame coordinates + // Assuming ballistic motion, this is given by: + // z0_lab = gamma*( z_boost*(1-beta*betaz_lab) - ct_boost*(betaz_lab-beta) ) + // where betaz_lab is the speed of the particle in the lab frame + // + // In order for this equation to be solvable, betaz_lab + // is explicitly assumed to have no dependency on z0_lab + plasma_injector->getMomentum(u, x, y, 0.); // No z0_lab dependency + // At this point u is the lab-frame momentum + // => Apply the above formula for z0_lab + Real gamma_lab = std::sqrt( 1 + (u[0]*u[0] + u[1]*u[1] + u[2]*u[2])/(c*c) ); + Real betaz_lab = u[2]/gamma_lab/c; + Real t = WarpX::GetInstance().gett_new(lev); + Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) ); + // If the particle is not within the lab-frame zmin, zmax, etc. + // go to the next generated particle. + if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue; + // call `getDensity` with lab-frame parameters + dens = plasma_injector->getDensity(x, y, z0_lab); + // At this point u and dens are the lab-frame quantities + // => Perform Lorentz transform + dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab ); + u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab ); } Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); #ifdef WARPX_RZ if (plasma_injector->radially_weighted) { - weight *= 2*MathConst::pi*xb; + weight *= 2*MathConst::pi*xb; } else { - // This is not correct since it might shift the particle - // out of the local grid - x = std::sqrt(xb*rmax); - weight *= dx[0]; + // This is not correct since it might shift the particle + // out of the local grid + x = std::sqrt(xb*rmax); + weight *= dx[0]; } #endif attribs[PIdx::w ] = weight; @@ -550,18 +550,18 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) particle_tile.push_back_real(particle_comps["uzold"], u[2]); } - AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); + AddOneParticle(lev, grid_id, tile_id, x, y, z, attribs); } } if (cost) { - wt = (amrex::second() - wt) / tile_box.d_numPts(); + wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); amrex::ParallelFor(tile_box, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - costarr(i,j,k) += wt; - }); + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + costarr(i,j,k) += wt; + }); } } } @@ -655,7 +655,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) // Count the number of cells in this direction in overlap_realbox overlap_box.setSmall( dir, 0 ); overlap_box.setBig( dir, - int( round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )) - 1); + int( round((overlap_realbox.hi(dir)-overlap_realbox.lo(dir))/dx[dir] )) - 1); } if (no_overlap == 1) { continue; // Go to the next tile @@ -664,8 +664,8 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) const int grid_id = mfi.index(); const int tile_id = mfi.LocalTileIndex(); - Cuda::HostVector host_particles; - std::array, PIdx::nattribs> host_attribs; + Cuda::HostVector host_particles; + std::array, PIdx::nattribs> host_attribs; // Loop through the cells of overlap_box and inject // the corresponding particles @@ -725,54 +725,54 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) Real dens; std::array u; if (WarpX::gamma_boost == 1.){ - // Lab-frame simulation - // If the particle is not within the species's - // xmin, xmax, ymin, ymax, zmin, zmax, go to - // the next generated particle. - if (!plasma_injector->insideBounds(xb, yb, z)) continue; - plasma_injector->getMomentum(u, x, y, z); - dens = plasma_injector->getDensity(x, y, z); + // Lab-frame simulation + // If the particle is not within the species's + // xmin, xmax, ymin, ymax, zmin, zmax, go to + // the next generated particle. + if (!plasma_injector->insideBounds(xb, yb, z)) continue; + plasma_injector->getMomentum(u, x, y, z); + dens = plasma_injector->getDensity(x, y, z); } else { - // Boosted-frame simulation - Real c = PhysConst::c; - Real gamma_boost = WarpX::gamma_boost; - Real beta_boost = WarpX::beta_boost; - // Since the user provides the density distribution - // at t_lab=0 and in the lab-frame coordinates, - // we need to find the lab-frame position of this - // particle at t_lab=0, from its boosted-frame coordinates - // Assuming ballistic motion, this is given by: - // z0_lab = gamma*( z_boost*(1-beta*betaz_lab) - ct_boost*(betaz_lab-beta) ) - // where betaz_lab is the speed of the particle in the lab frame - // - // In order for this equation to be solvable, betaz_lab - // is explicitly assumed to have no dependency on z0_lab - plasma_injector->getMomentum(u, x, y, 0.); // No z0_lab dependency - // At this point u is the lab-frame momentum - // => Apply the above formula for z0_lab - Real gamma_lab = std::sqrt( 1 + (u[0]*u[0] + u[1]*u[1] + u[2]*u[2])/(c*c) ); - Real betaz_lab = u[2]/gamma_lab/c; - Real t = WarpX::GetInstance().gett_new(lev); - Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) ); - // If the particle is not within the lab-frame zmin, zmax, etc. - // go to the next generated particle. - if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue; - // call `getDensity` with lab-frame parameters - dens = plasma_injector->getDensity(x, y, z0_lab); - // At this point u and dens are the lab-frame quantities - // => Perform Lorentz transform - dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab ); - u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab ); + // Boosted-frame simulation + Real c = PhysConst::c; + Real gamma_boost = WarpX::gamma_boost; + Real beta_boost = WarpX::beta_boost; + // Since the user provides the density distribution + // at t_lab=0 and in the lab-frame coordinates, + // we need to find the lab-frame position of this + // particle at t_lab=0, from its boosted-frame coordinates + // Assuming ballistic motion, this is given by: + // z0_lab = gamma*( z_boost*(1-beta*betaz_lab) - ct_boost*(betaz_lab-beta) ) + // where betaz_lab is the speed of the particle in the lab frame + // + // In order for this equation to be solvable, betaz_lab + // is explicitly assumed to have no dependency on z0_lab + plasma_injector->getMomentum(u, x, y, 0.); // No z0_lab dependency + // At this point u is the lab-frame momentum + // => Apply the above formula for z0_lab + Real gamma_lab = std::sqrt( 1 + (u[0]*u[0] + u[1]*u[1] + u[2]*u[2])/(c*c) ); + Real betaz_lab = u[2]/gamma_lab/c; + Real t = WarpX::GetInstance().gett_new(lev); + Real z0_lab = gamma_boost * ( z*(1-beta_boost*betaz_lab) - c*t*(betaz_lab-beta_boost) ); + // If the particle is not within the lab-frame zmin, zmax, etc. + // go to the next generated particle. + if (!plasma_injector->insideBounds(xb, yb, z0_lab)) continue; + // call `getDensity` with lab-frame parameters + dens = plasma_injector->getDensity(x, y, z0_lab); + // At this point u and dens are the lab-frame quantities + // => Perform Lorentz transform + dens = gamma_boost * dens * ( 1 - beta_boost*betaz_lab ); + u[2] = gamma_boost * ( u[2] -beta_boost*c*gamma_lab ); } Real weight = dens * scale_fac / (AMREX_D_TERM(fac, *fac, *fac)); #ifdef WARPX_RZ if (plasma_injector->radially_weighted) { - weight *= 2*MathConst::pi*xb; + weight *= 2*MathConst::pi*xb; } else { - // This is not correct since it might shift the particle - // out of the local grid - x = std::sqrt(xb*rmax); - weight *= dx[0]; + // This is not correct since it might shift the particle + // out of the local grid + x = std::sqrt(xb*rmax); + weight *= dx[0]; } #endif attribs[PIdx::w ] = weight; @@ -793,50 +793,50 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) particle_tile.push_back_real(particle_comps["uzold"], u[2]); } - ParticleType p; - p.id() = ParticleType::NextID(); - p.cpu() = ParallelDescriptor::MyProc(); + ParticleType p; + p.id() = ParticleType::NextID(); + p.cpu() = ParallelDescriptor::MyProc(); #if (AMREX_SPACEDIM == 3) - p.pos(0) = x; - p.pos(1) = y; - p.pos(2) = z; + p.pos(0) = x; + p.pos(1) = y; + p.pos(2) = z; #elif (AMREX_SPACEDIM == 2) #ifdef WARPX_RZ attribs[PIdx::theta] = theta; #endif - p.pos(0) = xb; - p.pos(1) = z; + p.pos(0) = xb; + p.pos(1) = z; #endif - host_particles.push_back(p); - for (int kk = 0; kk < PIdx::nattribs; ++kk) - host_attribs[kk].push_back(attribs[kk]); + host_particles.push_back(p); + for (int kk = 0; kk < PIdx::nattribs; ++kk) + host_attribs[kk].push_back(attribs[kk]); } } - auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; + auto& particle_tile = GetParticles(lev)[std::make_pair(grid_id,tile_id)]; auto old_size = particle_tile.GetArrayOfStructs().size(); auto new_size = old_size + host_particles.size(); - particle_tile.resize(new_size); + particle_tile.resize(new_size); Cuda::thrust_copy(host_particles.begin(), host_particles.end(), particle_tile.GetArrayOfStructs().begin() + old_size); - for (int kk = 0; kk < PIdx::nattribs; ++kk) { + for (int kk = 0; kk < PIdx::nattribs; ++kk) { Cuda::thrust_copy(host_attribs[kk].begin(), host_attribs[kk].end(), particle_tile.GetStructOfArrays().GetRealData(kk).begin() + old_size); - } + } if (cost) { - wt = (amrex::second() - wt) / tile_box.d_numPts(); + wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); amrex::ParallelFor(tile_box, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - costarr(i,j,k) += wt; - }); + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + costarr(i,j,k) += wt; + }); } } } @@ -963,13 +963,13 @@ FieldGatherES (const amrex::Vector, WRPX_INTERPOLATE_CIC(particles.dataPtr(), nstride, np, Exp.dataPtr(), Eyp.dataPtr(), #if AMREX_SPACEDIM == 3 - Ezp.dataPtr(), + Ezp.dataPtr(), #endif - exfab.dataPtr(), eyfab.dataPtr(), + exfab.dataPtr(), eyfab.dataPtr(), #if AMREX_SPACEDIM == 3 - ezfab.dataPtr(), + ezfab.dataPtr(), #endif - box.loVect(), box.hiVect(), plo, dx, &ng); + box.loVect(), box.hiVect(), plo, dx, &ng); } else { const FArrayBox& exfab_coarse = coarse_Ex[pti]; @@ -1004,7 +1004,7 @@ FieldGatherES (const amrex::Vector, void PhysicalParticleContainer::EvolveES (const Vector, 3> >& E, - Vector >& rho, + Vector >& rho, Real t, Real dt) { BL_PROFILE("PPC::EvolveES()"); @@ -1014,7 +1014,7 @@ PhysicalParticleContainer::EvolveES (const VectorGeom(lev); const RealBox& prob_domain = gm.ProbDomain(); - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { // Particle structs auto& particles = pti.GetArrayOfStructs(); int nstride = particles.dataShape().first; @@ -1071,11 +1071,11 @@ PhysicalParticleContainer::FieldGather (int lev, { Cuda::ManagedDeviceVector xp, yp, zp; - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) - { + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { Real wt = amrex::second(); - const Box& box = pti.validbox(); + const Box& box = pti.validbox(); auto& attribs = pti.GetAttribs(); @@ -1088,63 +1088,63 @@ PhysicalParticleContainer::FieldGather (int lev, const long np = pti.numParticles(); - // Data on the grid - const FArrayBox& exfab = Ex[pti]; - const FArrayBox& eyfab = Ey[pti]; - const FArrayBox& ezfab = Ez[pti]; - const FArrayBox& bxfab = Bx[pti]; - const FArrayBox& byfab = By[pti]; - const FArrayBox& bzfab = Bz[pti]; - - Exp.assign(np,0.0); - Eyp.assign(np,0.0); - Ezp.assign(np,0.0); - Bxp.assign(np,0.0); - Byp.assign(np,0.0); - Bzp.assign(np,0.0); - - // - // copy data from particle container to temp arrays - // + // Data on the grid + const FArrayBox& exfab = Ex[pti]; + const FArrayBox& eyfab = Ey[pti]; + const FArrayBox& ezfab = Ez[pti]; + const FArrayBox& bxfab = Bx[pti]; + const FArrayBox& byfab = By[pti]; + const FArrayBox& bzfab = Bz[pti]; + + Exp.assign(np,0.0); + Eyp.assign(np,0.0); + Ezp.assign(np,0.0); + Bxp.assign(np,0.0); + Byp.assign(np,0.0); + Bzp.assign(np,0.0); + + // + // copy data from particle container to temp arrays + // pti.GetPosition(xp, yp, zp); const std::array& xyzmin = WarpX::LowerCorner(box, lev); const int* ixyzmin = box.loVect(); - // - // Field Gather - // - const int ll4symtry = false; + // + // Field Gather + // + const int ll4symtry = false; long lvect_fieldgathe = 64; - warpx_geteb_energy_conserving( - &np, - xp.dataPtr(), - yp.dataPtr(), - zp.dataPtr(), - Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(), - Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(), - ixyzmin, - &xyzmin[0], &xyzmin[1], &xyzmin[2], - &dx[0], &dx[1], &dx[2], - &WarpX::nox, &WarpX::noy, &WarpX::noz, - BL_TO_FORTRAN_ANYD(exfab), - BL_TO_FORTRAN_ANYD(eyfab), - BL_TO_FORTRAN_ANYD(ezfab), - BL_TO_FORTRAN_ANYD(bxfab), - BL_TO_FORTRAN_ANYD(byfab), - BL_TO_FORTRAN_ANYD(bzfab), - &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, - &lvect_fieldgathe, &WarpX::field_gathering_algo); + warpx_geteb_energy_conserving( + &np, + xp.dataPtr(), + yp.dataPtr(), + zp.dataPtr(), + Exp.dataPtr(),Eyp.dataPtr(),Ezp.dataPtr(), + Bxp.dataPtr(),Byp.dataPtr(),Bzp.dataPtr(), + ixyzmin, + &xyzmin[0], &xyzmin[1], &xyzmin[2], + &dx[0], &dx[1], &dx[2], + &WarpX::nox, &WarpX::noy, &WarpX::noz, + BL_TO_FORTRAN_ANYD(exfab), + BL_TO_FORTRAN_ANYD(eyfab), + BL_TO_FORTRAN_ANYD(ezfab), + BL_TO_FORTRAN_ANYD(bxfab), + BL_TO_FORTRAN_ANYD(byfab), + BL_TO_FORTRAN_ANYD(bzfab), + &ll4symtry, &WarpX::l_lower_order_in_v, &WarpX::do_nodal, + &lvect_fieldgathe, &WarpX::field_gathering_algo); if (cost) { const Box& tbx = pti.tilebox(); wt = (amrex::second() - wt) / tbx.d_numPts(); Array4 const& costarr = cost->array(pti); amrex::ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - costarr(i,j,k) += wt; - }); + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + costarr(i,j,k) += wt; + }); } } } @@ -1152,9 +1152,9 @@ PhysicalParticleContainer::FieldGather (int lev, void PhysicalParticleContainer::Evolve (int lev, - const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, - const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, - MultiFab& jx, MultiFab& jy, MultiFab& jz, + const MultiFab& Ex, const MultiFab& Ey, const MultiFab& Ez, + const MultiFab& Bx, const MultiFab& By, const MultiFab& Bz, + MultiFab& jx, MultiFab& jy, MultiFab& jz, MultiFab* cjx, MultiFab* cjy, MultiFab* cjz, MultiFab* rho, MultiFab* crho, const MultiFab* cEx, const MultiFab* cEy, const MultiFab* cEz, @@ -1200,8 +1200,8 @@ PhysicalParticleContainer::Evolve (int lev, RealVector tmp; ParticleVector particle_tmp; - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) - { + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) + { Real wt = amrex::second(); const Box& box = pti.validbox(); @@ -1282,14 +1282,14 @@ PhysicalParticleContainer::Evolve (int lev, #endif } - Exp.assign(np,0.0); - Eyp.assign(np,0.0); - Ezp.assign(np,0.0); - Bxp.assign(np,WarpX::B_external[0]); - Byp.assign(np,WarpX::B_external[1]); - Bzp.assign(np,WarpX::B_external[2]); + Exp.assign(np,0.0); + Eyp.assign(np,0.0); + Ezp.assign(np,0.0); + Bxp.assign(np,WarpX::B_external[0]); + Byp.assign(np,WarpX::B_external[1]); + Bzp.assign(np,WarpX::B_external[2]); - m_giv[thread_num].resize(np); + m_giv[thread_num].resize(np); long nfine_current = np; long nfine_gather = np; @@ -1384,12 +1384,12 @@ PhysicalParticleContainer::Evolve (int lev, const long np_current = (cjx) ? nfine_current : np; - // - // copy data from particle container to temp arrays - // - BL_PROFILE_VAR_START(blp_copy); + // + // copy data from particle container to temp arrays + // + BL_PROFILE_VAR_START(blp_copy); pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); - BL_PROFILE_VAR_STOP(blp_copy); + BL_PROFILE_VAR_STOP(blp_copy); if (rho) DepositCharge(pti, wp, rho, crho, 0, np_current, np, thread_num, lev); @@ -1568,10 +1568,10 @@ PhysicalParticleContainer::Evolve (int lev, wt = (amrex::second() - wt) / tbx.d_numPts(); Array4 const& costarr = cost->array(pti); amrex::ParallelFor(tbx, - [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept - { - costarr(i,j,k) += wt; - }); + [=] AMREX_GPU_DEVICE (int i, int j, int k) noexcept + { + costarr(i,j,k) += wt; + }); } } } @@ -1603,9 +1603,9 @@ PhysicalParticleContainer::SplitParticles(int lev) { pti.GetPosition(xp, yp, zp); const std::array& dx = WarpX::CellSize(lev); - // particle Array Of Structs data + // particle Array Of Structs data auto& particles = pti.GetArrayOfStructs(); - // particle Struct Of Arrays data + // particle Struct Of Arrays data auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w ]; auto& uxp = attribs[PIdx::ux]; @@ -1615,13 +1615,13 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; i& xp, + Cuda::ManagedDeviceVector& xp, Cuda::ManagedDeviceVector& yp, Cuda::ManagedDeviceVector& zp, Cuda::ManagedDeviceVector& giv, @@ -1795,8 +1795,8 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = 0; #endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) - { - const Box& box = pti.validbox(); + { + const Box& box = pti.validbox(); auto& attribs = pti.GetAttribs(); @@ -1812,26 +1812,26 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const long np = pti.numParticles(); - // Data on the grid - const FArrayBox& exfab = Ex[pti]; - const FArrayBox& eyfab = Ey[pti]; - const FArrayBox& ezfab = Ez[pti]; - const FArrayBox& bxfab = Bx[pti]; - const FArrayBox& byfab = By[pti]; - const FArrayBox& bzfab = Bz[pti]; - - Exp.assign(np,0.0); - Eyp.assign(np,0.0); - Ezp.assign(np,0.0); - Bxp.assign(np,WarpX::B_external[0]); - Byp.assign(np,WarpX::B_external[1]); - Bzp.assign(np,WarpX::B_external[2]); - - m_giv[thread_num].resize(np); - - // - // copy data from particle container to temp arrays - // + // Data on the grid + const FArrayBox& exfab = Ex[pti]; + const FArrayBox& eyfab = Ey[pti]; + const FArrayBox& ezfab = Ez[pti]; + const FArrayBox& bxfab = Bx[pti]; + const FArrayBox& byfab = By[pti]; + const FArrayBox& bzfab = Bz[pti]; + + Exp.assign(np,0.0); + Eyp.assign(np,0.0); + Ezp.assign(np,0.0); + Bxp.assign(np,WarpX::B_external[0]); + Byp.assign(np,WarpX::B_external[1]); + Bzp.assign(np,WarpX::B_external[2]); + + m_giv[thread_num].resize(np); + + // + // copy data from particle container to temp arrays + // pti.GetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); @@ -1875,7 +1875,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, } void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, - const Real* yp, const Real* zp) + const Real* yp, const Real* zp) { auto& attribs = pti.GetAttribs(); @@ -1894,16 +1894,16 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, const long np = pti.numParticles(); ParallelFor( np, - [=] AMREX_GPU_DEVICE (long i) { - xpold[i]=xp[i]; - ypold[i]=yp[i]; - zpold[i]=zp[i]; + [=] AMREX_GPU_DEVICE (long i) { + xpold[i]=xp[i]; + ypold[i]=yp[i]; + zpold[i]=zp[i]; - uxpold[i]=uxp[i]; - uypold[i]=uyp[i]; - uzpold[i]=uzp[i]; - } - ); + uxpold[i]=uxp[i]; + uypold[i]=uyp[i]; + uzpold[i]=uzp[i]; + } + ); } void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real z_old, -- cgit v1.2.3