From 747fdc9b046a00a3d48996ca6a1c86404edabe6c Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 5 Sep 2019 18:44:21 +0200 Subject: Initial work to add a lambda component --- Source/Particles/PhysicalParticleContainer.cpp | 141 ++++++++++++++----------- 1 file changed, 80 insertions(+), 61 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 015482e3f..5cefe7419 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -51,7 +51,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); if (not do_user_plot_vars){ // By default, all particle variables are dumped to plotfiles, - // including {x,y,z,ux,uy,uz}old variables when running in a + // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); @@ -68,9 +68,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ - // Return error if var not in PIdx. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ParticleStringNames::to_index.count(var), + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } @@ -89,6 +89,11 @@ void PhysicalParticleContainer::InitData() // Init ionization module here instead of in the PhysicalParticleContainer // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} + +#ifdef WARPX_QED + if(do_qed){InitLambda();} +#endif + AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute } @@ -142,7 +147,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -154,7 +159,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution distz(z_m, z_rms); if (ParallelDescriptor::IOProcessor()) { - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -386,11 +391,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int tile_id = mfi.LocalTileIndex(); // Max number of new particles, if particles are created in the whole - // overlap_box. All of them are created, and invalid ones are then + // overlap_box. All of them are created, and invalid ones are then // discaded int max_new_particles = overlap_box.numPts() * num_ppc; - // If refine injection, build pointer dp_cellid that holds pointer to + // If refine injection, build pointer dp_cellid that holds pointer to // array of refined cell IDs. Vector cellid_v; if (refine_injection and lev == 0) @@ -445,7 +450,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } - + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -457,9 +462,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; - // Loop over all new particles and inject them (creates too many + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). - // The invalid ones are given negative ID and are deleted during the + // The invalid ones are given negative ID and are deleted during the // next redistribute. amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept { @@ -614,7 +619,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = z; #endif }, shared_mem_bytes); - + if (cost) { wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); @@ -852,7 +857,7 @@ PhysicalParticleContainer::FieldGather (int lev, MultiFab* cost = WarpX::getCosts(lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -902,7 +907,7 @@ PhysicalParticleContainer::FieldGather (int lev, // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); if (cost) { @@ -935,11 +940,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // 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; @@ -964,9 +969,9 @@ PhysicalParticleContainer::Evolve (int lev, tmp_particle_data[lev][index][i].resize(np); } } - + #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1165,7 +1170,7 @@ PhysicalParticleContainer::Evolve (int lev, } const long np_current = (cjx) ? nfine_current : np; - + // // copy data from particle container to temp arrays // @@ -1181,14 +1186,14 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (has_buffer){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } } - + if (! do_not_push) { const long np_gather = (cEx) ? nfine_gather : np; @@ -1200,7 +1205,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - exfab, eyfab, ezfab, bxfab, byfab, bzfab, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); if (np_gather < np) @@ -1215,7 +1220,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbxfab = &(*cBx)[pti]; FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; - + if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1253,13 +1258,13 @@ PhysicalParticleContainer::Evolve (int lev, 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)); 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)); bzeli = filtered_Bz.elixir(); @@ -1267,14 +1272,14 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + // Field gather for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, - cEx->nGrow(), e_is_nodal, - nfine_gather, np-nfine_gather, + cEx->nGrow(), e_is_nodal, + nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1284,7 +1289,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_ppc_pp); @@ -1298,7 +1303,7 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - + // Deposit inside domains DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -1318,7 +1323,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) { // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1387,7 +1392,7 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; icharge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; @@ -1594,7 +1599,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1634,7 +1639,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); // This wraps the momentum advance so that inheritors can modify the call. @@ -1681,7 +1686,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - + const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); @@ -1697,7 +1702,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, xpold[i]=xp[i]; ypold[i]=yp[i]; zpold[i]=zp[i]; - + uxpold[i]=uxp[i]; uypold[i]=uyp[i]; uzpold[i]=uzp[i]; @@ -1737,7 +1742,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1784,7 +1789,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -1819,7 +1824,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); @@ -1844,7 +1849,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) AddPlasma(lev, injection_box); } -/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, +/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp. * \param Exp-Bzp: fields on particles. * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti @@ -1854,7 +1859,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param np_to_gather: number of particles onto which fields are gathered * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located - * \param gather_lev: level from which particles gather fields (lev-1) for + * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. */ void @@ -1881,14 +1886,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - + // If no particles, do not do anything if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); // Set staggering shift depending on e_is_nodal const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; - + // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. Box box; @@ -1898,26 +1903,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(gather_lev); box = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Add guard cells to the box. box.grow(ngE); - + const Array4& ex_arr = exfab->array(); const Array4& ey_arr = eyfab->array(); const Array4& ez_arr = ezfab->array(); const Array4& bx_arr = bxfab->array(); const Array4& by_arr = byfab->array(); const Array4& bz_arr = bzfab->array(); - + const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; - + // Lower corner of tile box physical domain const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev); - + const Dim3 lo = lbound(box); - + // Depending on l_lower_in_v and WarpX::nox, call // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ @@ -2000,10 +2005,10 @@ void PhysicalParticleContainer::InitIonizationModule () } // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) // For now, we assume l=0 and m=0. - // The approximate expressions are used, + // The approximate expressions are used, // without Gamma function Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; - Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * std::pow(PhysConst::alpha,4)/PhysConst::r_e; Real UH = table_ionization_energies[0]; Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; @@ -2018,18 +2023,18 @@ void PhysicalParticleContainer::InitIonizationModule () Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; - adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * + * * \param mfi: tile or grid * \param lev: MR level * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 + * This function initialized the array, and set each element to 1 or 0 * depending on whether the particle is ionized or not. */ void @@ -2075,7 +2080,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Loop over all particles in grid/tile. If ionized, set mask to 1 // and increment ionization level. - ParallelFor( + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level @@ -2107,3 +2112,17 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const } ); } + + +#ifdef WARPX_QED + +// initializes an additional real component if QED effects are enabled +void +PhysicalParticleContainer::InitLambda() +{ + if(!do_qed) return; + + // Add runtime real component for ionization level + AddRealComp("lambda"); +} +#endif -- cgit v1.2.3 From a4ee8fc2c33b36ec60784e2b312818295cf1860e Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 6 Sep 2019 17:35:33 +0200 Subject: Add lamba if qed effects are enabled --- Examples/Physics_applications/QED/inputs.2d | 3 +++ Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 7 +++++-- 3 files changed, 9 insertions(+), 3 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Physics_applications/QED/inputs.2d b/Examples/Physics_applications/QED/inputs.2d index e2f3aec34..2bf6f1f4d 100644 --- a/Examples/Physics_applications/QED/inputs.2d +++ b/Examples/Physics_applications/QED/inputs.2d @@ -82,6 +82,9 @@ photons.uz_m = 500. photons.ux_th = 2. photons.uy_th = 2. photons.uz_th = 50. +##########QED#################### +photons.do_qed = 1 +################################# ################################# diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index af4b8b2db..a9a2610dc 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -133,7 +133,7 @@ public: #ifdef WARPX_QED // a flag to enable/disable QED effects - bool do_qed=true; //QED_TO_CHANGE + bool do_qed=false; // initializes an additional real component if QED effects are enabled void InitLambda(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ebf923732..51567f026 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -46,6 +46,10 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp "version of Redistribute in AMReX does not work with runtime parameters"); #endif +#ifdef WARPX_QED + pp.query("do_qed", do_qed); +#endif + pp.query("plot_species", plot_species); int do_user_plot_vars; do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); @@ -2130,14 +2134,13 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const #ifdef WARPX_QED - // initializes an additional real component if QED effects are enabled void PhysicalParticleContainer::InitLambda() { if(!do_qed) return; - // Add runtime real component for ionization level AddRealComp("lambda"); + } #endif -- cgit v1.2.3 From 6d3b11d4e2349979916907afbdbf903ba9eae186 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 6 Sep 2019 19:07:10 +0200 Subject: fixed name --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 51567f026..eb12ad97d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2140,7 +2140,7 @@ PhysicalParticleContainer::InitLambda() { if(!do_qed) return; // Add runtime real component for ionization level - AddRealComp("lambda"); + AddRealComp("tau"); } #endif -- cgit v1.2.3 From 8b5122ced93b9ec90f9df12ce21cceb5f22c33bf Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 6 Sep 2019 19:30:57 +0200 Subject: Added lambda (has a bug) --- Source/Diagnostics/ParticleIO.cpp | 26 +++++---- Source/Particles/PhysicalParticleContainer.H | 4 +- Source/Particles/PhysicalParticleContainer.cpp | 1 + Source/Particles/WarpXParticleContainer.H | 76 ++++++++++++++------------ 4 files changed, 58 insertions(+), 49 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 8bfa45a59..3785d41d2 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -22,8 +22,8 @@ RigidInjectedParticleContainer::ReadHeader (std::istream& is) int zinject_plane_tmp; is >> zinject_plane_tmp; zinject_plane_levels.push_back(zinject_plane_tmp); - WarpX::GotoNextLine(is); - } + WarpX::GotoNextLine(is); + } for (int i = 0; i < nlevs; ++i) { @@ -31,7 +31,7 @@ RigidInjectedParticleContainer::ReadHeader (std::istream& is) is >> done_injecting_tmp; done_injecting.push_back(done_injecting_tmp); WarpX::GotoNextLine(is); - } + } } void @@ -78,7 +78,7 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const { for (unsigned i = 0, n = species_names.size(); i < n; ++i) { - auto& pc = allcontainers[i]; + auto& pc = allcontainers[i]; if (pc->plot_species) { Vector real_names; @@ -90,19 +90,19 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); real_names.push_back("momentum_z"); - + real_names.push_back("Ex"); real_names.push_back("Ey"); real_names.push_back("Ez"); - + real_names.push_back("Bx"); real_names.push_back("By"); real_names.push_back("Bz"); - + #ifdef WARPX_DIM_RZ real_names.push_back("theta"); #endif - + if(pc->do_field_ionization){ int_names.push_back("ionization_level"); // int_flags specifies, for each integer attribs, whether it is @@ -112,6 +112,12 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const int_flags.resize(1, 1); } +#ifdef WARPX_QED + if(pc->do_qed){ + real_names.push_back("tau"); + } +#endif + // Convert momentum to SI pc->ConvertUnits(ConvertDirection::WarpX_to_SI); // real_names contains a list of all particle attributes. @@ -134,7 +140,7 @@ MultiParticleContainer::Restart (const std::string& dir) } void -MultiParticleContainer::ReadHeader (std::istream& is) +MultiParticleContainer::ReadHeader (std::istream& is) { for (auto& pc : allcontainers) { pc->ReadHeader(is); @@ -149,7 +155,7 @@ MultiParticleContainer::WriteHeader (std::ostream& os) const } } -// Particle momentum is defined as gamma*velocity, which is neither +// Particle momentum is defined as gamma*velocity, which is neither // SI mass*gamma*velocity nor normalized gamma*velocity/c. // This converts momentum to SI units (or vice-versa) to write SI data // to file. diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index a9a2610dc..45b92f0da 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -132,10 +132,8 @@ public: virtual void ConvertUnits (ConvertDirection convert_dir) override; #ifdef WARPX_QED - // a flag to enable/disable QED effects - bool do_qed=false; - // initializes an additional real component if QED effects are enabled + // a flag to enable/disable QED effects (do_qed) is defined in WarpXParticleContainer void InitLambda(); #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index eb12ad97d..b43facb98 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2141,6 +2141,7 @@ PhysicalParticleContainer::InitLambda() if(!do_qed) return; // Add runtime real component for ionization level AddRealComp("tau"); + plot_flags.resize(PIdx::nattribs + 1, 1); } #endif diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index c39eec9dc..cb295c516 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -73,19 +73,19 @@ public: const amrex::Cuda::ManagedDeviceVector& y, const amrex::Cuda::ManagedDeviceVector& z); #endif - const std::array& GetAttribs () const { - return GetStructOfArrays().GetRealData(); + const std::array& GetAttribs () const { + return GetStructOfArrays().GetRealData(); } - - std::array& GetAttribs () { - return GetStructOfArrays().GetRealData(); + + std::array& GetAttribs () { + return GetStructOfArrays().GetRealData(); } - const RealVector& GetAttribs (int comp) const { + const RealVector& GetAttribs (int comp) const { return GetStructOfArrays().GetRealData(comp); } - - RealVector& GetAttribs (int comp) { + + RealVector& GetAttribs (int comp) { return GetStructOfArrays().GetRealData(comp); } @@ -102,15 +102,15 @@ class WarpXParticleContainer public: friend MultiParticleContainer; - // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components // and 0 int components for the particle data. using DiagnosticParticleData = amrex::StructOfArrays; // DiagnosticParticles is a vector, with one element per MR level. - // DiagnosticParticles[lev] is typically a key-value pair where the key is - // a pair [grid_index, tile_index], and the value is the corresponding + // DiagnosticParticles[lev] is typically a key-value pair where the key is + // a pair [grid_index, tile_index], and the value is the corresponding // DiagnosticParticleData (see above) on this tile. using DiagnosticParticles = amrex::Vector, DiagnosticParticleData> >; - + WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); virtual ~WarpXParticleContainer() {} @@ -124,12 +124,12 @@ public: const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) {} -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void EvolveES (const amrex::Vector, 3> >& E, - amrex::Vector >& rho, + amrex::Vector >& rho, amrex::Real t, amrex::Real dt) = 0; #endif // WARPX_DO_ELECTROSTATIC - + virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, @@ -143,10 +143,10 @@ public: virtual void PostRestart () = 0; virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) {} - + void AllocData (); /// @@ -154,7 +154,7 @@ public: /// It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. /// This is the electrostatic version of the particle push. - /// + /// void PushXES (amrex::Real dt); /// @@ -162,13 +162,13 @@ public: /// It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. /// This is the electromagnetic version of the particle push. - /// + /// void PushX ( amrex::Real dt); void PushX (int lev, amrex::Real dt); /// /// This pushes the particle momenta by dt. - /// + /// virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -208,9 +208,9 @@ public: 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. + // 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. // Current status: // PhysicalParticleContainer: implemented. // LaserParticleContainer: implemented. @@ -219,7 +219,7 @@ public: // Update optional sub-class-specific injection location. virtual void UpdateContinuousInjectionPosition(amrex::Real dt) {} - /// + /// /// This returns the total charge for all the particles in this ParticleContainer. /// This is needed when solving Poisson's equation with periodic boundary conditions. /// @@ -259,9 +259,9 @@ public: // split along axes (0) or diagonals (1) int split_type = 0; - using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddRealComp; + using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddRealComp; using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddIntComp; - + void AddRealComp (const std::string& name, bool comm=true) { particle_comps[name] = NumRealComps(); @@ -274,7 +274,7 @@ public: AddIntComp(comm); } - int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector& ionization_mask) @@ -286,7 +286,7 @@ protected: std::map particle_comps; std::map particle_icomps; - + int species_id; amrex::Real charge; @@ -296,9 +296,9 @@ protected: static int do_not_push; - // Whether to allow particles outside of the simulation domain to be + // Whether to allow particles outside of the simulation domain to be // initialized when they enter the domain. - // This is currently required because continuous injection does not + // This is currently required because continuous injection does not // support all features allowed by direct injection. int do_continuous_injection = 0; @@ -312,9 +312,13 @@ protected: amrex::Gpu::ManagedVector adk_prefactor; amrex::Gpu::ManagedVector adk_exp_prefactor; std::string physical_element; - + int do_boosted_frame_diags = 1; +#ifdef WARPX_QED + bool do_qed = false; +#endif + amrex::Vector local_rho; amrex::Vector local_jx; amrex::Vector local_jy; @@ -322,20 +326,20 @@ protected: using DataContainer = amrex::Gpu::ManagedDeviceVector; using PairIndex = std::pair; - + amrex::Vector m_xp, m_yp, m_zp, m_giv; - // Whether to dump particle quantities. + // Whether to dump particle quantities. // If true, particle position is always dumped. int plot_species = 1; - // For all particle attribs (execept position), whether or not + // For all particle attribs (execept position), whether or not // to dump to file. amrex::Vector plot_flags; // list of names of attributes to dump. amrex::Vector plot_vars; - + amrex::Vector > > tmp_particle_data; - + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; -- cgit v1.2.3 From 20d7f6f4d06665895e7e935859a3b53e3c0466b5 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 11 Sep 2019 16:42:39 +0200 Subject: Revert style change --- Source/Diagnostics/ParticleIO.cpp | 20 ++-- Source/Particles/PhysicalParticleContainer.H | 10 +- Source/Particles/PhysicalParticleContainer.cpp | 122 ++++++++++++------------- Source/Particles/WarpXParticleContainer.H | 72 +++++++-------- 4 files changed, 112 insertions(+), 112 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 3785d41d2..936febe42 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -22,8 +22,8 @@ RigidInjectedParticleContainer::ReadHeader (std::istream& is) int zinject_plane_tmp; is >> zinject_plane_tmp; zinject_plane_levels.push_back(zinject_plane_tmp); - WarpX::GotoNextLine(is); - } + WarpX::GotoNextLine(is); + } for (int i = 0; i < nlevs; ++i) { @@ -31,7 +31,7 @@ RigidInjectedParticleContainer::ReadHeader (std::istream& is) is >> done_injecting_tmp; done_injecting.push_back(done_injecting_tmp); WarpX::GotoNextLine(is); - } + } } void @@ -78,7 +78,7 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const { for (unsigned i = 0, n = species_names.size(); i < n; ++i) { - auto& pc = allcontainers[i]; + auto& pc = allcontainers[i]; if (pc->plot_species) { Vector real_names; @@ -90,19 +90,19 @@ MultiParticleContainer::WritePlotFile (const std::string& dir) const real_names.push_back("momentum_x"); real_names.push_back("momentum_y"); real_names.push_back("momentum_z"); - + real_names.push_back("Ex"); real_names.push_back("Ey"); real_names.push_back("Ez"); - + real_names.push_back("Bx"); real_names.push_back("By"); real_names.push_back("Bz"); - + #ifdef WARPX_DIM_RZ real_names.push_back("theta"); #endif - + if(pc->do_field_ionization){ int_names.push_back("ionization_level"); // int_flags specifies, for each integer attribs, whether it is @@ -140,7 +140,7 @@ MultiParticleContainer::Restart (const std::string& dir) } void -MultiParticleContainer::ReadHeader (std::istream& is) +MultiParticleContainer::ReadHeader (std::istream& is) { for (auto& pc : allcontainers) { pc->ReadHeader(is); @@ -155,7 +155,7 @@ MultiParticleContainer::WriteHeader (std::ostream& os) const } } -// Particle momentum is defined as gamma*velocity, which is neither +// Particle momentum is defined as gamma*velocity, which is neither // SI mass*gamma*velocity nor normalized gamma*velocity/c. // This converts momentum to SI units (or vice-versa) to write SI data // to file. diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 45b92f0da..b7b33f654 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -24,7 +24,7 @@ public: void InitIonizationModule (); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector, 3> >& E, const amrex::Vector > > >& masks) override; @@ -32,7 +32,7 @@ public: amrex::Vector >& rho, amrex::Real t, amrex::Real dt) override; #endif // WARPX_DO_ELECTROSTATIC - + virtual void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -99,14 +99,14 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} void SplitParticles(int lev); - + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector& ionization_mask) override; @@ -125,7 +125,7 @@ public: std::array u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b43facb98..46223d347 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -55,7 +55,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); if (not do_user_plot_vars){ // By default, all particle variables are dumped to plotfiles, - // including {x,y,z,ux,uy,uz}old variables when running in a + // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); @@ -72,9 +72,9 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ - // Return error if var not in PIdx. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ParticleStringNames::to_index.count(var), + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } @@ -151,7 +151,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -163,7 +163,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution distz(z_m, z_rms); if (ParallelDescriptor::IOProcessor()) { - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -396,11 +396,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int tile_id = mfi.LocalTileIndex(); // Max number of new particles, if particles are created in the whole - // overlap_box. All of them are created, and invalid ones are then + // overlap_box. All of them are created, and invalid ones are then // discaded int max_new_particles = overlap_box.numPts() * num_ppc; - // If refine injection, build pointer dp_cellid that holds pointer to + // If refine injection, build pointer dp_cellid that holds pointer to // array of refined cell IDs. Vector cellid_v; if (refine_injection and lev == 0) @@ -455,7 +455,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } - + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -467,9 +467,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; - // Loop over all new particles and inject them (creates too many + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). - // The invalid ones are given negative ID and are deleted during the + // The invalid ones are given negative ID and are deleted during the // next redistribute. amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept { @@ -636,7 +636,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = z; #endif }, shared_mem_bytes); - + if (cost) { wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); @@ -874,7 +874,7 @@ PhysicalParticleContainer::FieldGather (int lev, MultiFab* cost = WarpX::getCosts(lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -924,7 +924,7 @@ PhysicalParticleContainer::FieldGather (int lev, // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); if (cost) { @@ -957,11 +957,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // 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; @@ -986,9 +986,9 @@ PhysicalParticleContainer::Evolve (int lev, tmp_particle_data[lev][index][i].resize(np); } } - + #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1187,7 +1187,7 @@ PhysicalParticleContainer::Evolve (int lev, } const long np_current = (cjx) ? nfine_current : np; - + // // copy data from particle container to temp arrays // @@ -1203,14 +1203,14 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (has_buffer){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } } - + if (! do_not_push) { const long np_gather = (cEx) ? nfine_gather : np; @@ -1222,7 +1222,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - exfab, eyfab, ezfab, bxfab, byfab, bzfab, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); @@ -1238,7 +1238,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbxfab = &(*cBx)[pti]; FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; - + if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1276,13 +1276,13 @@ PhysicalParticleContainer::Evolve (int lev, 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)); 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)); bzeli = filtered_Bz.elixir(); @@ -1290,14 +1290,14 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + // Field gather for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, - cEx->nGrow(), e_is_nodal, - nfine_gather, np-nfine_gather, + cEx->nGrow(), e_is_nodal, + nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1307,7 +1307,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_ppc_pp); @@ -1321,7 +1321,7 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - + // Deposit inside domains DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -1341,7 +1341,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) { // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1410,7 +1410,7 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; icharge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; @@ -1617,7 +1617,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1657,7 +1657,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -1705,7 +1705,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - + const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); @@ -1721,7 +1721,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, xpold[i]=xp[i]; ypold[i]=yp[i]; zpold[i]=zp[i]; - + uxpold[i]=uxp[i]; uypold[i]=uyp[i]; uzpold[i]=uzp[i]; @@ -1761,7 +1761,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1808,7 +1808,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -1843,7 +1843,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); @@ -1868,7 +1868,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) AddPlasma(lev, injection_box); } -/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, +/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp. * \param Exp-Bzp: fields on particles. * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti @@ -1878,7 +1878,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param np_to_gather: number of particles onto which fields are gathered * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located - * \param gather_lev: level from which particles gather fields (lev-1) for + * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. */ void @@ -1905,14 +1905,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - + // If no particles, do not do anything if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); // Set staggering shift depending on e_is_nodal const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; - + // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. Box box; @@ -1922,26 +1922,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(gather_lev); box = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Add guard cells to the box. box.grow(ngE); - + const Array4& ex_arr = exfab->array(); const Array4& ey_arr = eyfab->array(); const Array4& ez_arr = ezfab->array(); const Array4& bx_arr = bxfab->array(); const Array4& by_arr = byfab->array(); const Array4& bz_arr = bzfab->array(); - + const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; - + // Lower corner of tile box physical domain const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev); - + const Dim3 lo = lbound(box); - + // Depending on l_lower_in_v and WarpX::nox, call // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ @@ -2024,10 +2024,10 @@ void PhysicalParticleContainer::InitIonizationModule () } // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) // For now, we assume l=0 and m=0. - // The approximate expressions are used, + // The approximate expressions are used, // without Gamma function Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; - Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * std::pow(PhysConst::alpha,4)/PhysConst::r_e; Real UH = table_ionization_energies[0]; Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; @@ -2042,18 +2042,18 @@ void PhysicalParticleContainer::InitIonizationModule () Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; - adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * + * * \param mfi: tile or grid * \param lev: MR level * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 + * This function initialized the array, and set each element to 1 or 0 * depending on whether the particle is ionized or not. */ void @@ -2099,7 +2099,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Loop over all particles in grid/tile. If ionized, set mask to 1 // and increment ionization level. - ParallelFor( + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index cb295c516..b30c81cb1 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -73,19 +73,19 @@ public: const amrex::Cuda::ManagedDeviceVector& y, const amrex::Cuda::ManagedDeviceVector& z); #endif - const std::array& GetAttribs () const { - return GetStructOfArrays().GetRealData(); + const std::array& GetAttribs () const { + return GetStructOfArrays().GetRealData(); } - - std::array& GetAttribs () { - return GetStructOfArrays().GetRealData(); + + std::array& GetAttribs () { + return GetStructOfArrays().GetRealData(); } - const RealVector& GetAttribs (int comp) const { + const RealVector& GetAttribs (int comp) const { return GetStructOfArrays().GetRealData(comp); } - - RealVector& GetAttribs (int comp) { + + RealVector& GetAttribs (int comp) { return GetStructOfArrays().GetRealData(comp); } @@ -102,15 +102,15 @@ class WarpXParticleContainer public: friend MultiParticleContainer; - // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components // and 0 int components for the particle data. using DiagnosticParticleData = amrex::StructOfArrays; // DiagnosticParticles is a vector, with one element per MR level. - // DiagnosticParticles[lev] is typically a key-value pair where the key is - // a pair [grid_index, tile_index], and the value is the corresponding + // DiagnosticParticles[lev] is typically a key-value pair where the key is + // a pair [grid_index, tile_index], and the value is the corresponding // DiagnosticParticleData (see above) on this tile. using DiagnosticParticles = amrex::Vector, DiagnosticParticleData> >; - + WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); virtual ~WarpXParticleContainer() {} @@ -124,12 +124,12 @@ public: const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) {} -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void EvolveES (const amrex::Vector, 3> >& E, - amrex::Vector >& rho, + amrex::Vector >& rho, amrex::Real t, amrex::Real dt) = 0; #endif // WARPX_DO_ELECTROSTATIC - + virtual void Evolve (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz, @@ -143,10 +143,10 @@ public: virtual void PostRestart () = 0; virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) {} - + void AllocData (); /// @@ -154,7 +154,7 @@ public: /// It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. /// This is the electrostatic version of the particle push. - /// + /// void PushXES (amrex::Real dt); /// @@ -162,13 +162,13 @@ public: /// It is used to desynchronize the particles after initializaton /// or when restarting from a checkpoint. /// This is the electromagnetic version of the particle push. - /// + /// void PushX ( amrex::Real dt); void PushX (int lev, amrex::Real dt); /// /// This pushes the particle momenta by dt. - /// + /// virtual void PushP (int lev, amrex::Real dt, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -208,9 +208,9 @@ public: 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. + // 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. // Current status: // PhysicalParticleContainer: implemented. // LaserParticleContainer: implemented. @@ -219,7 +219,7 @@ public: // Update optional sub-class-specific injection location. virtual void UpdateContinuousInjectionPosition(amrex::Real dt) {} - /// + /// /// This returns the total charge for all the particles in this ParticleContainer. /// This is needed when solving Poisson's equation with periodic boundary conditions. /// @@ -259,9 +259,9 @@ public: // split along axes (0) or diagonals (1) int split_type = 0; - using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddRealComp; + using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddRealComp; using amrex::ParticleContainer<0, 0, PIdx::nattribs>::AddIntComp; - + void AddRealComp (const std::string& name, bool comm=true) { particle_comps[name] = NumRealComps(); @@ -274,7 +274,7 @@ public: AddIntComp(comm); } - int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } + int DoBoostedFrameDiags () const { return do_boosted_frame_diags; } virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector& ionization_mask) @@ -286,7 +286,7 @@ protected: std::map particle_comps; std::map particle_icomps; - + int species_id; amrex::Real charge; @@ -296,9 +296,9 @@ protected: static int do_not_push; - // Whether to allow particles outside of the simulation domain to be + // Whether to allow particles outside of the simulation domain to be // initialized when they enter the domain. - // This is currently required because continuous injection does not + // This is currently required because continuous injection does not // support all features allowed by direct injection. int do_continuous_injection = 0; @@ -312,7 +312,7 @@ protected: amrex::Gpu::ManagedVector adk_prefactor; amrex::Gpu::ManagedVector adk_exp_prefactor; std::string physical_element; - + int do_boosted_frame_diags = 1; #ifdef WARPX_QED @@ -326,20 +326,20 @@ protected: using DataContainer = amrex::Gpu::ManagedDeviceVector; using PairIndex = std::pair; - + amrex::Vector m_xp, m_yp, m_zp, m_giv; - // Whether to dump particle quantities. + // Whether to dump particle quantities. // If true, particle position is always dumped. int plot_species = 1; - // For all particle attribs (execept position), whether or not + // For all particle attribs (execept position), whether or not // to dump to file. amrex::Vector plot_flags; // list of names of attributes to dump. amrex::Vector plot_vars; - + amrex::Vector > > tmp_particle_data; - + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; -- cgit v1.2.3 From 8943396e2775221899a628ed876066f776651e22 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 12 Sep 2019 12:13:51 +0200 Subject: Fixed issue with tau output --- Source/Particles/MultiParticleContainer.cpp | 1 - Source/Particles/PhysicalParticleContainer.H | 2 +- Source/Particles/PhysicalParticleContainer.cpp | 45 ++++++++++++++++---------- 3 files changed, 29 insertions(+), 19 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index e8125588d..74e59867e 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -17,7 +17,6 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) { - ReadParameters(); allcontainers.resize(nspecies + nlasers); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index b7b33f654..3826423ad 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -132,7 +132,7 @@ public: virtual void ConvertUnits (ConvertDirection convert_dir) override; #ifdef WARPX_QED - // initializes an additional real component if QED effects are enabled + // initializes an additional real component (the optical depth) if QED effects are enabled // a flag to enable/disable QED effects (do_qed) is defined in WarpXParticleContainer void InitLambda(); #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 46223d347..4abfb8fbd 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -46,35 +46,46 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp "version of Redistribute in AMReX does not work with runtime parameters"); #endif + #ifdef WARPX_QED + //Add real component if QED is enabled pp.query("do_qed", do_qed); + if(do_qed){ + AddRealComp("tau"); + } #endif + //variable to set plot_flags size + int plot_flag_size = PIdx::nattribs; + if(WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + plot_flag_size += 6; + +#ifdef WARPX_QED + if(do_qed){ + // plot_flag will have an entry for the optical depth + plot_flag_size ++; + } +#endif + //_______________________________ + pp.query("plot_species", plot_species); int do_user_plot_vars; do_user_plot_vars = pp.queryarr("plot_vars", plot_vars); if (not do_user_plot_vars){ // By default, all particle variables are dumped to plotfiles, - // including {x,y,z,ux,uy,uz}old variables when running in a + // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ - plot_flags.resize(PIdx::nattribs + 6, 1); - } else { - plot_flags.resize(PIdx::nattribs, 1); - } + plot_flags.resize(plot_flag_size, 1); } else { // Set plot_flag to 0 for all attribs - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ - plot_flags.resize(PIdx::nattribs + 6, 0); - } else { - plot_flags.resize(PIdx::nattribs, 0); - } + plot_flags.resize(plot_flag_size, 0); + // If not none, set plot_flags values to 1 for elements in plot_vars. if (plot_vars[0] != "none"){ for (const auto& var : plot_vars){ - // Return error if var not in PIdx. - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - ParticleStringNames::to_index.count(var), + // Return error if var not in PIdx. + AMREX_ALWAYS_ASSERT_WITH_MESSAGE( + ParticleStringNames::to_index.count(var), "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } @@ -927,6 +938,7 @@ PhysicalParticleContainer::FieldGather (int lev, &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); + if (cost) { const Box& tbx = pti.tilebox(); wt = (amrex::second() - wt) / tbx.d_numPts(); @@ -2139,9 +2151,8 @@ void PhysicalParticleContainer::InitLambda() { if(!do_qed) return; - // Add runtime real component for ionization level - AddRealComp("tau"); - plot_flags.resize(PIdx::nattribs + 1, 1); + // Add runtime real component for optical depth + //TODO } #endif -- cgit v1.2.3 From d16ab16c18524ecd776ec8d35e2c48f4b0148f1f Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 12 Sep 2019 12:19:45 +0200 Subject: removed InitTau from PhysicalParticleContainer --- Source/Particles/PhysicalParticleContainer.H | 16 +-- Source/Particles/PhysicalParticleContainer.cpp | 133 +++++++++++-------------- 2 files changed, 63 insertions(+), 86 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 3826423ad..7946b4650 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -24,7 +24,7 @@ public: void InitIonizationModule (); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector, 3> >& E, const amrex::Vector > > >& masks) override; @@ -32,7 +32,7 @@ public: amrex::Vector >& rho, amrex::Real t, amrex::Real dt) override; #endif // WARPX_DO_ELECTROSTATIC - + virtual void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -99,14 +99,14 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} void SplitParticles(int lev); - + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector& ionization_mask) override; @@ -125,18 +125,12 @@ public: std::array u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; virtual void ConvertUnits (ConvertDirection convert_dir) override; -#ifdef WARPX_QED - // initializes an additional real component (the optical depth) if QED effects are enabled - // a flag to enable/disable QED effects (do_qed) is defined in WarpXParticleContainer - void InitLambda(); -#endif - protected: std::string species_name; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4abfb8fbd..fe69ef5dc 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -75,7 +75,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp // By default, all particle variables are dumped to plotfiles, // including {x,y,z,ux,uy,uz}old variables when running in a // boosted frame - plot_flags.resize(plot_flag_size, 1); + plot_flags.resize(plot_flag_size, 1); } else { // Set plot_flag to 0 for all attribs plot_flags.resize(plot_flag_size, 0); @@ -105,10 +105,6 @@ void PhysicalParticleContainer::InitData() // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} -#ifdef WARPX_QED - if(do_qed){InitLambda();} -#endif - AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute } @@ -162,7 +158,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -174,7 +170,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution distz(z_m, z_rms); if (ParallelDescriptor::IOProcessor()) { - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -407,11 +403,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int tile_id = mfi.LocalTileIndex(); // Max number of new particles, if particles are created in the whole - // overlap_box. All of them are created, and invalid ones are then + // overlap_box. All of them are created, and invalid ones are then // discaded int max_new_particles = overlap_box.numPts() * num_ppc; - // If refine injection, build pointer dp_cellid that holds pointer to + // If refine injection, build pointer dp_cellid that holds pointer to // array of refined cell IDs. Vector cellid_v; if (refine_injection and lev == 0) @@ -466,7 +462,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } - + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -478,9 +474,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; - // Loop over all new particles and inject them (creates too many + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). - // The invalid ones are given negative ID and are deleted during the + // The invalid ones are given negative ID and are deleted during the // next redistribute. amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept { @@ -647,7 +643,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = z; #endif }, shared_mem_bytes); - + if (cost) { wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); @@ -885,7 +881,7 @@ PhysicalParticleContainer::FieldGather (int lev, MultiFab* cost = WarpX::getCosts(lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -935,7 +931,7 @@ PhysicalParticleContainer::FieldGather (int lev, // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -969,11 +965,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // 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; @@ -998,9 +994,9 @@ PhysicalParticleContainer::Evolve (int lev, tmp_particle_data[lev][index][i].resize(np); } } - + #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1199,7 +1195,7 @@ PhysicalParticleContainer::Evolve (int lev, } const long np_current = (cjx) ? nfine_current : np; - + // // copy data from particle container to temp arrays // @@ -1215,14 +1211,14 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (has_buffer){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } } - + if (! do_not_push) { const long np_gather = (cEx) ? nfine_gather : np; @@ -1234,7 +1230,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - exfab, eyfab, ezfab, bxfab, byfab, bzfab, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); @@ -1250,7 +1246,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbxfab = &(*cBx)[pti]; FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; - + if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1288,13 +1284,13 @@ PhysicalParticleContainer::Evolve (int lev, 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)); 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)); bzeli = filtered_Bz.elixir(); @@ -1302,14 +1298,14 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + // Field gather for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, - cEx->nGrow(), e_is_nodal, - nfine_gather, np-nfine_gather, + cEx->nGrow(), e_is_nodal, + nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1319,7 +1315,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_ppc_pp); @@ -1333,7 +1329,7 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - + // Deposit inside domains DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -1353,7 +1349,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) { // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1422,7 +1418,7 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; icharge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; @@ -1629,7 +1625,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1669,7 +1665,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -1717,7 +1713,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - + const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); @@ -1733,7 +1729,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, xpold[i]=xp[i]; ypold[i]=yp[i]; zpold[i]=zp[i]; - + uxpold[i]=uxp[i]; uypold[i]=uyp[i]; uzpold[i]=uzp[i]; @@ -1773,7 +1769,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1820,7 +1816,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -1855,7 +1851,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); @@ -1880,7 +1876,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) AddPlasma(lev, injection_box); } -/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, +/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp. * \param Exp-Bzp: fields on particles. * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti @@ -1890,7 +1886,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param np_to_gather: number of particles onto which fields are gathered * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located - * \param gather_lev: level from which particles gather fields (lev-1) for + * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. */ void @@ -1917,14 +1913,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - + // If no particles, do not do anything if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); // Set staggering shift depending on e_is_nodal const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; - + // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. Box box; @@ -1934,26 +1930,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(gather_lev); box = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Add guard cells to the box. box.grow(ngE); - + const Array4& ex_arr = exfab->array(); const Array4& ey_arr = eyfab->array(); const Array4& ez_arr = ezfab->array(); const Array4& bx_arr = bxfab->array(); const Array4& by_arr = byfab->array(); const Array4& bz_arr = bzfab->array(); - + const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; - + // Lower corner of tile box physical domain const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev); - + const Dim3 lo = lbound(box); - + // Depending on l_lower_in_v and WarpX::nox, call // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ @@ -2036,10 +2032,10 @@ void PhysicalParticleContainer::InitIonizationModule () } // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) // For now, we assume l=0 and m=0. - // The approximate expressions are used, + // The approximate expressions are used, // without Gamma function Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; - Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * std::pow(PhysConst::alpha,4)/PhysConst::r_e; Real UH = table_ionization_energies[0]; Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; @@ -2054,18 +2050,18 @@ void PhysicalParticleContainer::InitIonizationModule () Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; - adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * + * * \param mfi: tile or grid * \param lev: MR level * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 + * This function initialized the array, and set each element to 1 or 0 * depending on whether the particle is ionized or not. */ void @@ -2111,7 +2107,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Loop over all particles in grid/tile. If ionized, set mask to 1 // and increment ionization level. - ParallelFor( + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level @@ -2143,16 +2139,3 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const } ); } - - -#ifdef WARPX_QED -// initializes an additional real component if QED effects are enabled -void -PhysicalParticleContainer::InitLambda() -{ - if(!do_qed) return; - // Add runtime real component for optical depth - //TODO - -} -#endif -- cgit v1.2.3 From a7057e0309c023a3f1a3c64ff4339c8cd3da9be4 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 12 Sep 2019 12:23:18 +0200 Subject: reverted style changes --- Source/Particles/PhysicalParticleContainer.H | 10 +-- Source/Particles/PhysicalParticleContainer.cpp | 115 ++++++++++++------------- 2 files changed, 62 insertions(+), 63 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 7946b4650..c7494fbdf 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -24,7 +24,7 @@ public: void InitIonizationModule (); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector, 3> >& E, const amrex::Vector > > >& masks) override; @@ -32,7 +32,7 @@ public: amrex::Vector >& rho, amrex::Real t, amrex::Real dt) override; #endif // WARPX_DO_ELECTROSTATIC - + virtual void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -99,14 +99,14 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} void SplitParticles(int lev); - + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector& ionization_mask) override; @@ -125,7 +125,7 @@ public: std::array u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index fe69ef5dc..da4e3afc2 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -104,7 +104,6 @@ void PhysicalParticleContainer::InitData() // Init ionization module here instead of in the PhysicalParticleContainer // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} - AddParticles(0); // Note - add on level 0 Redistribute(); // We then redistribute } @@ -158,7 +157,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -170,7 +169,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution distz(z_m, z_rms); if (ParallelDescriptor::IOProcessor()) { - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -403,11 +402,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int tile_id = mfi.LocalTileIndex(); // Max number of new particles, if particles are created in the whole - // overlap_box. All of them are created, and invalid ones are then + // overlap_box. All of them are created, and invalid ones are then // discaded int max_new_particles = overlap_box.numPts() * num_ppc; - // If refine injection, build pointer dp_cellid that holds pointer to + // If refine injection, build pointer dp_cellid that holds pointer to // array of refined cell IDs. Vector cellid_v; if (refine_injection and lev == 0) @@ -462,7 +461,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } - + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -474,9 +473,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; - // Loop over all new particles and inject them (creates too many + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). - // The invalid ones are given negative ID and are deleted during the + // The invalid ones are given negative ID and are deleted during the // next redistribute. amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept { @@ -643,7 +642,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = z; #endif }, shared_mem_bytes); - + if (cost) { wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); @@ -881,7 +880,7 @@ PhysicalParticleContainer::FieldGather (int lev, MultiFab* cost = WarpX::getCosts(lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -931,7 +930,7 @@ PhysicalParticleContainer::FieldGather (int lev, // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -965,11 +964,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // 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; @@ -994,9 +993,9 @@ PhysicalParticleContainer::Evolve (int lev, tmp_particle_data[lev][index][i].resize(np); } } - + #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1195,7 +1194,7 @@ PhysicalParticleContainer::Evolve (int lev, } const long np_current = (cjx) ? nfine_current : np; - + // // copy data from particle container to temp arrays // @@ -1211,14 +1210,14 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (has_buffer){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } } - + if (! do_not_push) { const long np_gather = (cEx) ? nfine_gather : np; @@ -1230,7 +1229,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - exfab, eyfab, ezfab, bxfab, byfab, bzfab, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); @@ -1246,7 +1245,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbxfab = &(*cBx)[pti]; FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; - + if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1284,13 +1283,13 @@ PhysicalParticleContainer::Evolve (int lev, 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)); 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)); bzeli = filtered_Bz.elixir(); @@ -1298,14 +1297,14 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + // Field gather for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, - cEx->nGrow(), e_is_nodal, - nfine_gather, np-nfine_gather, + cEx->nGrow(), e_is_nodal, + nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1315,7 +1314,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_ppc_pp); @@ -1329,7 +1328,7 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - + // Deposit inside domains DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -1349,7 +1348,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) { // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1418,7 +1417,7 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; icharge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; @@ -1625,7 +1624,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1665,7 +1664,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -1713,7 +1712,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - + const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); @@ -1729,7 +1728,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, xpold[i]=xp[i]; ypold[i]=yp[i]; zpold[i]=zp[i]; - + uxpold[i]=uxp[i]; uypold[i]=uyp[i]; uzpold[i]=uzp[i]; @@ -1769,7 +1768,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1816,7 +1815,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -1851,7 +1850,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); @@ -1876,7 +1875,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) AddPlasma(lev, injection_box); } -/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, +/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp. * \param Exp-Bzp: fields on particles. * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti @@ -1886,7 +1885,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param np_to_gather: number of particles onto which fields are gathered * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located - * \param gather_lev: level from which particles gather fields (lev-1) for + * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. */ void @@ -1913,14 +1912,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - + // If no particles, do not do anything if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); // Set staggering shift depending on e_is_nodal const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; - + // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. Box box; @@ -1930,26 +1929,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(gather_lev); box = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Add guard cells to the box. box.grow(ngE); - + const Array4& ex_arr = exfab->array(); const Array4& ey_arr = eyfab->array(); const Array4& ez_arr = ezfab->array(); const Array4& bx_arr = bxfab->array(); const Array4& by_arr = byfab->array(); const Array4& bz_arr = bzfab->array(); - + const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; - + // Lower corner of tile box physical domain const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev); - + const Dim3 lo = lbound(box); - + // Depending on l_lower_in_v and WarpX::nox, call // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ @@ -2032,10 +2031,10 @@ void PhysicalParticleContainer::InitIonizationModule () } // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) // For now, we assume l=0 and m=0. - // The approximate expressions are used, + // The approximate expressions are used, // without Gamma function Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; - Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * std::pow(PhysConst::alpha,4)/PhysConst::r_e; Real UH = table_ionization_energies[0]; Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; @@ -2050,18 +2049,18 @@ void PhysicalParticleContainer::InitIonizationModule () Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; - adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * + * * \param mfi: tile or grid * \param lev: MR level * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 + * This function initialized the array, and set each element to 1 or 0 * depending on whether the particle is ionized or not. */ void @@ -2107,7 +2106,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Loop over all particles in grid/tile. If ionized, set mask to 1 // and increment ionization level. - ParallelFor( + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level -- cgit v1.2.3 From 85ca7fb2a3786956af222af095f500c9a05a4b9d Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 12 Sep 2019 15:34:44 +0200 Subject: Now physical particles have Tau initialization via QS engine --- Examples/Physics_applications/QED/inputs.2d | 4 + Source/Particles/MultiParticleContainer.H | 11 +- Source/Particles/PhotonParticleContainer.cpp | 22 ++-- Source/Particles/PhysicalParticleContainer.H | 19 +++- Source/Particles/PhysicalParticleContainer.cpp | 152 +++++++++++++++---------- Source/QED/Make.package | 2 + Source/QED/amrex_rng_wrapper.h | 4 +- Source/QED/breit_wheeler_engine_wrapper.cpp | 16 +++ 8 files changed, 147 insertions(+), 83 deletions(-) create mode 100644 Source/QED/breit_wheeler_engine_wrapper.cpp (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Physics_applications/QED/inputs.2d b/Examples/Physics_applications/QED/inputs.2d index f08670736..cb7b41703 100644 --- a/Examples/Physics_applications/QED/inputs.2d +++ b/Examples/Physics_applications/QED/inputs.2d @@ -61,6 +61,10 @@ electrons.uz_m = 500. electrons.ux_th = 2. electrons.uy_th = 2. electrons.uz_th = 50. +##########QED#################### +electrons.do_qed = 1 +electrons.do_qed_quantum_sync = 1 +################################# photons.charge = -q_e photons.mass = m_e diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 67392c449..e377351bd 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -24,6 +24,7 @@ // #ifdef WARPX_QED #include "breit_wheeler_engine_wrapper.h" + #include "quantum_sync_engine_wrapper.h" #endif class MultiParticleContainer @@ -92,7 +93,7 @@ public: void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, - const amrex::MultiFab& By, const amrex::MultiFab& Bz); + const amrex::MultiFab& By, const amrex::MultiFab& Bz); /// /// This evolves all the particles by one PIC time step, including current deposition, the @@ -129,7 +130,7 @@ public: /// /// This deposits the particle charge onto a node-centered MultiFab and returns a unique ptr /// to it. The charge density is accumulated over all the particles in the MultiParticleContainer - /// + /// std::unique_ptr GetChargeDensity(int lev, bool local = false); void doFieldIonization (); @@ -137,7 +138,7 @@ public: void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; - + void Restart (const std::string& dir); void PostRestart (); @@ -191,9 +192,9 @@ public: PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } - //bw_engine is a public member of the MultiParticleContainer object #ifdef WARPX_QED warpx_breit_wheeler_engine bw_engine; + warpx_quantum_sync_engine qs_engine; #endif protected: @@ -220,7 +221,7 @@ private: void mapSpeciesProduct (); int getSpeciesID (std::string product_str); - + // Number of species dumped in BoostedFrameDiagnostics int nspecies_boosted_frame_diags = 0; // map_species_boosted_frame_diags[i] is the species ID in diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 5fecaf525..1da8cace2 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -129,17 +129,17 @@ PhotonParticleContainer::Evolve (int lev, void PhotonParticleContainer::InitTauBreitWheeler() { BL_PROFILE("PhotonParticleContainer::InitOpticalDepth"); -//Looping over all the particles -int num_levels = finestLevel() + 1; -for (int lev=0; lev < num_levels; ++lev) - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti){ - auto taus = pti.GetAttribs(particle_comps["tau"]).dataPtr(); - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - taus[i] = warpx_breit_wheeler_engine::get_optical_depth(); - } - ); + //Looping over all the particles + int num_levels = finestLevel() + 1; + for (int lev=0; lev < num_levels; ++lev) + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti){ + auto taus = pti.GetAttribs(particle_comps["tau"]).dataPtr(); + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + taus[i] = warpx_breit_wheeler_engine::get_optical_depth(); + } + ); } } #endif diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index c7494fbdf..578668fdb 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -24,7 +24,7 @@ public: void InitIonizationModule (); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector, 3> >& E, const amrex::Vector > > >& masks) override; @@ -32,7 +32,7 @@ public: amrex::Vector >& rho, amrex::Real t, amrex::Real dt) override; #endif // WARPX_DO_ELECTROSTATIC - + virtual void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -99,14 +99,14 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} void SplitParticles(int lev); - + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector& ionization_mask) override; @@ -125,7 +125,7 @@ public: std::array u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; @@ -145,6 +145,15 @@ protected: // Inject particles during the whole simulation void ContinuousInjection (const amrex::RealBox& injection_box) override; +private: +#ifdef WARPX_QED + // A flag to enable quantum_synchrotron process for leptons + bool do_qed_quantum_sync = false; + + // A function to initialize the Tau component according to the QS engine + void InitTauQuantumSync(); +#endif + }; #endif diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index da4e3afc2..ae226015d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -50,9 +50,14 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp #ifdef WARPX_QED //Add real component if QED is enabled pp.query("do_qed", do_qed); - if(do_qed){ + if(do_qed) AddRealComp("tau"); - } + + //IF do_qed is enabled, find out if Quantum Synchrotron process is enabled + if(do_qed) + pp.query("do_qed_quantum_sync", do_qed_quantum_sync); + + //TODO: SHOULD CHECK IF SPECIES IS EITHER ELECTRONS OR POSITRONS!! #endif //variable to set plot_flags size @@ -63,7 +68,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp #ifdef WARPX_QED if(do_qed){ // plot_flag will have an entry for the optical depth - plot_flag_size ++; + plot_flag_size++; } #endif //_______________________________ @@ -104,7 +109,14 @@ void PhysicalParticleContainer::InitData() // Init ionization module here instead of in the PhysicalParticleContainer // constructor because dt is required if (do_field_ionization) {InitIonizationModule();} + AddParticles(0); // Note - add on level 0 + +#ifdef WARPX_QED + if(do_qed_quantum_sync) + InitTauQuantumSync(); +#endif + Redistribute(); // We then redistribute } @@ -157,7 +169,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -169,7 +181,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution distz(z_m, z_rms); if (ParallelDescriptor::IOProcessor()) { - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -402,11 +414,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int tile_id = mfi.LocalTileIndex(); // Max number of new particles, if particles are created in the whole - // overlap_box. All of them are created, and invalid ones are then + // overlap_box. All of them are created, and invalid ones are then // discaded int max_new_particles = overlap_box.numPts() * num_ppc; - // If refine injection, build pointer dp_cellid that holds pointer to + // If refine injection, build pointer dp_cellid that holds pointer to // array of refined cell IDs. Vector cellid_v; if (refine_injection and lev == 0) @@ -461,7 +473,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } - + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -473,9 +485,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; - // Loop over all new particles and inject them (creates too many + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). - // The invalid ones are given negative ID and are deleted during the + // The invalid ones are given negative ID and are deleted during the // next redistribute. amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept { @@ -642,7 +654,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = z; #endif }, shared_mem_bytes); - + if (cost) { wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); @@ -880,7 +892,7 @@ PhysicalParticleContainer::FieldGather (int lev, MultiFab* cost = WarpX::getCosts(lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -930,7 +942,7 @@ PhysicalParticleContainer::FieldGather (int lev, // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -964,11 +976,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // 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; @@ -993,9 +1005,9 @@ PhysicalParticleContainer::Evolve (int lev, tmp_particle_data[lev][index][i].resize(np); } } - + #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1194,7 +1206,7 @@ PhysicalParticleContainer::Evolve (int lev, } const long np_current = (cjx) ? nfine_current : np; - + // // copy data from particle container to temp arrays // @@ -1210,14 +1222,14 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (has_buffer){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } } - + if (! do_not_push) { const long np_gather = (cEx) ? nfine_gather : np; @@ -1229,7 +1241,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - exfab, eyfab, ezfab, bxfab, byfab, bzfab, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); @@ -1245,7 +1257,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbxfab = &(*cBx)[pti]; FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; - + if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1283,13 +1295,13 @@ PhysicalParticleContainer::Evolve (int lev, 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)); 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)); bzeli = filtered_Bz.elixir(); @@ -1297,14 +1309,14 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + // Field gather for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, - cEx->nGrow(), e_is_nodal, - nfine_gather, np-nfine_gather, + cEx->nGrow(), e_is_nodal, + nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1314,7 +1326,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_ppc_pp); @@ -1328,7 +1340,7 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - + // Deposit inside domains DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -1348,7 +1360,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) { // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1417,7 +1429,7 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; icharge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; @@ -1624,7 +1636,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1664,7 +1676,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -1712,7 +1724,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - + const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); @@ -1728,7 +1740,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, xpold[i]=xp[i]; ypold[i]=yp[i]; zpold[i]=zp[i]; - + uxpold[i]=uxp[i]; uypold[i]=uyp[i]; uzpold[i]=uzp[i]; @@ -1768,7 +1780,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1815,7 +1827,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -1850,7 +1862,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); @@ -1875,7 +1887,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) AddPlasma(lev, injection_box); } -/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, +/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp. * \param Exp-Bzp: fields on particles. * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti @@ -1885,7 +1897,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param np_to_gather: number of particles onto which fields are gathered * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located - * \param gather_lev: level from which particles gather fields (lev-1) for + * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. */ void @@ -1912,14 +1924,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - + // If no particles, do not do anything if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); // Set staggering shift depending on e_is_nodal const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; - + // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. Box box; @@ -1929,26 +1941,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(gather_lev); box = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Add guard cells to the box. box.grow(ngE); - + const Array4& ex_arr = exfab->array(); const Array4& ey_arr = eyfab->array(); const Array4& ez_arr = ezfab->array(); const Array4& bx_arr = bxfab->array(); const Array4& by_arr = byfab->array(); const Array4& bz_arr = bzfab->array(); - + const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; - + // Lower corner of tile box physical domain const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev); - + const Dim3 lo = lbound(box); - + // Depending on l_lower_in_v and WarpX::nox, call // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ @@ -2031,10 +2043,10 @@ void PhysicalParticleContainer::InitIonizationModule () } // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) // For now, we assume l=0 and m=0. - // The approximate expressions are used, + // The approximate expressions are used, // without Gamma function Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; - Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * std::pow(PhysConst::alpha,4)/PhysConst::r_e; Real UH = table_ionization_energies[0]; Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; @@ -2049,18 +2061,18 @@ void PhysicalParticleContainer::InitIonizationModule () Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; - adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * + * * \param mfi: tile or grid * \param lev: MR level * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 + * This function initialized the array, and set each element to 1 or 0 * depending on whether the particle is ionized or not. */ void @@ -2106,7 +2118,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Loop over all particles in grid/tile. If ionized, set mask to 1 // and increment ionization level. - ParallelFor( + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level @@ -2138,3 +2150,23 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const } ); } + +#ifdef WARPX_QED +// A function to initialize the Tau component according to the QS engine +void PhysicalParticleContainer::InitTauQuantumSync() +{ + BL_PROFILE("PhysicalParticleContainer::InitTauQuantumSync"); +//Looping over all the particles + int num_levels = finestLevel() + 1; + for (int lev=0; lev < num_levels; ++lev) + for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti){ + auto taus = pti.GetAttribs(particle_comps["tau"]).dataPtr(); + amrex::ParallelFor( + pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + taus[i] = warpx_quantum_sync_engine::get_optical_depth(); + } + ); + } +} +#endif diff --git a/Source/QED/Make.package b/Source/QED/Make.package index 73ac89228..5c71d44a2 100644 --- a/Source/QED/Make.package +++ b/Source/QED/Make.package @@ -1,8 +1,10 @@ CEXE_headers += qed_wrapper_commons.h CEXE_headers += amrex_rng_wrapper.h CEXE_headers += breit_wheeler_engine_wrapper.h +CEXE_headers += quantum_syncs_engine_wrapper.h CEXE_sources += amrex_rng_wrapper.cpp CEXE_sources += breit_wheeler_engine_wrapper.cpp +CEXE_sources += quantum_sync_engine_wrapper.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/QED VPATH_LOCATIONS += $(WARPX_HOME)/Source/QED diff --git a/Source/QED/amrex_rng_wrapper.h b/Source/QED/amrex_rng_wrapper.h index f9818c589..1e6b7c259 100644 --- a/Source/QED/amrex_rng_wrapper.h +++ b/Source/QED/amrex_rng_wrapper.h @@ -1,5 +1,5 @@ #ifndef WARPX_amrex_rng_wrapper_h_ -#define WARPX_rng_wrapper_h_ +#define WARPX_amrex_rng_wrapper_h_ //This file provides a wrapper aroud the RNG //provided by the amrex library @@ -17,4 +17,4 @@ public: amrex::Real exp(amrex::Real l); }; -#endif //amrex_rng_wrapper_hpp_ +#endif //WARPX_amrex_rng_wrapper_h_ diff --git a/Source/QED/breit_wheeler_engine_wrapper.cpp b/Source/QED/breit_wheeler_engine_wrapper.cpp new file mode 100644 index 000000000..acef9026d --- /dev/null +++ b/Source/QED/breit_wheeler_engine_wrapper.cpp @@ -0,0 +1,16 @@ +#include "breit_wheeler_engine_wrapper.h" +//This file provides a wrapper aroud the breit_wheeler engine +//provided by the standard template library + +using namespace picsar::multi_physics; + +warpx_breit_wheeler_engine::warpx_breit_wheeler_engine(): + breit_wheeler_engine(std::move(amrex_rng_wrapper{})) +{}; + +//Interface for the get_optical_depth method of the BW engine +amrex::Real +AMREX_GPU_HOST_DEVICE +warpx_breit_wheeler_engine::get_optical_depth(){ + return internal_get_optical_depth(amrex::Random()); +} -- cgit v1.2.3 From bd1e4b6801095e7813bcf418c28bb784ee730075 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Thu, 12 Sep 2019 15:44:49 +0200 Subject: Revert style change --- Source/Particles/MultiParticleContainer.H | 8 +- Source/Particles/PhysicalParticleContainer.H | 10 +-- Source/Particles/PhysicalParticleContainer.cpp | 114 ++++++++++++------------- 3 files changed, 66 insertions(+), 66 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index e377351bd..acf1be4c3 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -93,7 +93,7 @@ public: void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, const amrex::MultiFab& Ez, const amrex::MultiFab& Bx, - const amrex::MultiFab& By, const amrex::MultiFab& Bz); + const amrex::MultiFab& By, const amrex::MultiFab& Bz); /// /// This evolves all the particles by one PIC time step, including current deposition, the @@ -130,7 +130,7 @@ public: /// /// This deposits the particle charge onto a node-centered MultiFab and returns a unique ptr /// to it. The charge density is accumulated over all the particles in the MultiParticleContainer - /// + /// std::unique_ptr GetChargeDensity(int lev, bool local = false); void doFieldIonization (); @@ -138,7 +138,7 @@ public: void Checkpoint (const std::string& dir) const; void WritePlotFile (const std::string& dir) const; - + void Restart (const std::string& dir); void PostRestart (); @@ -221,7 +221,7 @@ private: void mapSpeciesProduct (); int getSpeciesID (std::string product_str); - + // Number of species dumped in BoostedFrameDiagnostics int nspecies_boosted_frame_diags = 0; // map_species_boosted_frame_diags[i] is the species ID in diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 578668fdb..200a724f2 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -24,7 +24,7 @@ public: void InitIonizationModule (); -#ifdef WARPX_DO_ELECTROSTATIC +#ifdef WARPX_DO_ELECTROSTATIC virtual void FieldGatherES(const amrex::Vector, 3> >& E, const amrex::Vector > > >& masks) override; @@ -32,7 +32,7 @@ public: amrex::Vector >& rho, amrex::Real t, amrex::Real dt) override; #endif // WARPX_DO_ELECTROSTATIC - + virtual void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -99,14 +99,14 @@ public: const amrex::MultiFab& Bx, const amrex::MultiFab& By, const amrex::MultiFab& Bz) override; - + void copy_attribs(WarpXParIter& pti,const amrex::Real* xp, const amrex::Real* yp, const amrex::Real* zp); virtual void PostRestart () final {} void SplitParticles(int lev); - + virtual void buildIonizationMask (const amrex::MFIter& mfi, const int lev, amrex::Gpu::ManagedDeviceVector& ionization_mask) override; @@ -125,7 +125,7 @@ public: std::array u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, - const amrex::Real z_new, const amrex::Real t_boost, + const amrex::Real z_new, const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ae226015d..12af2714d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -34,7 +34,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); pp.query("do_continuous_injection", do_continuous_injection); - // Whether to plot back-transformed (lab-frame) diagnostics + // Whether to plot back-transformed (lab-frame) diagnostics // for this species. pp.query("do_boosted_frame_diags", do_boosted_frame_diags); @@ -169,7 +169,7 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart, + Real q_tot, long npart, int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); @@ -181,7 +181,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution distz(z_m, z_rms); if (ParallelDescriptor::IOProcessor()) { - // If do_symmetrize, create 4x fewer particles, and + // If do_symmetrize, create 4x fewer particles, and // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) if (do_symmetrize){ npart /= 4; @@ -414,11 +414,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) const int tile_id = mfi.LocalTileIndex(); // Max number of new particles, if particles are created in the whole - // overlap_box. All of them are created, and invalid ones are then + // overlap_box. All of them are created, and invalid ones are then // discaded int max_new_particles = overlap_box.numPts() * num_ppc; - // If refine injection, build pointer dp_cellid that holds pointer to + // If refine injection, build pointer dp_cellid that holds pointer to // array of refined cell IDs. Vector cellid_v; if (refine_injection and lev == 0) @@ -473,7 +473,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if (do_field_ionization) { pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } - + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -485,9 +485,9 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) bool loc_do_field_ionization = do_field_ionization; int loc_ionization_initial_level = ionization_initial_level; - // Loop over all new particles and inject them (creates too many + // Loop over all new particles and inject them (creates too many // particles, in particular does not consider xmin, xmax etc.). - // The invalid ones are given negative ID and are deleted during the + // The invalid ones are given negative ID and are deleted during the // next redistribute. amrex::For(max_new_particles, [=] AMREX_GPU_DEVICE (int ip) noexcept { @@ -654,7 +654,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) p.pos(1) = z; #endif }, shared_mem_bytes); - + if (cost) { wt = (amrex::second() - wt) / tile_box.d_numPts(); Array4 const& costarr = cost->array(mfi); @@ -892,7 +892,7 @@ PhysicalParticleContainer::FieldGather (int lev, MultiFab* cost = WarpX::getCosts(lev); #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -942,7 +942,7 @@ PhysicalParticleContainer::FieldGather (int lev, // int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -976,11 +976,11 @@ PhysicalParticleContainer::Evolve (int lev, BL_PROFILE_VAR_NS("PPC::FieldGather", blp_fg); BL_PROFILE_VAR_NS("PPC::ParticlePush", blp_ppc_pp); BL_PROFILE_VAR_NS("PPC::Evolve::partition", blp_partition); - + const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - // Get instances of NCI Godfrey filters + // 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; @@ -1005,9 +1005,9 @@ PhysicalParticleContainer::Evolve (int lev, tmp_particle_data[lev][index][i].resize(np); } } - + #ifdef _OPENMP -#pragma omp parallel +#pragma omp parallel #endif { #ifdef _OPENMP @@ -1206,7 +1206,7 @@ PhysicalParticleContainer::Evolve (int lev, } const long np_current = (cjx) ? nfine_current : np; - + // // copy data from particle container to temp arrays // @@ -1222,14 +1222,14 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - DepositCharge(pti, wp, ion_lev, rho, 0, 0, + DepositCharge(pti, wp, ion_lev, rho, 0, 0, np_current, thread_num, lev, lev); if (has_buffer){ DepositCharge(pti, wp, ion_lev, crho, 0, np_current, np-np_current, thread_num, lev, lev-1); } } - + if (! do_not_push) { const long np_gather = (cEx) ? nfine_gather : np; @@ -1241,7 +1241,7 @@ PhysicalParticleContainer::Evolve (int lev, // BL_PROFILE_VAR_START(blp_fg); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - exfab, eyfab, ezfab, bxfab, byfab, bzfab, + exfab, eyfab, ezfab, bxfab, byfab, bzfab, Ex.nGrow(), e_is_nodal, 0, np_gather, thread_num, lev, lev); @@ -1257,7 +1257,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* cbxfab = &(*cBx)[pti]; FArrayBox const* cbyfab = &(*cBy)[pti]; FArrayBox const* cbzfab = &(*cBz)[pti]; - + if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1295,13 +1295,13 @@ PhysicalParticleContainer::Evolve (int lev, 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)); 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)); bzeli = filtered_Bz.elixir(); @@ -1309,14 +1309,14 @@ PhysicalParticleContainer::Evolve (int lev, cbzfab = &filtered_Bz; #endif } - + // Field gather for particles in gather buffers e_is_nodal = cEx->is_nodal() and cEy->is_nodal() and cEz->is_nodal(); - FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, cexfab, ceyfab, cezfab, cbxfab, cbyfab, cbzfab, - cEx->nGrow(), e_is_nodal, - nfine_gather, np-nfine_gather, + cEx->nGrow(), e_is_nodal, + nfine_gather, np-nfine_gather, thread_num, lev, lev-1); } @@ -1326,7 +1326,7 @@ PhysicalParticleContainer::Evolve (int lev, // Particle Push // BL_PROFILE_VAR_START(blp_ppc_pp); - PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], + PushPX(pti, m_xp[thread_num], m_yp[thread_num], m_zp[thread_num], m_giv[thread_num], dt); BL_PROFILE_VAR_STOP(blp_ppc_pp); @@ -1340,7 +1340,7 @@ PhysicalParticleContainer::Evolve (int lev, } else { ion_lev = nullptr; } - + // Deposit inside domains DepositCurrent(pti, wp, uxp, uyp, uzp, ion_lev, &jx, &jy, &jz, 0, np_current, thread_num, @@ -1360,7 +1360,7 @@ PhysicalParticleContainer::Evolve (int lev, pti.SetPosition(m_xp[thread_num], m_yp[thread_num], m_zp[thread_num]); BL_PROFILE_VAR_STOP(blp_copy); } - + if (rho) { // Deposit charge after particle push, in component 1 of MultiFab rho. int* AMREX_RESTRICT ion_lev; @@ -1429,7 +1429,7 @@ PhysicalParticleContainer::SplitParticles(int lev) for(int i=0; icharge; const Real m = this-> mass; if (WarpX::particle_pusher_algo == ParticlePusherAlgo::Boris){ - amrex::ParallelFor( + amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { Real qp = q; @@ -1636,7 +1636,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int thread_num = omp_get_thread_num(); #else int thread_num = 0; -#endif +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { const Box& box = pti.validbox(); @@ -1676,7 +1676,7 @@ PhysicalParticleContainer::PushP (int lev, Real dt, int e_is_nodal = Ex.is_nodal() and Ey.is_nodal() and Ez.is_nodal(); FieldGather(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, - &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, + &exfab, &eyfab, &ezfab, &bxfab, &byfab, &bzfab, Ex.nGrow(), e_is_nodal, 0, np, thread_num, lev, lev); @@ -1724,7 +1724,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, Real* AMREX_RESTRICT uxp = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uyp = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uzp = attribs[PIdx::uz].dataPtr(); - + const auto np = pti.numParticles(); const auto lev = pti.GetLevel(); const auto index = pti.GetPairIndex(); @@ -1740,7 +1740,7 @@ void PhysicalParticleContainer::copy_attribs(WarpXParIter& pti,const Real* xp, xpold[i]=xp[i]; ypold[i]=yp[i]; zpold[i]=zp[i]; - + uxpold[i]=uxp[i]; uypold[i]=uyp[i]; uzpold[i]=uzp[i]; @@ -1780,7 +1780,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real slice_box.setHi(direction, z_max); diagnostic_particles.resize(finestLevel()+1); - + for (int lev = 0; lev < nlevs; ++lev) { const Real* dx = Geom(lev).CellSize(); @@ -1827,7 +1827,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real auto& uzp_old = tmp_particle_data[lev][index][TmpIdx::uzold]; const long np = pti.numParticles(); - + Real uzfrm = -WarpX::gamma_boost*WarpX::beta_boost*PhysConst::c; Real inv_c2 = 1.0/PhysConst::c/PhysConst::c; @@ -1862,7 +1862,7 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real Real uzp = uz_old_p *weight_old + uz_new_p *weight_new; diagnostic_particles[lev][index].GetRealData(DiagIdx::w).push_back(wp[i]); - + diagnostic_particles[lev][index].GetRealData(DiagIdx::x).push_back(xp); diagnostic_particles[lev][index].GetRealData(DiagIdx::y).push_back(yp); diagnostic_particles[lev][index].GetRealData(DiagIdx::z).push_back(zp); @@ -1887,7 +1887,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) AddPlasma(lev, injection_box); } -/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, +/* \brief Gather fields from FArrayBox exfab, eyfab, ezfab, bxfab, byfab, * bzfab into arrays of fields on particles Exp, Eyp, Ezp, Bxp, Byp, Bzp. * \param Exp-Bzp: fields on particles. * \param exfab-bzfab: FAB of electric and magnetic fields for particles in pti @@ -1897,7 +1897,7 @@ PhysicalParticleContainer::ContinuousInjection(const RealBox& injection_box) * \param np_to_gather: number of particles onto which fields are gathered * \param thread_num: if using OpenMP, thread number * \param lev: level on which particles are located - * \param gather_lev: level from which particles gather fields (lev-1) for + * \param gather_lev: level from which particles gather fields (lev-1) for particles in buffers. */ void @@ -1924,14 +1924,14 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - + // If no particles, do not do anything if (np_to_gather == 0) return; // Get cell size on gather_lev const std::array& dx = WarpX::CellSize(std::max(gather_lev,0)); // Set staggering shift depending on e_is_nodal const Real stagger_shift = e_is_nodal ? 0.0 : 0.5; - + // Get box from which field is gathered. // If not gathering from the finest level, the box is coarsened. Box box; @@ -1941,26 +1941,26 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, const IntVect& ref_ratio = WarpX::RefRatio(gather_lev); box = amrex::coarsen(pti.tilebox(),ref_ratio); } - + // Add guard cells to the box. box.grow(ngE); - + const Array4& ex_arr = exfab->array(); const Array4& ey_arr = eyfab->array(); const Array4& ez_arr = ezfab->array(); const Array4& bx_arr = bxfab->array(); const Array4& by_arr = byfab->array(); const Array4& bz_arr = bzfab->array(); - + const Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr() + offset; const Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr() + offset; - + // Lower corner of tile box physical domain const std::array& xyzmin = WarpX::LowerCorner(box, gather_lev); - + const Dim3 lo = lbound(box); - + // Depending on l_lower_in_v and WarpX::nox, call // different versions of template function doGatherShapeN if (WarpX::l_lower_order_in_v){ @@ -2043,10 +2043,10 @@ void PhysicalParticleContainer::InitIonizationModule () } // Compute ADK prefactors (See Chen, JCP 236 (2013), equation (2)) // For now, we assume l=0 and m=0. - // The approximate expressions are used, + // The approximate expressions are used, // without Gamma function Real wa = std::pow(PhysConst::alpha,3) * PhysConst::c / PhysConst::r_e; - Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * + Real Ea = PhysConst::m_e * PhysConst::c*PhysConst::c /PhysConst::q_e * std::pow(PhysConst::alpha,4)/PhysConst::r_e; Real UH = table_ionization_energies[0]; Real l_eff = std::sqrt(UH/ionization_energies[0]) - 1.; @@ -2061,18 +2061,18 @@ void PhysicalParticleContainer::InitIonizationModule () Real C2 = std::pow(2,2*n_eff)/(n_eff*tgamma(n_eff+l_eff+1)*tgamma(n_eff-l_eff)); adk_power[i] = -(2*n_eff - 1); Real Uion = ionization_energies[i]; - adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) + adk_prefactor[i] = dt * wa * C2 * ( Uion/(2*UH) ) * std::pow(2*std::pow((Uion/UH),3./2)*Ea,2*n_eff - 1); adk_exp_prefactor[i] = -2./3 * std::pow( Uion/UH,3./2) * Ea; } } /* \brief create mask of ionized particles (1 if ionized, 0 otherwise) - * + * * \param mfi: tile or grid * \param lev: MR level * \param ionization_mask: Array with as many elements as particles in mfi. - * This function initialized the array, and set each element to 1 or 0 + * This function initialized the array, and set each element to 1 or 0 * depending on whether the particle is ionized or not. */ void @@ -2118,7 +2118,7 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const // Loop over all particles in grid/tile. If ionized, set mask to 1 // and increment ionization level. - ParallelFor( + ParallelFor( np, [=] AMREX_GPU_DEVICE (long i) { // Get index of ionization_level -- cgit v1.2.3 From 309cbaf7dafc84d2ba721b94bbf190b29acaa069 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 17 Sep 2019 09:59:59 +0200 Subject: Optical depths is always plotted if QED is on --- Source/Particles/PhysicalParticleContainer.cpp | 7 +++++++ 1 file changed, 7 insertions(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index cb51f6d96..bbd14bb53 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -96,6 +96,13 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp } } } + + #ifdef WARPX_QED + if(do_qed){ + //Optical depths is always plotted if QED is on + plot_flags[plot_flag_size-1] = 1; + } + #endif } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) -- cgit v1.2.3 From 7f178a55df4d0973f6788d9c27fd10b128ad4571 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Tue, 8 Oct 2019 18:25:10 +0200 Subject: Heavy refactoring --- Source/Particles/MultiParticleContainer.H | 21 +++++++----- Source/Particles/MultiParticleContainer.cpp | 22 +++++++++++++ Source/Particles/PhotonParticleContainer.H | 4 --- Source/Particles/PhotonParticleContainer.cpp | 5 ++- Source/Particles/PhysicalParticleContainer.H | 25 ++++++++++++++- Source/Particles/PhysicalParticleContainer.cpp | 34 ++++++++++++++++++-- Source/Particles/WarpXParticleContainer.H | 15 ++++++++- Source/QED/BreitWheelerEngineWrapper.cpp | 28 ++++++++++++++++ Source/QED/BreitWheelerEngineWrapper.h | 44 ++++++++++++++++++++++++++ Source/QED/Make.package | 12 +++---- Source/QED/QedWrapperCommons.h | 16 ++++++++++ Source/QED/QuantumSyncEngineWrapper.cpp | 29 +++++++++++++++++ Source/QED/QuantumSyncEngineWrapper.h | 44 ++++++++++++++++++++++++++ Source/QED/amrex_rng_wrapper.cpp | 20 ------------ Source/QED/amrex_rng_wrapper.h | 22 ------------- Source/QED/breit_wheeler_engine_wrapper.cpp | 16 ---------- Source/QED/breit_wheeler_engine_wrapper.h | 27 ---------------- Source/QED/qed_wrapper_commons.h | 13 -------- Source/QED/quantum_sync_engine_wrapper.cpp | 16 ---------- Source/QED/quantum_sync_engine_wrapper.h | 27 ---------------- 20 files changed, 275 insertions(+), 165 deletions(-) create mode 100644 Source/QED/BreitWheelerEngineWrapper.cpp create mode 100644 Source/QED/BreitWheelerEngineWrapper.h create mode 100644 Source/QED/QedWrapperCommons.h create mode 100644 Source/QED/QuantumSyncEngineWrapper.cpp create mode 100644 Source/QED/QuantumSyncEngineWrapper.h delete mode 100644 Source/QED/amrex_rng_wrapper.cpp delete mode 100644 Source/QED/amrex_rng_wrapper.h delete mode 100644 Source/QED/breit_wheeler_engine_wrapper.cpp delete mode 100644 Source/QED/breit_wheeler_engine_wrapper.h delete mode 100644 Source/QED/qed_wrapper_commons.h delete mode 100644 Source/QED/quantum_sync_engine_wrapper.cpp delete mode 100644 Source/QED/quantum_sync_engine_wrapper.h (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 8869e4d4a..6bf53d6ec 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -23,8 +23,8 @@ // #ifdef WARPX_QED - #include "breit_wheeler_engine_wrapper.h" - #include "quantum_sync_engine_wrapper.h" + #include "BreitWheelerEngineWrapper.h" + #include "QuantumSyncEngineWrapper.h" #endif class MultiParticleContainer @@ -196,12 +196,6 @@ public: PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } - -#ifdef WARPX_QED - warpx_breit_wheeler_engine bw_engine; - warpx_quantum_sync_engine qs_engine; -#endif - protected: // Particle container types @@ -219,6 +213,17 @@ protected: std::vector species_types; +#ifdef WARPX_QED + // The QED engines + BreitWheelerEngine bw_engine; + QuantumSynchrotronEngine qs_engine; + //_______________________________ + + //Initialize QED engines and provides smart pointers + //to species who need QED processes + void InitQED (); +#endif + private: // physical particles (+ laser) diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 715c97b99..a5632dede 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -1,6 +1,7 @@ #include #include #include +#include #include #include @@ -149,6 +150,11 @@ MultiParticleContainer::InitData () // For each species, get the ID of its product species. // This is used for ionization and pair creation processes. mapSpeciesProduct(); + +#ifdef WARPX_QED + InitQED(); +#endif + } @@ -726,3 +732,19 @@ MultiParticleContainer::doFieldIonization () } // lev } // pc_source } + +#ifdef WARPX_QED +void MultiParticleContainer::InitQED () +{ + for (auto& pc : allcontainers) { + if(pc->has_quantum_sync()){ + pc->set_quantum_sync_engine_ptr + (std::make_shared(qs_engine)); + } + if(pc->has_breit_wheeler()){ + pc->set_breit_wheeler_engine_ptr + (std::make_shared(bw_engine)); + } + } +} +#endif diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 98e8e0aab..bda539419 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -68,13 +68,9 @@ public: private: #ifdef WARPX_QED - // A flag to enable breit_wheeler process for photons - bool do_qed_breit_wheeler = false; - // A function to initialize the Tau component according to the BW engine void InitTauBreitWheeler(); #endif - }; #endif // #ifndef WARPX_PhotonParticleContainer_H_ diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index a8a51723e..9ed2e10e6 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -128,6 +128,9 @@ PhotonParticleContainer::Evolve (int lev, void PhotonParticleContainer::InitTauBreitWheeler() { BL_PROFILE("PhotonParticleContainer::InitOpticalDepth"); + //Get functor + auto get_opt = shr_ptr_bw_engine->build_optical_depth_functor(); + //Looping over all the particles int num_levels = finestLevel() + 1; for (int lev=0; lev < num_levels; ++lev) @@ -136,7 +139,7 @@ void PhotonParticleContainer::InitTauBreitWheeler() amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - taus[i] = warpx_breit_wheeler_engine::get_optical_depth(); + taus[i] = get_opt(); } ); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 79685c08b..4fa621514 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -9,6 +9,11 @@ #include +#ifdef WARPX_QED +#include +#include +#endif + class PhysicalParticleContainer : public WarpXParticleContainer { @@ -180,6 +185,16 @@ public: amrex::FArrayBox const * & ezfab, amrex::FArrayBox const * & bxfab, amrex::FArrayBox const * & byfab, amrex::FArrayBox const * & bzfab); +#ifdef WARPX_QED + bool has_quantum_sync() override; + bool has_breit_wheeler() override; + + void set_breit_wheeler_engine_ptr + (std::shared_ptr) override; + void set_quantum_sync_engine_ptr + (std::shared_ptr) override; +#endif + protected: std::string species_name; @@ -194,13 +209,21 @@ protected: // Inject particles during the whole simulation void ContinuousInjection (const amrex::RealBox& injection_box) override; -private: #ifdef WARPX_QED // A flag to enable quantum_synchrotron process for leptons bool do_qed_quantum_sync = false; + // A flag to enable breit_wheeler process [photons only!!] + bool do_qed_breit_wheeler = false; + // A function to initialize the Tau component according to the QS engine void InitTauQuantumSync(); + + // A smart pointer to an instance of a Quantum Synchrotron engine + std::shared_ptr shr_ptr_qs_engine; + + // A smart pointer to an instance of a Breit Wheeler engine [photons only!] + std::shared_ptr shr_ptr_bw_engine; #endif }; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 9e12eb26e..ac5b98dae 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1,6 +1,8 @@ #include #include +#include + #include #include #include @@ -2110,11 +2112,39 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const } #ifdef WARPX_QED + +bool PhysicalParticleContainer::has_quantum_sync() +{ + return do_qed_quantum_sync; +} + +bool PhysicalParticleContainer::has_breit_wheeler() +{ + return do_qed_breit_wheeler; +} + +void +PhysicalParticleContainer:: +set_breit_wheeler_engine_ptr(std::shared_ptr ptr) +{ + shr_ptr_bw_engine = ptr; +} + +void +PhysicalParticleContainer:: +set_quantum_sync_engine_ptr(std::shared_ptr ptr) +{ + shr_ptr_qs_engine = ptr; +} + // A function to initialize the Tau component according to the QS engine void PhysicalParticleContainer::InitTauQuantumSync() { BL_PROFILE("PhysicalParticleContainer::InitTauQuantumSync"); -//Looping over all the particles + //Get functor + auto get_opt = shr_ptr_qs_engine->build_optical_depth_functor(); + + //Looping over all the particles int num_levels = finestLevel() + 1; for (int lev=0; lev < num_levels; ++lev) for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti){ @@ -2122,7 +2152,7 @@ void PhysicalParticleContainer::InitTauQuantumSync() amrex::ParallelFor( pti.numParticles(), [=] AMREX_GPU_DEVICE (long i) { - taus[i] = warpx_quantum_sync_engine::get_optical_depth(); + taus[i] = get_opt(); } ); } diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 4b36b12f6..3f327c969 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -1,12 +1,17 @@ #ifndef WARPX_WarpXParticleContainer_H_ #define WARPX_WarpXParticleContainer_H_ +#include + #include "WarpXDtType.H" #include #include -#include +#ifdef WARPX_QED +#include +#include +#endif enum struct ConvertDirection{WarpX_to_SI, SI_to_WarpX}; @@ -315,6 +320,14 @@ protected: #ifdef WARPX_QED bool do_qed = false; + + virtual bool has_quantum_sync(){return false;}; + virtual bool has_breit_wheeler(){return false;}; + + virtual void + set_breit_wheeler_engine_ptr(std::shared_ptr){}; + virtual void + set_quantum_sync_engine_ptr(std::shared_ptr){}; #endif amrex::Vector local_rho; diff --git a/Source/QED/BreitWheelerEngineWrapper.cpp b/Source/QED/BreitWheelerEngineWrapper.cpp new file mode 100644 index 000000000..b654cc9b0 --- /dev/null +++ b/Source/QED/BreitWheelerEngineWrapper.cpp @@ -0,0 +1,28 @@ +#include "BreitWheelerEngineWrapper.h" +//This file provides a wrapper aroud the breit_wheeler engine +//provided by the PICSAR library + +using namespace picsar::multi_physics; + +// Functors ================================== + +// Initialization of the optical depth +AMREX_GPU_DEVICE +amrex::Real +BreitWheelerGetOpticalDepth::operator() () const +{ + return WarpXBreitWheelerWrapper:: + internal_get_optical_depth(amrex::Random()); +} +//____________________________________________ + +// Factory class ============================= + +BreitWheelerEngine::BreitWheelerEngine(){} + +BreitWheelerGetOpticalDepth BreitWheelerEngine::build_optical_depth_functor() +{ + return BreitWheelerGetOpticalDepth(); +} + +//============================================ diff --git a/Source/QED/BreitWheelerEngineWrapper.h b/Source/QED/BreitWheelerEngineWrapper.h new file mode 100644 index 000000000..6aacf4f4e --- /dev/null +++ b/Source/QED/BreitWheelerEngineWrapper.h @@ -0,0 +1,44 @@ +#ifndef WARPX_breit_wheeler_engine_wrapper_h_ +#define WARPX_breit_wheeler_engine_wrapper_h_ + +//This file provides a wrapper aroud the breit_wheeler engine +//provided by the QED modules of the PICSAR library + +#include "QedWrapperCommons.h" + +//BW ENGINE from PICSAR +#include "breit_wheeler_engine.hpp" + +using WarpXBreitWheelerWrapper = + picsar::multi_physics::breit_wheeler_engine; + +// Functors ================================== + +// These functors provide the core elementary functions of the library +// Can be included in GPU kernels + +// Initialization of the optical depth +class BreitWheelerGetOpticalDepth +{ +public: + BreitWheelerGetOpticalDepth() + {}; + + AMREX_GPU_DEVICE + amrex::Real operator() () const; +}; +//____________________________________________ + +// Factory class ============================= +class BreitWheelerEngine +{ +public: + BreitWheelerEngine(); + + //Builds the functor to initialize the optical depth + BreitWheelerGetOpticalDepth build_optical_depth_functor(); +}; + +//============================================ + +#endif //WARPX_breit_wheeler_engine_wrapper_H_ diff --git a/Source/QED/Make.package b/Source/QED/Make.package index 5c71d44a2..7af5860ee 100644 --- a/Source/QED/Make.package +++ b/Source/QED/Make.package @@ -1,10 +1,8 @@ -CEXE_headers += qed_wrapper_commons.h -CEXE_headers += amrex_rng_wrapper.h -CEXE_headers += breit_wheeler_engine_wrapper.h -CEXE_headers += quantum_syncs_engine_wrapper.h -CEXE_sources += amrex_rng_wrapper.cpp -CEXE_sources += breit_wheeler_engine_wrapper.cpp -CEXE_sources += quantum_sync_engine_wrapper.cpp +CEXE_headers += QedWrapperCommons.h +CEXE_headers += BreitWheelerEngineWrapper.h +CEXE_headers += QuantumSyncsEngineWrapper.h +CEXE_sources += BreitWheelerEngineWrapper.cpp +CEXE_sources += QuantumSyncEngineWrapper.cpp INCLUDE_LOCATIONS += $(WARPX_HOME)/Source/QED VPATH_LOCATIONS += $(WARPX_HOME)/Source/QED diff --git a/Source/QED/QedWrapperCommons.h b/Source/QED/QedWrapperCommons.h new file mode 100644 index 000000000..821034c06 --- /dev/null +++ b/Source/QED/QedWrapperCommons.h @@ -0,0 +1,16 @@ +#ifndef WARPX_amrex_qed_wrapper_commons_h_ +#define WARPX_amrex_qed_wrapper_commons_h_ + +//Common definitions for the QED library wrappers + +#include + +//Sets the decorator for GPU +#define PXRMP_GPU AMREX_GPU_DEVICE +//Sets SI units in the library +#define PXRMP_WITH_SI_UNITS + +//An empty data type +struct DummyStruct{}; + +#endif //WARPX_amrex_qed_wrapper_commons_h_ diff --git a/Source/QED/QuantumSyncEngineWrapper.cpp b/Source/QED/QuantumSyncEngineWrapper.cpp new file mode 100644 index 000000000..ce130bb62 --- /dev/null +++ b/Source/QED/QuantumSyncEngineWrapper.cpp @@ -0,0 +1,29 @@ +#include "QuantumSyncEngineWrapper.h" +//This file provides a wrapper aroud the quantum_sync engine +//provided by the PICSAR library + +using namespace picsar::multi_physics; + +// Functors ================================== + +// Initialization of the optical depth + +AMREX_GPU_DEVICE +amrex::Real +QuantumSynchrotronGetOpticalDepth::operator() () const +{ + return WarpXQuantumSynchrotronWrapper:: + internal_get_optical_depth(amrex::Random()); +} +//____________________________________________ + +// Factory class ============================= + +QuantumSynchrotronEngine::QuantumSynchrotronEngine(){} + +QuantumSynchrotronGetOpticalDepth QuantumSynchrotronEngine::build_optical_depth_functor() +{ + return QuantumSynchrotronGetOpticalDepth(); +} + +//============================================ diff --git a/Source/QED/QuantumSyncEngineWrapper.h b/Source/QED/QuantumSyncEngineWrapper.h new file mode 100644 index 000000000..0cb444c07 --- /dev/null +++ b/Source/QED/QuantumSyncEngineWrapper.h @@ -0,0 +1,44 @@ +#ifndef WARPX_quantum_sync_engine_wrapper_h_ +#define WARPX_quantum_sync_engine_wrapper_h_ + +//This file provides a wrapper aroud the breit_wheeler engine +//provided by the QED modules of the PICSAR library + +#include "QedWrapperCommons.h" + +//QS ENGINE from PICSAR +#include "quantum_sync_engine.hpp" + +using WarpXQuantumSynchrotronWrapper = + picsar::multi_physics::quantum_synchrotron_engine; + +// Functors ================================== + +// These functors provide the core elementary functions of the library +// Can be included in GPU kernels + +// Initialization of the optical depth +class QuantumSynchrotronGetOpticalDepth +{ +public: + QuantumSynchrotronGetOpticalDepth() + {}; + + AMREX_GPU_DEVICE + amrex::Real operator() () const; +}; +//____________________________________________ + +// Factory class ============================= +class QuantumSynchrotronEngine +{ +public: + QuantumSynchrotronEngine(); + + //Builds the functor to initialize the optical depth + QuantumSynchrotronGetOpticalDepth build_optical_depth_functor(); +}; + +//============================================ + +#endif //WARPX_quantum_sync_engine_wrapper_h_ diff --git a/Source/QED/amrex_rng_wrapper.cpp b/Source/QED/amrex_rng_wrapper.cpp deleted file mode 100644 index 10fe02c1d..000000000 --- a/Source/QED/amrex_rng_wrapper.cpp +++ /dev/null @@ -1,20 +0,0 @@ -//This file provides a wrapper aroud the RNG -//provided by the amrex library - -#include "amrex_rng_wrapper.h" - -//RNG wrapper BW engine -//Get rnd number uniformly distributed in [a,b) -amrex::Real AMREX_GPU_DEVICE -amrex_rng_wrapper::unf(amrex::Real a, amrex::Real b) -{ - return (b-a)*amrex::Random() + a; -} - -//Get rnd number with exponential distribution -amrex::Real AMREX_GPU_DEVICE -amrex_rng_wrapper::exp(amrex::Real l) -{ - amrex::Real zero_plus_to_one = 1.0 - unf(0.0, 1.0); - return -log(zero_plus_to_one)/l; -} diff --git a/Source/QED/amrex_rng_wrapper.h b/Source/QED/amrex_rng_wrapper.h deleted file mode 100644 index e846cb44f..000000000 --- a/Source/QED/amrex_rng_wrapper.h +++ /dev/null @@ -1,22 +0,0 @@ -#ifndef WARPX_amrex_rng_wrapper_h_ -#define WARPX_amrex_rng_wrapper_h_ - -//This file provides a wrapper aroud the RNG -//provided by the amrex library - -#include - -//RNG wrapper BW engine -class amrex_rng_wrapper -{ -public: - //Get rnd number uniformly distributed in [a,b) - amrex::Real AMREX_GPU_DEVICE - unf(amrex::Real a, amrex::Real b); - - //Get rnd number with exponential distribution - amrex::Real AMREX_GPU_DEVICE - exp(amrex::Real l); -}; - -#endif //WARPX_amrex_rng_wrapper_h_ diff --git a/Source/QED/breit_wheeler_engine_wrapper.cpp b/Source/QED/breit_wheeler_engine_wrapper.cpp deleted file mode 100644 index 846c76099..000000000 --- a/Source/QED/breit_wheeler_engine_wrapper.cpp +++ /dev/null @@ -1,16 +0,0 @@ -#include "breit_wheeler_engine_wrapper.h" -//This file provides a wrapper aroud the breit_wheeler engine -//provided by the standard template library - -using namespace picsar::multi_physics; - -warpx_breit_wheeler_engine::warpx_breit_wheeler_engine(): - breit_wheeler_engine(std::move(amrex_rng_wrapper{})) -{}; - -//Interface for the get_optical_depth method of the BW engine -amrex::Real -AMREX_GPU_DEVICE -warpx_breit_wheeler_engine::get_optical_depth(){ - return internal_get_optical_depth(amrex::Random()); -} diff --git a/Source/QED/breit_wheeler_engine_wrapper.h b/Source/QED/breit_wheeler_engine_wrapper.h deleted file mode 100644 index 390717b09..000000000 --- a/Source/QED/breit_wheeler_engine_wrapper.h +++ /dev/null @@ -1,27 +0,0 @@ -#ifndef WARPX_breit_wheeler_engine_wrapper_h_ -#define WARPX_breit_wheeler_engine_wrapper_h_ - -//This file provides a wrapper aroud the breit_wheeler engine -//provided by the QED modules of the PICSAR library - -//BW ENGINE -#include "qed_wrapper_commons.h" -#include "breit_wheeler_engine.hpp" - -#include "amrex_rng_wrapper.h" - -class warpx_breit_wheeler_engine : - public picsar::multi_physics:: - breit_wheeler_engine -{ -public: - warpx_breit_wheeler_engine(); - - //Interface for the get_optical_depth method of the BW engine - static AMREX_GPU_DEVICE - amrex::Real get_optical_depth(); -}; - -//___________________________________________ - -#endif //WARPX_breit_wheeler_engine_wrapper_H_ diff --git a/Source/QED/qed_wrapper_commons.h b/Source/QED/qed_wrapper_commons.h deleted file mode 100644 index 49f82ad5b..000000000 --- a/Source/QED/qed_wrapper_commons.h +++ /dev/null @@ -1,13 +0,0 @@ -#ifndef WARPX_amrex_qed_wrapper_commons_h_ -#define WARPX_amrex_qed_wrapper_commons_h_ - -//Common definitions for the QED library wrappers - -#include - -//Sets the decorator for GPU -#define PXRMP_GPU AMREX_GPU_DEVICE -//Sets SI units in the library -#define PXRMP_WITH_SI_UNITS - -#endif diff --git a/Source/QED/quantum_sync_engine_wrapper.cpp b/Source/QED/quantum_sync_engine_wrapper.cpp deleted file mode 100644 index c41cd4b61..000000000 --- a/Source/QED/quantum_sync_engine_wrapper.cpp +++ /dev/null @@ -1,16 +0,0 @@ -#include "quantum_sync_engine_wrapper.h" -//This file provides a wrapper aroud the quantum_sync engine -//provided by the standard template library - -using namespace picsar::multi_physics; - -warpx_quantum_sync_engine::warpx_quantum_sync_engine(): - quantum_synchrotron_engine(std::move(amrex_rng_wrapper{})) -{}; - -//Interface for the get_optical_depth method of the QS engine -amrex::Real -AMREX_GPU_DEVICE -warpx_quantum_sync_engine::get_optical_depth(){ - return internal_get_optical_depth(amrex::Random()); -} diff --git a/Source/QED/quantum_sync_engine_wrapper.h b/Source/QED/quantum_sync_engine_wrapper.h deleted file mode 100644 index 0af2c2189..000000000 --- a/Source/QED/quantum_sync_engine_wrapper.h +++ /dev/null @@ -1,27 +0,0 @@ -#ifndef WARPX_quantum_sync_engine_wrapper_h_ -#define WARPX_quantum_sync_engine_wrapper_h_ - -//This file provides a wrapper aroud the breit_wheeler engine -//provided by the QED modules of the PICSAR library - -//BW ENGINE -#include "qed_wrapper_commons.h" -#include "quantum_sync_engine.hpp" - -#include "amrex_rng_wrapper.h" - -class warpx_quantum_sync_engine : - public picsar::multi_physics:: - quantum_synchrotron_engine -{ -public: - warpx_quantum_sync_engine(); - - //Interface for the get_optical_depth method of the BW engine - static AMREX_GPU_DEVICE - amrex::Real get_optical_depth(); -}; - -//___________________________________________ - -#endif //WARPX_quantum_sync_engine_wrapper_h_ -- cgit v1.2.3 From 5817971d808a76cbfd79c95ff964af0186615472 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Wed, 9 Oct 2019 18:08:19 +0200 Subject: Added tests for the initialization of the optical depth (disabled since they would need PICSAR on the QED branch) --- .../Modules/qed/breit_wheeler/check_2d_tau_init.py | 30 +++++++ .../qed/breit_wheeler/inputs.2d_test_tau_init | 68 ++++++++++++++++ .../qed/quantum_synchrotron/check_2d_tau_init.py | 38 +++++++++ .../quantum_synchrotron/inputs.2d_test_tau_init | 92 ++++++++++++++++++++++ Regression/WarpX-tests.ini | 32 ++++++++ Source/Particles/PhysicalParticleContainer.cpp | 1 - 6 files changed, 260 insertions(+), 1 deletion(-) create mode 100755 Examples/Modules/qed/breit_wheeler/check_2d_tau_init.py create mode 100644 Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init create mode 100755 Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py create mode 100644 Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Modules/qed/breit_wheeler/check_2d_tau_init.py b/Examples/Modules/qed/breit_wheeler/check_2d_tau_init.py new file mode 100755 index 000000000..98f37e5e3 --- /dev/null +++ b/Examples/Modules/qed/breit_wheeler/check_2d_tau_init.py @@ -0,0 +1,30 @@ +#! /usr/bin/env python3 +import yt +import numpy as np +import scipy.stats as st +import sys + +# This script checks if photons initialized with Breit Wheeler process enabled +# do actually have an exponentially distributed optical depth + +# Tolerance +tol = 1e-2 + +def check(): + filename = sys.argv[1] + data_set = yt.load(filename) + + all_data = data_set.all_data() + res_tau = all_data["photons", 'particle_tau'] + + loc, scale = st.expon.fit(res_tau) + + # loc should be very close to 0, scale should be very close to 1 + assert(np.abs(loc - 0) < tol) + assert(np.abs(scale - 1) < tol) + +def main(): + check() + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init b/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init new file mode 100644 index 000000000..9eb29ad1b --- /dev/null +++ b/Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init @@ -0,0 +1,68 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 0 +amr.n_cell = 128 128 +amr.max_grid_size = 128 # maximum size of each AMReX box, used to decompose the domain +amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain +amr.plot_int = 10 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 # Is periodic? +geometry.prob_lo = -32.e-6 -32.e-6 # physical domain +geometry.prob_hi = 32.e-6 32.e-6 +amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported) +warpx.fine_tag_lo = -5.e-6 -35.e-6 +warpx.fine_tag_hi = 5.e-6 -25.e-6 + +################################# +############ NUMERICS ########### +################################# +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = standard +algo.particle_pusher = boris +interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z +interpolation.noy = 3 +interpolation.noz = 3 +warpx.verbose = 1 +warpx.do_dive_cleaning = 0 +warpx.plot_raw_fields = 0 +warpx.plot_raw_fields_guards = 0 +warpx.plot_finepatch = 0 +warpx.plot_crsepatch = 0 +warpx.use_filter = 1 +warpx.cfl = 1. # if 1., the time step is set to its CFL limit +warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition + +################################# +############ PLASMA ############# +################################# +particles.nspecies = 1 # number of species +particles.species_names = photons +particles.photon_species = photons +################################# + +photons.charge = -q_e +photons.mass = m_e +photons.injection_style = "NUniformPerCell" +photons.profile = "constant" +photons.xmin = -30e-6 +photons.ymin = -30e-6 +photons.zmin = -30e-6 +photons.xmax = 30e-6 +photons.ymax = 30e-6 +photons.zmax = 30e-6 +photons.num_particles_per_cell_each_dim = 2 2 +photons.density = 1e19 +photons.profile = "constant" +photons.momentum_distribution_type = "gaussian" +photons.ux_m = 0.0 +photons.uy_m = 0.0 +photons.uz_m = 0.0 +photons.ux_th = 100. +photons.uy_th = 100. +photons.uz_th = 100. +##########QED#################### +photons.do_qed = 1 +photons.do_qed_breit_wheeler = 1 +################################# diff --git a/Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py b/Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py new file mode 100755 index 000000000..65b441719 --- /dev/null +++ b/Examples/Modules/qed/quantum_synchrotron/check_2d_tau_init.py @@ -0,0 +1,38 @@ +#! /usr/bin/env python3 +import yt +import numpy as np +import scipy.stats as st +import sys + +# This script checks if electrons and positrons initialized with +# Quantum Synchrotron process enabled +# do actually have an exponentially distributed optical depth + +# Tolerance +tol = 1e-2 + +def check(): + filename = sys.argv[1] + data_set = yt.load(filename) + + all_data = data_set.all_data() + res_ele_tau = all_data["electrons", 'particle_tau'] + res_pos_tau = all_data["positrons", 'particle_tau'] + + loc_ele, scale_ele = st.expon.fit(res_ele_tau) + loc_pos, scale_pos = st.expon.fit(res_pos_tau) + + print(loc_ele, scale_ele) + print(loc_pos, scale_pos) + + # loc should be very close to 0, scale should be very close to 1 + assert(np.abs(loc_ele - 0) < tol) + assert(np.abs(loc_pos - 0) < tol) + assert(np.abs(scale_ele - 1) < tol) + assert(np.abs(scale_pos - 1) < tol) + +def main(): + check() + +if __name__ == "__main__": + main() \ No newline at end of file diff --git a/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init b/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init new file mode 100644 index 000000000..7d26aee0c --- /dev/null +++ b/Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init @@ -0,0 +1,92 @@ +################################# +####### GENERAL PARAMETERS ###### +################################# +max_step = 0 +amr.n_cell = 128 128 +amr.max_grid_size = 128 # maximum size of each AMReX box, used to decompose the domain +amr.blocking_factor = 32 # minimum size of each AMReX box, used to decompose the domain +amr.plot_int = 10 +geometry.coord_sys = 0 # 0: Cartesian +geometry.is_periodic = 0 0 # Is periodic? +geometry.prob_lo = -32.e-6 -32.e-6 # physical domain +geometry.prob_hi = 32.e-6 32.e-6 +amr.max_level = 0 # Maximum level in hierarchy (1 might be unstable, >1 is not supported) +warpx.fine_tag_lo = -5.e-6 -35.e-6 +warpx.fine_tag_hi = 5.e-6 -25.e-6 + +################################# +############ NUMERICS ########### +################################# +algo.current_deposition = esirkepov +algo.charge_deposition = standard +algo.field_gathering = standard +algo.particle_pusher = boris +interpolation.nox = 3 # Particle interpolation order. Must be the same in x, y, and z +interpolation.noy = 3 +interpolation.noz = 3 +warpx.verbose = 1 +warpx.do_dive_cleaning = 0 +warpx.plot_raw_fields = 0 +warpx.plot_raw_fields_guards = 0 +warpx.plot_finepatch = 0 +warpx.plot_crsepatch = 0 +warpx.use_filter = 1 +warpx.cfl = 1. # if 1., the time step is set to its CFL limit +warpx.do_pml = 1 # use Perfectly Matched Layer as boundary condition + +################################# +############ PLASMA ############# +################################# +particles.nspecies = 2 # number of species +particles.species_names = electrons positrons +################################# + +electrons.charge = -q_e +electrons.mass = m_e +electrons.injection_style = "NUniformPerCell" +electrons.profile = "constant" +electrons.xmin = -30e-6 +electrons.ymin = -30e-6 +electrons.zmin = -30e-6 +electrons.xmax = 30e-6 +electrons.ymax = 30e-6 +electrons.zmax = 30e-6 +electrons.num_particles_per_cell_each_dim = 2 2 +electrons.density = 1e19 +electrons.profile = "constant" +electrons.momentum_distribution_type = "gaussian" +electrons.ux_m = 0.0 +electrons.uy_m = 0.0 +electrons.uz_m = 0.0 +electrons.ux_th = 100. +electrons.uy_th = 100. +electrons.uz_th = 100. +##########QED#################### +electrons.do_qed = 1 +electrons.do_qed_quantum_sync = 1 +################################# + +positrons.charge = q_e +positrons.mass = m_e +positrons.injection_style = "NUniformPerCell" +positrons.profile = "constant" +positrons.xmin = -30e-6 +positrons.ymin = -30e-6 +positrons.zmin = -30e-6 +positrons.xmax = 30e-6 +positrons.ymax = 30e-6 +positrons.zmax = 30e-6 +positrons.num_particles_per_cell_each_dim = 2 2 +positrons.density = 1e19 +positrons.profile = "constant" +positrons.momentum_distribution_type = "gaussian" +positrons.ux_m = 0.0 +positrons.uy_m = 0.0 +positrons.uz_m = 0.0 +positrons.ux_th = 100. +positrons.uy_th = 100. +positrons.uz_th = 100. +##########QED#################### +positrons.do_qed = 1 +positrons.do_qed_quantum_sync = 1 +################################# diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini index e40f7ee54..a5805f566 100644 --- a/Regression/WarpX-tests.ini +++ b/Regression/WarpX-tests.ini @@ -632,3 +632,35 @@ compileTest = 0 doVis = 0 compareParticles = 0 analysisRoutine = Examples/Tests/photon_pusher/check.py + +# NEEDS PICSAR on the QED branch +#[breit_wheeler_tau_init] +#buildDir = . +#inputFile = Examples/Modules/qed/breit_wheeler/inputs.2d_test_tau_init +#dim = 2 +#addToCompileString = QED=TRUE +#restartTest = 0 +#useMPI = 1 +#numprocs = 2 +#useOMP = 1 +#numthreads = 2 +#compileTest = 0 +#doVis = 0 +#compareParticles = 0 +#analysisRoutine = Examples/Modules/qed/breit_wheeler/check.py + +# NEEDS PICSAR on the QED branch +#[quantum_sync_tau_init] +#buildDir = . +#inputFile = Examples/Modules/qed/quantum_synchrotron/inputs.2d_test_tau_init +#dim = 2 +#addToCompileString = QED=TRUE +#restartTest = 0 +#useMPI = 1 +#numprocs = 2 +#useOMP = 1 +#numthreads = 2 +#compileTest = 0 +#doVis = 0 +#compareParticles = 0 +#analysisRoutine = Examples/Modules/qed/quantum_synchrotron/check.py \ No newline at end of file diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ac5b98dae..3d3d3ff74 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2143,7 +2143,6 @@ void PhysicalParticleContainer::InitTauQuantumSync() BL_PROFILE("PhysicalParticleContainer::InitTauQuantumSync"); //Get functor auto get_opt = shr_ptr_qs_engine->build_optical_depth_functor(); - //Looping over all the particles int num_levels = finestLevel() + 1; for (int lev=0; lev < num_levels; ++lev) -- cgit v1.2.3 From 7df331a7c81899600240aab5954b8950d01ba4f5 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 11 Oct 2019 15:00:57 +0200 Subject: Moved tau initialization in AddPlasma for both photons and physical particles --- Source/Particles/PhotonParticleContainer.H | 7 --- Source/Particles/PhotonParticleContainer.cpp | 28 ------------ Source/Particles/PhysicalParticleContainer.H | 3 -- Source/Particles/PhysicalParticleContainer.cpp | 59 +++++++++++++++----------- 4 files changed, 34 insertions(+), 63 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhotonParticleContainer.H b/Source/Particles/PhotonParticleContainer.H index 7c0da9331..407cf26f3 100644 --- a/Source/Particles/PhotonParticleContainer.H +++ b/Source/Particles/PhotonParticleContainer.H @@ -73,13 +73,6 @@ public: int lev, int depos_lev, amrex::Real dt) override {}; - -private: - -#ifdef WARPX_QED - // A function to initialize the Tau component according to the BW engine - void InitTauBreitWheeler(); -#endif }; #endif // #ifndef WARPX_PhotonParticleContainer_H_ diff --git a/Source/Particles/PhotonParticleContainer.cpp b/Source/Particles/PhotonParticleContainer.cpp index 9ed2e10e6..3c70a957f 100644 --- a/Source/Particles/PhotonParticleContainer.cpp +++ b/Source/Particles/PhotonParticleContainer.cpp @@ -45,11 +45,6 @@ void PhotonParticleContainer::InitData() { AddParticles(0); // Note - add on level 0 -#ifdef WARPX_QED - if(do_qed_breit_wheeler) - InitTauBreitWheeler(); -#endif - if (maxLevel() > 0) { Redistribute(); // We then redistribute } @@ -122,26 +117,3 @@ PhotonParticleContainer::Evolve (int lev, t, dt); } - -#ifdef WARPX_QED -// A function to initialize the Tau component according to the BW engine -void PhotonParticleContainer::InitTauBreitWheeler() -{ - BL_PROFILE("PhotonParticleContainer::InitOpticalDepth"); - //Get functor - auto get_opt = shr_ptr_bw_engine->build_optical_depth_functor(); - - //Looping over all the particles - int num_levels = finestLevel() + 1; - for (int lev=0; lev < num_levels; ++lev) - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti){ - auto taus = pti.GetAttribs(particle_comps["tau"]).dataPtr(); - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - taus[i] = get_opt(); - } - ); - } -} -#endif diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ff98d1de9..84f46de5c 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -216,9 +216,6 @@ protected: // A flag to enable breit_wheeler process [photons only!!] bool do_qed_breit_wheeler = false; - // A function to initialize the Tau component according to the QS engine - void InitTauQuantumSync(); - // A smart pointer to an instance of a Quantum Synchrotron engine std::shared_ptr shr_ptr_qs_engine; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 3d3d3ff74..357fb1f68 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -124,11 +124,6 @@ void PhysicalParticleContainer::InitData() AddParticles(0); // Note - add on level 0 -#ifdef WARPX_QED - if(do_qed_quantum_sync) - InitTauQuantumSync(); -#endif - Redistribute(); // We then redistribute } @@ -513,6 +508,30 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pi = soa.GetIntData(particle_icomps["ionization_level"]).data() + old_size; } +#ifdef WARPX_QED + //Pointer to the optical depth component + amrex::Real* p_tau; + + //If a QED effect is enabled, tau has to be initialized + bool loc_has_quantum_sync = has_quantum_sync(); + bool loc_has_breit_wheeler = has_breit_wheeler(); + if(loc_has_quantum_sync || loc_has_breit_wheeler){ + p_tau = soa.GetRealData(particle_comps["tau"]).data() + old_size; + } + + //If needed, get the appropriate functors from the engines + QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; + BreitWheelerGetOpticalDepth breit_wheeler_get_opt; + if(loc_has_quantum_sync){ + quantum_sync_get_opt = + shr_ptr_qs_engine->build_optical_depth_functor(); + } + if(loc_has_breit_wheeler){ + breit_wheeler_get_opt = + shr_ptr_bw_engine->build_optical_depth_functor(); + } +#endif + const GpuArray overlap_corner {AMREX_D_DECL(overlap_realbox.lo(0), overlap_realbox.lo(1), @@ -660,6 +679,16 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) pi[ip] = loc_ionization_initial_level; } +#ifdef WARPX_QED + if(loc_has_quantum_sync){ + p_tau[ip] = quantum_sync_get_opt(); + } + + if(loc_has_breit_wheeler){ + p_tau[ip] = breit_wheeler_get_opt(); + } +#endif + u.x *= PhysConst::c; u.y *= PhysConst::c; u.z *= PhysConst::c; @@ -2136,24 +2165,4 @@ set_quantum_sync_engine_ptr(std::shared_ptr ptr) { shr_ptr_qs_engine = ptr; } - -// A function to initialize the Tau component according to the QS engine -void PhysicalParticleContainer::InitTauQuantumSync() -{ - BL_PROFILE("PhysicalParticleContainer::InitTauQuantumSync"); - //Get functor - auto get_opt = shr_ptr_qs_engine->build_optical_depth_functor(); - //Looping over all the particles - int num_levels = finestLevel() + 1; - for (int lev=0; lev < num_levels; ++lev) - for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti){ - auto taus = pti.GetAttribs(particle_comps["tau"]).dataPtr(); - amrex::ParallelFor( - pti.numParticles(), - [=] AMREX_GPU_DEVICE (long i) { - taus[i] = get_opt(); - } - ); - } -} #endif -- cgit v1.2.3 From 2cbb00e200bf52fdfc982d5c71ab8e54624bf9a1 Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 11 Oct 2019 15:21:57 +0200 Subject: Style changes --- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 357fb1f68..54e3f409f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -523,11 +523,11 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) QuantumSynchrotronGetOpticalDepth quantum_sync_get_opt; BreitWheelerGetOpticalDepth breit_wheeler_get_opt; if(loc_has_quantum_sync){ - quantum_sync_get_opt = + quantum_sync_get_opt = shr_ptr_qs_engine->build_optical_depth_functor(); } if(loc_has_breit_wheeler){ - breit_wheeler_get_opt = + breit_wheeler_get_opt = shr_ptr_bw_engine->build_optical_depth_functor(); } #endif -- cgit v1.2.3 From 178103ad97018f2566225c180f631e04edd01cec Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 11 Oct 2019 15:26:35 +0200 Subject: Style changes --- Source/Particles/PhysicalParticleContainer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 54e3f409f..32c5e5283 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -686,7 +686,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if(loc_has_breit_wheeler){ p_tau[ip] = breit_wheeler_get_opt(); - } + } #endif u.x *= PhysConst::c; -- cgit v1.2.3 From c6ba1eb9bdae6914f235211b65f5245c5f5c07ef Mon Sep 17 00:00:00 2001 From: Luca Fedeli Date: Fri, 11 Oct 2019 15:35:27 +0200 Subject: Style changes --- Source/Particles/PhysicalParticleContainer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 32c5e5283..d10803320 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -511,7 +511,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) #ifdef WARPX_QED //Pointer to the optical depth component amrex::Real* p_tau; - + //If a QED effect is enabled, tau has to be initialized bool loc_has_quantum_sync = has_quantum_sync(); bool loc_has_breit_wheeler = has_breit_wheeler(); @@ -529,7 +529,7 @@ PhysicalParticleContainer::AddPlasma (int lev, RealBox part_realbox) if(loc_has_breit_wheeler){ breit_wheeler_get_opt = shr_ptr_bw_engine->build_optical_depth_functor(); - } + } #endif const GpuArray overlap_corner -- cgit v1.2.3 From 3e3670639e3f57cdc9f1a43a6c69d7e8308f100d Mon Sep 17 00:00:00 2001 From: lucafedeli88 Date: Fri, 11 Oct 2019 19:14:03 -0500 Subject: Removed merge conflict from file --- Source/Particles/PhysicalParticleContainer.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 7f6b2daa9..c3ea3b763 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -2205,14 +2205,13 @@ PhysicalParticleContainer::buildIonizationMask (const amrex::MFIter& mfi, const ); } -<<<<<<< HEAD //This function return true if the PhysicalParticleContainer contains electrons //or positrons, false otherwise bool PhysicalParticleContainer::AmIALepton(){ return (this-> mass == PhysConst::m_e); } -======= + #ifdef WARPX_QED bool PhysicalParticleContainer::has_quantum_sync() @@ -2239,4 +2238,3 @@ set_quantum_sync_engine_ptr(std::shared_ptr ptr) shr_ptr_qs_engine = ptr; } #endif ->>>>>>> upstream/dev -- cgit v1.2.3