From 38b802c0ba0d47fd179a2b384e8867248de078e9 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 7 May 2019 08:35:22 -0700 Subject: consistent continuous injection for PhysicalParticle and LaserParticle --- Source/Particles/PhysicalParticleContainer.cpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 17e6d98d9..c28aacd28 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -24,7 +24,8 @@ NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + // if (injected) { + if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; Real y = overlap_corner[1] + (iv[1] + 0.5)*dx[1]; @@ -81,6 +82,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp pp.query("do_backward_propagation", do_backward_propagation); pp.query("do_splitting", do_splitting); pp.query("split_type", split_type); + pp.query("do_continuous_injection", do_continuous_injection); } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) @@ -361,7 +363,9 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - if (injected) { + // if (injected) { + if (do_continuous_injection) { + Print()<<"in AddPlasmaCPU if do_continuous_injection"< Date: Tue, 7 May 2019 11:48:48 -0700 Subject: continuous injection working + comments --- Source/Evolve/WarpXEvolveEM.cpp | 4 - Source/Laser/LaserParticleContainer.H | 8 +- Source/Laser/LaserParticleContainer.cpp | 134 ++++++++++--------------- Source/Particles/MultiParticleContainer.H | 8 +- Source/Particles/MultiParticleContainer.cpp | 32 ++---- Source/Particles/PhysicalParticleContainer.H | 6 +- Source/Particles/PhysicalParticleContainer.cpp | 13 ++- Source/Particles/WarpXParticleContainer.H | 7 +- Source/Utils/WarpXMovingWindow.cpp | 29 +++--- Source/WarpX.cpp | 12 --- 10 files changed, 98 insertions(+), 155 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 7fe47d261..e98561be1 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -88,10 +88,6 @@ WarpX::EvolveEM (int numsteps) UpdateAuxilaryData(); } - // Performs continuous injection of all WarpXParticleContainer - // in mypc. - // mypc->ContinuousInjection(dt[0], geom[0].ProbDomain()); - if (do_subcycling == 0 || finest_level == 0) { OneStep_nosub(cur_time); } else if (do_subcycling == 1 && finest_level == 1) { diff --git a/Source/Laser/LaserParticleContainer.H b/Source/Laser/LaserParticleContainer.H index bb36f4795..ccd48f3b5 100644 --- a/Source/Laser/LaserParticleContainer.H +++ b/Source/Laser/LaserParticleContainer.H @@ -84,10 +84,11 @@ private: std::string field_function; // laser particle domain - amrex::RealBox laser_prob_domain; + amrex::RealBox laser_injection_box; // Theoretical position of the antenna. Used if do_continuous_injection=1. // Track the position of the antenna until it enters the simulation domain. - amrex::Real z_antenna_th; + amrex::Vector updated_position; + // Avoid injecting the antenna several times. int done_injecting=1; void ComputeSpacing (int lev, amrex::Real& Sx, amrex::Real& Sy) const; @@ -95,7 +96,8 @@ private: void InitData (int lev); // Inject the laser antenna during the simulation, if it started // outside of the simulation domain and enters it. - void ContinuousInjection(amrex::Real dt, const amrex::RealBox& prob_domain) override; + void ContinuousInjection(const amrex::RealBox& injection_box) override; + // Update position of the antenna void UpdateContinuousInjectionPosition(amrex::Real dt) override; }; diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 8499455de..a7fe9042d 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -149,42 +149,42 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, u_Y = {0., 1., 0.}; #endif - laser_prob_domain = Geometry::ProbDomain(); + laser_injection_box= Geometry::ProbDomain(); { Vector lo, hi; if (pp.queryarr("prob_lo", lo)) { - laser_prob_domain.setLo(lo); + laser_injection_box.setLo(lo); } if (pp.queryarr("prob_hi", hi)) { - laser_prob_domain.setHi(hi); + laser_injection_box.setHi(hi); } } if (do_continuous_injection){ // If laser antenna initially outside of the box, store its theoretical // position in z_antenna_th, and set done_injecting to 0. - z_antenna_th = position[2]; - const Real prob_lo_z = laser_prob_domain.lo()[AMREX_SPACEDIM-1]; - const Real prob_hi_z = laser_prob_domain.hi()[AMREX_SPACEDIM-1]; - if ( z_antenna_thprob_hi_z ){ + updated_position = position; + // Convert updated position to Real* to use RealBox.contains() +#if (AMREX_SPACEDIM == 3) + const Real* p_pos = updated_position.dataPtr(); +#else + const Real p_pos[2] = {updated_position[0], updated_position[2]}; +#endif + if ( not laser_injection_box.contains(p_pos) ){ done_injecting = 0; } - // Sanity checks: do_continuous_injection can be used only if the - // laser, the moving window and the boost are all in z direction. + + // Sanity checks + std::Vector windir(3, 0.0); +#if (AMREX_SPACEDIM==2) + windir[2*dir] = 1.0; +#else + windir[dir] = 1.0; +#endif AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - (nvec[0]-0)*(nvec[0]-0) + - (nvec[1]-0)*(nvec[1]-0) + - (nvec[2]-1)*(nvec[2]-1) < 1.e-12, - "do_continous_injection for laser particle only works if " + - "laser.direction = 0 0 1 (laser along z). TODO: all directions."); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - WarpX::moving_window_dir == AMREX_SPACEDIM-1, - "do_continous_injection for laser particle only works if " + - "moving window along z. TODO: all directions."); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE( - maxLevel() == 0, - "do_continous_injection for laser particle only works if " + - "max level = 0."); + (nvec[0]-windir[0]) + (nvec[1]-windor[1]) + (nvec[2]-windor[2]) + < 1.e-12, "do_continous_injection for laser particle only works" + + " if moving window direction and laser propagation direction are the same"); if ( WarpX::gamma_boost>1 ){ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( (WarpX::boost_direction[0]-0)*(WarpX::boost_direction[0]-0) + @@ -197,86 +197,60 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, } /* \brief Check if laser particles enter the box, and inject if necessary. - * \param dt: time step (assume no MR for now) - * \param prob_domain: a RealBox that contains current simulation boundaries. - * This function checks if the laser antenna should be injected at this - * iteration. If so, it injects it and set done_injecting to 1. + * \param injection_box: a RealBox where particles should be injected. */ void -LaserParticleContainer::ContinuousInjection (Real dt, const RealBox& prob_domain) +LaserParticleContainer::ContinuousInjection (const RealBox& injection_box) { - Print()<<" --- In LaserParticleContainer::ContinuousInjection"<=prob_lo_z && z_antenna_th 1)){ - // In boosted-frame simulations, the plasma has moved since the last + // In boosted-frame simulations, the antenna has moved since the last // call to this function, and injection position needs to be updated - z_antenna_th -= WarpX::beta_boost * #if ( AMREX_SPACEDIM == 3 ) + updated_position[dir] -= WarpX::beta_boost * WarpX::boost_direction[dir] * PhysConst::c * dt; #elif ( AMREX_SPACEDIM == 2 ) - // In 2D, dir=0 corresponds to x and dir=1 corresponds to z - // This needs to be converted in order to index `boost_direction` - // which has 3 components, for both 2D and 3D simulations. + // In 2D, dir=0 corresponds to x and dir=1 corresponds to z + // This needs to be converted in order to index `boost_direction` + // which has 3 components, for both 2D and 3D simulations. + updated_position[2*dir] -= WarpX::beta_boost * WarpX::boost_direction[2*dir] * PhysConst::c * dt; #endif } - - /* - if (WarpX::do_plasma_injection and (WarpX::gamma_boost > 1){ - z_antenna_th -= PhysConst::c * WarpX::beta_boost * dt; - } - */ } - - - - - - - - - - - - - - - void LaserParticleContainer::InitData () { @@ -294,9 +268,9 @@ LaserParticleContainer::InitData (int lev) // LaserParticleContainer::position contains the initial position of the // laser antenna. In the boosted frame, the antenna is moving. - // Update its position. + // Update its position with updated_position. if (do_continuous_injection){ - position[2] = z_antenna_th; + position = updated_position; } auto Transform = [&](int i, int j) -> Vector{ @@ -334,8 +308,8 @@ LaserParticleContainer::InitData (int lev) plane_hi[1] = std::max(plane_hi[1], j); }; - const Real* prob_lo = laser_prob_domain.lo(); - const Real* prob_hi = laser_prob_domain.hi(); + const Real* prob_lo = laser_injection_box.lo(); + const Real* prob_hi = laser_injection_box.hi(); #if (AMREX_SPACEDIM == 3) compute_min_max(prob_lo[0], prob_lo[1], prob_lo[2]); compute_min_max(prob_hi[0], prob_lo[1], prob_lo[2]); @@ -396,7 +370,7 @@ LaserParticleContainer::InitData (int lev) #else const Real x[2] = {pos[0], pos[2]}; #endif - if (laser_prob_domain.contains(x)) + if (laser_injection_box.contains(x)) { for (int k = 0; k<2; ++k) { particle_x.push_back(pos[0]); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index f48daf9bb..49a79231c 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -170,9 +170,11 @@ public: const amrex::Real t_boost, const amrex::Real t_lab, const amrex::Real dt, amrex::Vector& parts) const; - void ContinuousInjection(amrex::Real dt, - const amrex::RealBox& prob_domain) const; - + // Inject particles during the simulation (for particles entering the + // simulation domain after some iterations, due to flowing plasma and/or + // moving window). + void ContinuousInjection(const amrex::RealBox& injection_box) const; + // Update injection position for continuously-injected species. void UpdateContinuousInjectionPosition(amrex::Real dt) const; // diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 9810a168f..a5f695c69 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -416,45 +416,33 @@ MultiParticleContainer } /* \brief Continuous injection for particles initially outside of the domain. - * \param dt: timestep (so far, this only works without MR) - * \param prob_domain: current boundaries of the full domain. + * \param injection_box: Domain where new particles should be injected. * Loop over all WarpXParticleContainer in MultiParticleContainer and * calls virtual function ContinuousInjection. */ void -MultiParticleContainer::ContinuousInjection(Real dt, const RealBox& prob_domain) const +MultiParticleContainer::ContinuousInjection(const RealBox& injection_box) const { for (int i=0; ido_continuous_injection "<do_continuous_injection<do_continuous_injection) - { - Print()<<"i "<do_continuous_injection "<do_continuous_injection<ContinuousInjection(dt, prob_domain); + if (pc->do_continuous_injection){ + pc->ContinuousInjection(injection_box); } } } +/* \brief Update position of continuous injection parameters. + * \param dt: simulation time step (level 0) + * All classes inherited from WarpXParticleContainer do not have + * a position to update (PhysicalParticleContainer does not do anything). + */ void MultiParticleContainer::UpdateContinuousInjectionPosition(Real dt) const { for (int i=0; ido_continuous_injection "<do_continuous_injection<do_continuous_injection) - { - Print()<<"i "<do_continuous_injection "<do_continuous_injection<do_continuous_injection){ pc->UpdateContinuousInjectionPosition(dt); } } - /* - for (int i=nspecies; i(pc); - // auto& pc = allcontainers[i]; - // auto& lpc = dynamic_cast(pc); - lpc.UpdateContinuousInjectionPosition(dt); - } - */ } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 253aa58cd..4f365768b 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -101,9 +101,6 @@ public: const amrex::Real t_lab, const amrex::Real dt, DiagnosticParticles& diagnostic_particles) final; - // bool injected = false; - // int do_continuous_injection = false; - protected: std::string species_name; @@ -123,7 +120,8 @@ protected: int GetRefineFac(const amrex::Real x, const amrex::Real y, const amrex::Real z); std::unique_ptr m_refined_injection_mask = nullptr; - void ContinuousInjection(amrex::Real dt, const amrex::RealBox& prob_domain) override; + // Inject particles during the whole simulation + void ContinuousInjection(const amrex::RealBox& injection_box) override; }; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c28aacd28..2fa39d87d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -24,7 +24,6 @@ NumParticlesToAdd(const Box& overlap_box, const RealBox& overlap_realbox, for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - // if (injected) { if (do_continuous_injection) { #if ( AMREX_SPACEDIM == 3 ) Real x = overlap_corner[0] + (iv[0] + 0.5)*dx[0]; @@ -363,9 +362,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) for (IntVect iv = overlap_box.smallEnd(); iv <= overlap_box.bigEnd(); overlap_box.next(iv)) { int fac; - // if (injected) { if (do_continuous_injection) { - Print()<<"in AddPlasmaCPU if do_continuous_injection"< local_rho; amrex::Vector local_jx; diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 2ba1d2f59..38bde5e92 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -34,9 +34,16 @@ WarpX::MoveWindow (bool move_j) moving_window_x += moving_window_v * dt[0]; int dir = moving_window_dir; + // Update warpx.current_injection_position + // PhysicalParticleContainer uses this injection position UpdatePlasmaInjectionPosition( dt[0] ); - mypc->UpdateContinuousInjectionPosition( dt[0] ); - + if (WarpX::do_plasma_injection){ + // Update injection position for WarpXParticleContainer in mypc. + // Nothing to do for PhysicalParticleContainers + // For LaserParticleContainer, need to update the antenna position. + mypc->UpdateContinuousInjectionPosition( dt[0] ); + } + // compute the number of cells to shift on the base level Real new_lo[AMREX_SPACEDIM]; Real new_hi[AMREX_SPACEDIM]; @@ -164,21 +171,11 @@ WarpX::MoveWindow (bool move_j) particleBox.setLo( dir, new_injection_position ); particleBox.setHi( dir, current_injection_position ); } - // Perform the injection of new particles in particleBox - // Performs continuous injection of all WarpXParticleContainer - // in mypc. - + if (particleBox.ok() and (current_injection_position != new_injection_position)){ - mypc->ContinuousInjection(dt[0], particleBox); - /* - for (int i = 0; i < num_injected_species; ++i) { - int ispecies = injected_plasma_species[i]; - WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - auto& ppc = dynamic_cast(pc); - ppc.AddPlasma(lev, particleBox); - } - */ - // Update the injection position + // Performs continuous injection of all WarpXParticleContainer + // in mypc. + mypc->ContinuousInjection(particleBox); current_injection_position = new_injection_position; } } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 645e92161..bea91f1ba 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -146,18 +146,6 @@ WarpX::WarpX () // Particle Container mypc = std::unique_ptr (new MultiParticleContainer(this)); - /* - if (do_plasma_injection) { - for (int i = 0; i < num_injected_species; ++i) { - int ispecies = injected_plasma_species[i]; - WarpXParticleContainer& pc = mypc->GetParticleContainer(ispecies); - auto& ppc = dynamic_cast(pc); - // ppc.injected = true; - ppc.do_continuous_injection = 1; - } - } - */ - Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); -- cgit v1.2.3 From cae3e9ce6769fd5c626e27fad473da17a2b027fb Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 7 May 2019 17:47:56 -0700 Subject: NCI filter initialized. See how to apply it --- Source/Filter/NCIGodfreyFilter.H | 20 ++++----- Source/Filter/NCIGodfreyFilter.cpp | 59 +++++++++++++++++++++++++- Source/Initialization/WarpXInitData.cpp | 17 +++++--- Source/Particles/PhysicalParticleContainer.cpp | 40 +++++++++++++++++ Source/Utils/NCIGodfreyTables.H | 7 ++- 5 files changed, 123 insertions(+), 20 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H index 50907a241..2e6e65772 100644 --- a/Source/Filter/NCIGodfreyFilter.H +++ b/Source/Filter/NCIGodfreyFilter.H @@ -1,5 +1,4 @@ #include -#include #ifndef WARPX_GODFREY_FILTER_H_ #define WARPX_GODFREY_FILTER_H_ @@ -10,16 +9,15 @@ class NCIGodfreyFilter : public Filter { public: - NCIGodfreyFilter(godfrey_coeff_set coeff_set_in) - : Filter(), - coeff_set(coeff_set_in) + NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v); + /* { - if (coeff_set == godfrey_coeff_set::Ex_Ey_Bz){ - coeff_table = (amrex::Real*) table_nci_godfrey_Ex_Ey_Bz; - } else if (coeff_set == godfrey_coeff_set::Bx_By_Ez){ - coeff_table = (amrex::Real*) table_nci_godfrey_Bx_By_Ez; - } + : coeff_set( coeff_set_ ), + cdtodz( cdtodz_ ), + l_lower_order_in_v( l_lower_order_in_v_ ), + slen }; + */ virtual void ComputeStencils() override final; @@ -28,7 +26,9 @@ public: private: godfrey_coeff_set coeff_set; - amrex::Real* coeff_table; + amrex::Real cdtodz; + int l_lower_order_in_v; + // amrex::Real* coeff_table; }; diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 652bf9d85..e58cc5ac8 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -1,5 +1,6 @@ #include #include +#include #ifdef _OPENMP #include @@ -7,6 +8,62 @@ using namespace amrex; +NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_){ + coeff_set = coeff_set_; + cdtodz = cdtodz_; + l_lower_order_in_v = l_lower_order_in_v_; +#if (AMREX_SPACEDIM == 3) + slen = {1,1,5}; +#else + slen = {1,5,1}; +#endif +} + void NCIGodfreyFilter::ComputeStencils(){ - int i=0; + Print()<<"slen "<fdtd_nci_stencilz_by.resize(max_level+1); for (int lev = 0; lev <= max_level; ++lev) { - nci_godfrey_filter_exeybz[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Ex_Ey_Bz) ); - nci_godfrey_filter_bxbyez[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Bx_By_Ez) ); - const Geometry& gm = Geom(lev); const Real* dx = gm.CellSize(); amrex::Real dz, cdtodz; @@ -179,10 +176,16 @@ WarpX::InitNCICorrector () } cdtodz = PhysConst::c * dt[lev] / dz; - WRPX_PXR_NCI_CORR_INIT( (mypc->fdtd_nci_stencilz_ex)[lev].data(), - (mypc->fdtd_nci_stencilz_by)[lev].data(), - mypc->nstencilz_fdtd_nci_corr, cdtodz, - WarpX::l_lower_order_in_v); + nci_godfrey_filter_exeybz[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Ex_Ey_Bz, cdtodz, WarpX::l_lower_order_in_v) ); + nci_godfrey_filter_bxbyez[lev].reset( new NCIGodfreyFilter(godfrey_coeff_set::Bx_By_Ez, cdtodz, WarpX::l_lower_order_in_v) ); + + nci_godfrey_filter_exeybz[lev]->ComputeStencils(); + nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); + + //WRPX_PXR_NCI_CORR_INIT( (mypc->fdtd_nci_stencilz_ex)[lev].data(), + // (mypc->fdtd_nci_stencilz_by)[lev].data(), + // mypc->nstencilz_fdtd_nci_corr, cdtodz, + // WarpX::l_lower_order_in_v); } } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 17e6d98d9..4a60fc630 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1104,6 +1104,26 @@ PhysicalParticleContainer::Evolve (int lev, const auto& mypc = WarpX::GetInstance().GetPartContainer(); const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; + + + + + + + + + const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; + const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; + + + + + + + + + + BL_ASSERT(OnSameGrids(lev,jx)); MultiFab* cost = WarpX::getCosts(lev); @@ -1171,53 +1191,73 @@ PhysicalParticleContainer::Evolve (int lev, #endif // both 2d and 3d + // bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); + filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), BL_TO_FORTRAN_ANYD(filtered_Ex), BL_TO_FORTRAN_ANYD(Ex[pti]), mypc.fdtd_nci_stencilz_ex[lev].data(), &nstencilz_fdtd_nci_corr); + */ exfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), BL_TO_FORTRAN_ANYD(filtered_Ez), BL_TO_FORTRAN_ANYD(Ez[pti]), mypc.fdtd_nci_stencilz_by[lev].data(), &nstencilz_fdtd_nci_corr); + */ ezfab = &filtered_Ez; filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), BL_TO_FORTRAN_ANYD(filtered_By), BL_TO_FORTRAN_ANYD(By[pti]), mypc.fdtd_nci_stencilz_by[lev].data(), &nstencilz_fdtd_nci_corr); + */ byfab = &filtered_By; #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), BL_TO_FORTRAN_ANYD(filtered_Ey), BL_TO_FORTRAN_ANYD(Ey[pti]), mypc.fdtd_nci_stencilz_ex[lev].data(), &nstencilz_fdtd_nci_corr); + */ eyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), BL_TO_FORTRAN_ANYD(filtered_Bx), BL_TO_FORTRAN_ANYD(Bx[pti]), mypc.fdtd_nci_stencilz_by[lev].data(), &nstencilz_fdtd_nci_corr); + */ bxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), BL_TO_FORTRAN_ANYD(filtered_Bz), BL_TO_FORTRAN_ANYD(Bz[pti]), mypc.fdtd_nci_stencilz_ex[lev].data(), &nstencilz_fdtd_nci_corr); + */ bzfab = &filtered_Bz; #endif } diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H index dc712f95e..c181aead9 100644 --- a/Source/Utils/NCIGodfreyTables.H +++ b/Source/Utils/NCIGodfreyTables.H @@ -3,7 +3,10 @@ #ifndef WARPX_GODFREY_COEFF_TABLE_H_ #define WARPX_GODFREY_COEFF_TABLE_H_ -const amrex::Real table_nci_godfrey_Ex_Ey_Bz[4][100]{ +const int tab_width = 4; +const int tab_length = 100; + +const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ -2.47536,2.04288,-0.598163,0.0314711, -2.47545,2.04309,-0.598307,0.0315029, -2.4756,2.04342,-0.598549,0.0315558, @@ -106,7 +109,7 @@ const amrex::Real table_nci_godfrey_Ex_Ey_Bz[4][100]{ -2.98769,3.19374,-1.41463,0.208607 }; -const amrex::Real table_nci_godfrey_Bx_By_Ez[4][100]{ +const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ -2.80862,2.80104,-1.14615,0.154077, -2.80851,2.80078,-1.14595,0.154027, -2.80832,2.80034,-1.14561,0.153945, -- cgit v1.2.3 From bd9e80c56baf4bb6d8760befb3a42896f03fe912 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 11:04:43 -0700 Subject: Apply new C++ NCI filter. OK in 2D --- Source/Filter/Filter.H | 5 ++ Source/Filter/Filter.cpp | 69 +++++++++++++++++++++++++- Source/Filter/NCIGodfreyFilter.cpp | 2 + Source/Particles/PhysicalParticleContainer.cpp | 50 ++++++++++++++++--- 4 files changed, 119 insertions(+), 7 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index ef322d8e7..86218c913 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -14,6 +14,11 @@ public: const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); + void ApplyStencil (amrex::FArrayBox& dstfab, + const amrex::FArrayBox& srcfab, const amrex::Box& tbx, + const amrex::Box& srcbx, + int scomp=0, int dcomp=0, int ncomp=10000); + // public for cuda void DoFilter(const amrex::Box& tbx, amrex::Array4 const& tmp, diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 7e13f1cc9..ab1a6f18e 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -12,7 +12,7 @@ using namespace amrex; void Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) { - BL_PROFILE("BilinearFilter::ApplyStencil()"); + BL_PROFILE("BilinearFilter::ApplyStencil(MultiFab)"); ncomp = std::min(ncomp, srcmf.nComp()); for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) @@ -43,6 +43,43 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, + const Box& srcbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + + // for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) + //{ + const auto& src = srcfab.array(mfi); + const auto& dst = dstfab.array(mfi); + // const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); + + // Copy values in srcfab into tmpfab + // const Box& ibx = gbx & srcmf[mfi].box(); + const Box& ibx = gbx & srcbx; + + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + { + if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { + tmp(i,j,k,n) = src(i,j,k,n+scomp); + } else { + tmp(i,j,k,n) = 0.0; + } + }); + + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); + //} +} + void Filter::DoFilter (const Box& tbx, Array4 const& tmp, Array4 const& dst, @@ -116,6 +153,36 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +void +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, const Box& srcbx, int scomp, int dcomp, int ncomp) +{ + BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); + ncomp = std::min(ncomp, srcfab.nComp()); + //#ifdef _OPENMP + //#pragma omp parallel + //#endif + //{ + FArrayBox tmpfab; + //for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ + //const auto& srcfab = srcmf[mfi]; + // auto& dstfab = dstmf[mfi]; + //const auto& srcfab = srcmf[mfi]; + // auto& dstfab = dstmf[mfi]; + // const Box& tbx = mfi.growntilebox(); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + //const Box& ibx = gbx & srcfab.box(); + const Box& ibx = gbx & srcbx; + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); + //} + //} +} + void Filter::DoFilter (const Box& tbx, Array4 const& tmp, Array4 const& dst, diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index e58cc5ac8..095cb1d8b 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -13,8 +13,10 @@ NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdt cdtodz = cdtodz_; l_lower_order_in_v = l_lower_order_in_v_; #if (AMREX_SPACEDIM == 3) + stencil_length_each_dir = {1,1,5}; slen = {1,1,5}; #else + stencil_length_each_dir = {1,5}; slen = {1,5,1}; #endif } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 4a60fc630..771a41a67 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1179,6 +1179,7 @@ PhysicalParticleContainer::Evolve (int lev, FArrayBox const* byfab = &(By[pti]); FArrayBox const* bzfab = &(Bz[pti]); + Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1194,7 +1195,9 @@ PhysicalParticleContainer::Evolve (int lev, // bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) + exeli = filtered_Ex.elixir(); + //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), BL_TO_FORTRAN_ANYD(filtered_Ex), @@ -1205,7 +1208,9 @@ PhysicalParticleContainer::Evolve (int lev, exfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) + ezeli = filtered_Ez.elixir(); + //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), BL_TO_FORTRAN_ANYD(filtered_Ez), @@ -1216,7 +1221,9 @@ PhysicalParticleContainer::Evolve (int lev, ezfab = &filtered_Ez; filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) + byeli = filtered_By.elixir(); + //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), BL_TO_FORTRAN_ANYD(filtered_By), @@ -1228,7 +1235,9 @@ PhysicalParticleContainer::Evolve (int lev, #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) + eyeli = filtered_Ey.elixir(); + //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), BL_TO_FORTRAN_ANYD(filtered_Ey), @@ -1239,7 +1248,9 @@ PhysicalParticleContainer::Evolve (int lev, eyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) + bxeli = filtered_Bx.elixir(); + //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), BL_TO_FORTRAN_ANYD(filtered_Bx), @@ -1250,7 +1261,9 @@ PhysicalParticleContainer::Evolve (int lev, bxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) + bzeli = filtered_Bz.elixir(); + //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, box); /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), BL_TO_FORTRAN_ANYD(filtered_Bz), @@ -1422,6 +1435,7 @@ PhysicalParticleContainer::Evolve (int lev, const FArrayBox* cbyfab = &(*cBy)[pti]; const FArrayBox* cbzfab = &(*cBz)[pti]; + Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) @@ -1435,51 +1449,75 @@ PhysicalParticleContainer::Evolve (int lev, // both 2d and 3d filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + exeli = filtered_Ex.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), BL_TO_FORTRAN_ANYD(filtered_Ex), BL_TO_FORTRAN_ANYD((*cEx)[pti]), mypc.fdtd_nci_stencilz_ex[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cexfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); + ezeli = filtered_Ez.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), BL_TO_FORTRAN_ANYD(filtered_Ez), BL_TO_FORTRAN_ANYD((*cEz)[pti]), mypc.fdtd_nci_stencilz_by[lev-1].data(), &nstencilz_fdtd_nci_corr); cezfab = &filtered_Ez; + */ filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); + byeli = filtered_By.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), BL_TO_FORTRAN_ANYD(filtered_By), BL_TO_FORTRAN_ANYD((*cBy)[pti]), mypc.fdtd_nci_stencilz_by[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cbyfab = &filtered_By; #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); + eyeli = filtered_Ey.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), BL_TO_FORTRAN_ANYD(filtered_Ey), BL_TO_FORTRAN_ANYD((*cEy)[pti]), mypc.fdtd_nci_stencilz_ex[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ ceyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); + bxeli = filtered_Bx.elixir(); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), BL_TO_FORTRAN_ANYD(filtered_Bx), BL_TO_FORTRAN_ANYD((*cBx)[pti]), mypc.fdtd_nci_stencilz_by[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cbxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); + bzeli = filtered_Bz.elixir(); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, cbox); + /* WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), BL_TO_FORTRAN_ANYD(filtered_Bz), BL_TO_FORTRAN_ANYD((*cBz)[pti]), mypc.fdtd_nci_stencilz_ex[lev-1].data(), &nstencilz_fdtd_nci_corr); + */ cbzfab = &filtered_Bz; #endif } -- cgit v1.2.3 From e41e5eac1bccaa5c3732154aa052aee6c32f6518 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 15:24:22 -0700 Subject: fixes, and add calls to the filter --- Source/Filter/Filter.H | 4 +- Source/Filter/Filter.cpp | 70 +++++------- Source/Filter/NCIGodfreyFilter.H | 11 +- Source/Filter/NCIGodfreyFilter.cpp | 11 +- Source/Particles/PhysicalParticleContainer.cpp | 141 +++---------------------- Source/Utils/NCIGodfreyTables.H | 4 +- Source/WarpX.H | 2 +- 7 files changed, 57 insertions(+), 186 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index 86218c913..d728a16e6 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -10,13 +10,13 @@ public: Filter () = default; virtual void ComputeStencils() = 0; + void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); void ApplyStencil (amrex::FArrayBox& dstfab, const amrex::FArrayBox& srcfab, const amrex::Box& tbx, - const amrex::Box& srcbx, int scomp=0, int dcomp=0, int ncomp=10000); // public for cuda @@ -28,8 +28,8 @@ public: amrex::IntVect stencil_length_each_dir; protected: - amrex::Dim3 slen; amrex::Gpu::ManagedVector stencil_x, stencil_y, stencil_z; + amrex::Dim3 slen; private: diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index ab1a6f18e..9bbd230a1 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -44,29 +44,22 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, - const Box& srcbx, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.nComp()); + const auto& src = srcfab.array(mfi); + const auto& dst = dstfab.array(mfi); + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - // for (MFIter mfi(dstmf); mfi.isValid(); ++mfi) - //{ - const auto& src = srcfab.array(mfi); - const auto& dst = dstfab.array(mfi); - // const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - - // tmpfab has enough ghost cells for the stencil - FArrayBox tmp_fab(gbx,ncomp); - Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early - auto const& tmp = tmp_fab.array(); + // tmpfab has enough ghost cells for the stencil + FArrayBox tmp_fab(gbx,ncomp); + Elixir tmp_eli = tmp_fab.elixir(); // Prevent the tmp data from being deleted too early + auto const& tmp = tmp_fab.array(); - // Copy values in srcfab into tmpfab - // const Box& ibx = gbx & srcmf[mfi].box(); - const Box& ibx = gbx & srcbx; - - AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, + // Copy values in srcfab into tmpfab + AMREX_PARALLEL_FOR_4D ( gbx, ncomp, i, j, k, n, { if (ibx.contains(IntVect(AMREX_D_DECL(i,j,k)))) { tmp(i,j,k,n) = src(i,j,k,n+scomp); @@ -75,9 +68,8 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx } }); - // Apply filter - DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); - //} + // Apply filter + DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); } void Filter::DoFilter (const Box& tbx, @@ -154,33 +146,21 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } void -Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, const Box& srcbx, int scomp, int dcomp, int ncomp) +Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, + const Box& tbx, int scomp, int dcomp, int ncomp) { BL_PROFILE("BilinearFilter::ApplyStencil(FArrayBox)"); ncomp = std::min(ncomp, srcfab.nComp()); - //#ifdef _OPENMP - //#pragma omp parallel - //#endif - //{ - FArrayBox tmpfab; - //for (MFIter mfi(dstmf,true); mfi.isValid(); ++mfi){ - //const auto& srcfab = srcmf[mfi]; - // auto& dstfab = dstmf[mfi]; - //const auto& srcfab = srcmf[mfi]; - // auto& dstfab = dstmf[mfi]; - // const Box& tbx = mfi.growntilebox(); - const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); - // tmpfab has enough ghost cells for the stencil - tmpfab.resize(gbx,ncomp); - tmpfab.setVal(0.0, gbx, 0, ncomp); - // Copy values in srcfab into tmpfab - //const Box& ibx = gbx & srcfab.box(); - const Box& ibx = gbx & srcbx; - tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); - // Apply filter - DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); - //} - //} + FArrayBox tmpfab; + const Box& gbx = amrex::grow(tbx,stencil_length_each_dir-1); + // tmpfab has enough ghost cells for the stencil + tmpfab.resize(gbx,ncomp); + tmpfab.setVal(0.0, gbx, 0, ncomp); + // Copy values in srcfab into tmpfab + const Box& ibx = gbx; + tmpfab.copy(srcfab, ibx, scomp, ibx, 0, ncomp); + // Apply filter + DoFilter(tbx, tmpfab.array(), dstfab.array(), 0, dcomp, ncomp); } void Filter::DoFilter (const Box& tbx, diff --git a/Source/Filter/NCIGodfreyFilter.H b/Source/Filter/NCIGodfreyFilter.H index 2e6e65772..21e16ff74 100644 --- a/Source/Filter/NCIGodfreyFilter.H +++ b/Source/Filter/NCIGodfreyFilter.H @@ -9,15 +9,9 @@ class NCIGodfreyFilter : public Filter { public: + NCIGodfreyFilter () = default; + NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v); - /* - { - : coeff_set( coeff_set_ ), - cdtodz( cdtodz_ ), - l_lower_order_in_v( l_lower_order_in_v_ ), - slen - }; - */ virtual void ComputeStencils() override final; @@ -28,7 +22,6 @@ private: godfrey_coeff_set coeff_set; amrex::Real cdtodz; int l_lower_order_in_v; - // amrex::Real* coeff_table; }; diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 095cb1d8b..bab1cc33a 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -26,6 +26,7 @@ void NCIGodfreyFilter::ComputeStencils(){ Print()<<"slen.x "<(WarpX::noz)}); #endif - // both 2d and 3d - // bilinear_filter.ApplyStencil(*j_fp[lev][idim], *current_fp[lev][idim]); - filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); exeli = filtered_Ex.elixir(); - //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ex[pti], filtered_Ex) - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD(Ex[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox); exfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); - //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Ez[pti], filtered_Ez) - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD(Ez[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox); ezfab = &filtered_Ez; filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); - //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(By[pti], filtered_By) - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD(By[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox); byfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); - //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Ey[pti], filtered_Ey) - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD(Ey[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox); eyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); - //nci_godfrey_filter_bxbyez[lev]->ApplyStencil(Bx[pti], filtered_Bx) - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD(Bx[pti]), - mypc.fdtd_nci_stencilz_by[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox); bxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); - //nci_godfrey_filter_exeybz[lev]->ApplyStencil(Bz[pti], filtered_Bz) - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, box); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD(Bz[pti]), - mypc.fdtd_nci_stencilz_ex[lev].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox); bzfab = &filtered_Bz; #endif } @@ -1451,74 +1382,32 @@ PhysicalParticleContainer::Evolve (int lev, // both 2d and 3d filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); exeli = filtered_Ex.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ex), - BL_TO_FORTRAN_ANYD(filtered_Ex), - BL_TO_FORTRAN_ANYD((*cEx)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], tbox); cexfab = &filtered_Ex; filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ez), - BL_TO_FORTRAN_ANYD(filtered_Ez), - BL_TO_FORTRAN_ANYD((*cEz)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], tbox); cezfab = &filtered_Ez; - */ + filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_By), - BL_TO_FORTRAN_ANYD(filtered_By), - BL_TO_FORTRAN_ANYD((*cBy)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], tbox); cbyfab = &filtered_By; - #if (AMREX_SPACEDIM == 3) filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Ey), - BL_TO_FORTRAN_ANYD(filtered_Ey), - BL_TO_FORTRAN_ANYD((*cEy)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], tbox); ceyfab = &filtered_Ey; filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bx), - BL_TO_FORTRAN_ANYD(filtered_Bx), - BL_TO_FORTRAN_ANYD((*cBx)[pti]), - mypc.fdtd_nci_stencilz_by[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Bx, (*cBx)[pti], tbox); cbxfab = &filtered_Bx; filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox, cbox); - /* - WRPX_PXR_GODFREY_FILTER(BL_TO_FORTRAN_BOX(filtered_Bz), - BL_TO_FORTRAN_ANYD(filtered_Bz), - BL_TO_FORTRAN_ANYD((*cBz)[pti]), - mypc.fdtd_nci_stencilz_ex[lev-1].data(), - &nstencilz_fdtd_nci_corr); - */ + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], tbox); cbzfab = &filtered_Bz; #endif } diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H index c181aead9..b84f2cedb 100644 --- a/Source/Utils/NCIGodfreyTables.H +++ b/Source/Utils/NCIGodfreyTables.H @@ -4,9 +4,10 @@ #define WARPX_GODFREY_COEFF_TABLE_H_ const int tab_width = 4; -const int tab_length = 100; +const int tab_length = 101; const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ + -2.47536,2.04288,-0.598163,0.0314711, -2.47536,2.04288,-0.598163,0.0314711, -2.47545,2.04309,-0.598307,0.0315029, -2.4756,2.04342,-0.598549,0.0315558, @@ -110,6 +111,7 @@ const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ }; const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ + -2.80862,2.80104,-1.14615,0.154077, -2.80862,2.80104,-1.14615,0.154077, -2.80851,2.80078,-1.14595,0.154027, -2.80832,2.80034,-1.14561,0.153945, diff --git a/Source/WarpX.H b/Source/WarpX.H index 77bc47a7d..9ac64c6c6 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -146,7 +146,7 @@ public: BilinearFilter bilinear_filter; amrex::Vector< std::unique_ptr > nci_godfrey_filter_exeybz; amrex::Vector< std::unique_ptr > nci_godfrey_filter_bxbyez; - + static int num_mirrors; amrex::Vector mirror_z; amrex::Vector mirror_z_width; -- cgit v1.2.3 From 9d25c2b7c4aa473aa088a80f636086ed249b3a7a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 16:14:49 -0700 Subject: Add comments --- .../laser_acceleration/inputs.2d.boost | 14 +++++----- Source/Filter/Filter.H | 9 ++++++ Source/Filter/Filter.cpp | 32 ++++++++++++++++++++++ Source/Filter/NCIGodfreyFilter.cpp | 25 +++++++++++------ Source/Initialization/WarpXInitData.cpp | 3 +- Source/Particles/PhysicalParticleContainer.cpp | 22 +++++++++++++-- Source/Utils/NCIGodfreyTables.H | 7 +++++ Source/WarpX.cpp | 2 ++ 8 files changed, 96 insertions(+), 18 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost index b1fb28126..b8f0fff57 100644 --- a/Examples/Physics_applications/laser_acceleration/inputs.2d.boost +++ b/Examples/Physics_applications/laser_acceleration/inputs.2d.boost @@ -1,8 +1,8 @@ ################################# ######### BOX PARAMETERS ######## ################################# -# max_step = 2700 -stop_time = 1.9e-12 +max_step = 1000 +# stop_time = 1.9e-12 amr.n_cell = 128 1024 amr.max_grid_size = 64 amr.blocking_factor = 32 @@ -39,7 +39,7 @@ warpx.serialize_ics = 1 ################################# ####### BOOST PARAMETERS ######## ################################# -warpx.gamma_boost = 10. +warpx.gamma_boost = 30. warpx.boost_direction = z warpx.do_boosted_frame_diagnostic = 1 warpx.num_snapshots_lab = 7 @@ -61,11 +61,11 @@ electrons.momentum_distribution_type = "gaussian" electrons.xmin = -120.e-6 electrons.xmax = 120.e-6 electrons.zmin = 0.5e-3 -electrons.zmax = .0035 +electrons.zmax = 1. electrons.profile = "predefined" electrons.predefined_profile_name = "parabolic_channel" # predefined_profile_params = z_start ramp_up plateau ramp_down rc n0 -electrons.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 +electrons.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e25 electrons.do_continuous_injection = 1 ions.charge = q_e @@ -76,11 +76,11 @@ ions.momentum_distribution_type = "gaussian" ions.xmin = -120.e-6 ions.xmax = 120.e-6 ions.zmin = 0.5e-3 -ions.zmax = .0035 +ions.zmax = 1. ions.profile = "predefined" ions.predefined_profile_name = "parabolic_channel" # predefined_profile_params = z_start ramp_up plateau ramp_down rc n0 -ions.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e24 +ions.predefined_profile_params = .5e-3 .5e-3 2.e-3 .5e-3 50.e-6 3.5e25 ions.do_continuous_injection = 1 beam.charge = -q_e diff --git a/Source/Filter/Filter.H b/Source/Filter/Filter.H index d728a16e6..cf0848ee9 100644 --- a/Source/Filter/Filter.H +++ b/Source/Filter/Filter.H @@ -9,12 +9,16 @@ class Filter public: Filter () = default; + // This function has to be overriden by derived classes. virtual void ComputeStencils() = 0; + // Apply stencil on MultiFab. + // Guard cells are handled inside this function void ApplyStencil(amrex::MultiFab& dstmf, const amrex::MultiFab& srcmf, int scomp=0, int dcomp=0, int ncomp=10000); + // Apply stencil on a FabArray. void ApplyStencil (amrex::FArrayBox& dstfab, const amrex::FArrayBox& srcfab, const amrex::Box& tbx, int scomp=0, int dcomp=0, int ncomp=10000); @@ -25,10 +29,15 @@ public: amrex::Array4 const& dst, int scomp, int dcomp, int ncomp); + // In 2D, stencil_length_each_dir = {length(stencil_x), length(stencil_z)} amrex::IntVect stencil_length_each_dir; protected: + // Stencil along each direction. + // in 2D, stencil_y is not initialized. amrex::Gpu::ManagedVector stencil_x, stencil_y, stencil_z; + // Length of each stencil. + // In 2D, slen = {length(stencil_x), length(stencil_z), 1} amrex::Dim3 slen; private: diff --git a/Source/Filter/Filter.cpp b/Source/Filter/Filter.cpp index 9bbd230a1..25d935ef6 100644 --- a/Source/Filter/Filter.cpp +++ b/Source/Filter/Filter.cpp @@ -9,6 +9,13 @@ using namespace amrex; #ifdef AMREX_USE_CUDA +/* \brief Apply stencil on MultiFab (GPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) { @@ -43,6 +50,14 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +/* \brief Apply stencil on FArrayBox (GPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, int scomp, int dcomp, int ncomp) @@ -72,6 +87,8 @@ Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, DoFilter(tbx, tmp, dst, 0, dcomp, ncomp); } +/* \brief Apply stencil (2D/3D, CPU/GPU) + */ void Filter::DoFilter (const Box& tbx, Array4 const& tmp, Array4 const& dst, @@ -118,6 +135,13 @@ void Filter::DoFilter (const Box& tbx, #else +/* \brief Apply stencil on MultiFab (CPU version, 2D/3D). + * \param dstmf Destination MultiFab + * \param srcmf source MultiFab + * \param scomp first component of srcmf on which the filter is applied + * \param dcomp first component of dstmf on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dcomp, int ncomp) { @@ -145,6 +169,14 @@ Filter::ApplyStencil (MultiFab& dstmf, const MultiFab& srcmf, int scomp, int dco } } +/* \brief Apply stencil on FArrayBox (CPU version, 2D/3D). + * \param dstfab Destination FArrayBox + * \param srcmf source FArrayBox + * \param tbx Grown box on which srcfab is defined. + * \param scomp first component of srcfab on which the filter is applied + * \param dcomp first component of dstfab on which the filter is applied + * \param ncomp Number of components on which the filter is applied. + */ void Filter::ApplyStencil (FArrayBox& dstfab, const FArrayBox& srcfab, const Box& tbx, int scomp, int dcomp, int ncomp) diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index bab1cc33a..34fca7604 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -9,9 +9,11 @@ using namespace amrex; NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdtodz_, amrex::Real l_lower_order_in_v_){ + // Store parameters into class data members coeff_set = coeff_set_; cdtodz = cdtodz_; l_lower_order_in_v = l_lower_order_in_v_; + // NCI Godfrey filter has fixed size, and is applied along z only. #if (AMREX_SPACEDIM == 3) stencil_length_each_dir = {1,1,5}; slen = {1,1,5}; @@ -22,12 +24,6 @@ NCIGodfreyFilter::NCIGodfreyFilter(godfrey_coeff_set coeff_set_, amrex::Real cdt } void NCIGodfreyFilter::ComputeStencils(){ - Print()<<"slen "<ComputeStencils(); nci_godfrey_filter_bxbyez[lev]->ComputeStencils(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 0e003927d..1e0d68800 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1105,6 +1105,7 @@ PhysicalParticleContainer::Evolve (int lev, const auto& mypc = WarpX::GetInstance().GetPartContainer(); const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; + // Get instances of NCI Godfrey filters const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; @@ -1175,31 +1176,40 @@ PhysicalParticleContainer::Evolve (int lev, static_cast(WarpX::noz)}); #endif + // Filter Ex (Both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + // Safeguard for GPU exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox); + // Update exfab reference exfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox); ezfab = &filtered_Ez; + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox); byfab = &filtered_By; #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox); eyfab = &filtered_Ey; + // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox); bxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox); @@ -1378,33 +1388,41 @@ PhysicalParticleContainer::Evolve (int lev, static_cast(WarpX::noy), static_cast(WarpX::noz)}); #endif - - // both 2d and 3d + + // Filter Ex (both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); + // Safeguard for GPU exeli = filtered_Ex.elixir(); + // Apply filter on Ex, result stored in filtered_Ex nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], tbox); + // Update exfab reference cexfab = &filtered_Ex; + // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], tbox); cezfab = &filtered_Ez; + // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], tbox); cbyfab = &filtered_By; #if (AMREX_SPACEDIM == 3) + // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], tbox); 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], tbox); cbxfab = &filtered_Bx; + // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], tbox); diff --git a/Source/Utils/NCIGodfreyTables.H b/Source/Utils/NCIGodfreyTables.H index b84f2cedb..8cb105aa0 100644 --- a/Source/Utils/NCIGodfreyTables.H +++ b/Source/Utils/NCIGodfreyTables.H @@ -3,9 +3,14 @@ #ifndef WARPX_GODFREY_COEFF_TABLE_H_ #define WARPX_GODFREY_COEFF_TABLE_H_ +// Table width. This is related to the stencil length const int tab_width = 4; +// table length. Each line correspond to 1 value of cdtodz +// (here 101 values). const int tab_length = 101; +// Table of coefficient for Ex, Ey abd Bz +// We typically interpolate between two lines const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ -2.47536,2.04288,-0.598163,0.0314711, -2.47536,2.04288,-0.598163,0.0314711, @@ -110,6 +115,8 @@ const amrex::Real table_nci_godfrey_Ex_Ey_Bz[tab_length][tab_width]{ -2.98769,3.19374,-1.41463,0.208607 }; +// Table of coefficient for Bx, By and Ez +// We typically interpolate between two lines const amrex::Real table_nci_godfrey_Bx_By_Ez[tab_length][tab_width]{ -2.80862,2.80104,-1.14615,0.154077, -2.80862,2.80104,-1.14615,0.154077, diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 36dea14b1..5f63e2ad6 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -225,6 +225,8 @@ WarpX::WarpX () insitu_bridge = nullptr; #endif + // NCI Godfrey filters can have different stencils + // at different levels (the stencil depends on c*dt/dz) nci_godfrey_filter_exeybz.resize(nlevs_max); nci_godfrey_filter_bxbyez.resize(nlevs_max); } -- cgit v1.2.3 From 6bca0956123895d7641bf710d69aff04419470af Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 17:28:47 -0700 Subject: further cleaning and fix bug when using MR --- Source/Filter/NCIGodfreyFilter.cpp | 1 - Source/Initialization/WarpXInitData.cpp | 2 -- Source/Particles/MultiParticleContainer.H | 11 ---------- Source/Particles/MultiParticleContainer.cpp | 2 -- Source/Particles/PhysicalParticleContainer.cpp | 29 ++++++++++++-------------- Source/WarpX.cpp | 2 +- 6 files changed, 14 insertions(+), 33 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Filter/NCIGodfreyFilter.cpp b/Source/Filter/NCIGodfreyFilter.cpp index 34fca7604..28073725a 100644 --- a/Source/Filter/NCIGodfreyFilter.cpp +++ b/Source/Filter/NCIGodfreyFilter.cpp @@ -55,7 +55,6 @@ void NCIGodfreyFilter::ComputeStencils(){ Real weight_right = cdtodz - index/tab_length; Real prestencil[4]; for(int i=0; ifdtd_nci_stencilz_ex.resize(max_level+1); - mypc->fdtd_nci_stencilz_by.resize(max_level+1); for (int lev = 0; lev <= max_level; ++lev) { const Geometry& gm = Geom(lev); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 0c5e49c04..cc9dc1f59 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -178,17 +178,6 @@ public: void UpdateContinuousInjectionPosition(amrex::Real dt) const; int doContinuousInjection() const; - // - // Parameters for the Cherenkov corrector in the FDTD solver. - // Both stencils are calculated ar runtime. - // - // Number of coefficients for the stencil of the NCI corrector. - // The stencil is applied in the z direction only. - static constexpr int nstencilz_fdtd_nci_corr=5; - - amrex::Vector > fdtd_nci_stencilz_ex; - amrex::Vector > fdtd_nci_stencilz_by; - std::vector GetSpeciesNames() const { return species_names; } PhysicalParticleContainer& GetPCtmp () { return *pc_tmp; } diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 440906348..983530569 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -8,8 +8,6 @@ using namespace amrex; -constexpr int MultiParticleContainer::nstencilz_fdtd_nci_corr; - MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) { ReadParameters(); diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 1e0d68800..d98957c28 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1102,9 +1102,6 @@ PhysicalParticleContainer::Evolve (int lev, const std::array& dx = WarpX::CellSize(lev); const std::array& cdx = WarpX::CellSize(std::max(lev-1,0)); - const auto& mypc = WarpX::GetInstance().GetPartContainer(); - const int nstencilz_fdtd_nci_corr = mypc.nstencilz_fdtd_nci_corr; - // Get instances of NCI Godfrey filters const auto& nci_godfrey_filter_exeybz = WarpX::GetInstance().nci_godfrey_filter_exeybz; const auto& nci_godfrey_filter_bxbyez = WarpX::GetInstance().nci_godfrey_filter_bxbyez; @@ -1181,38 +1178,38 @@ PhysicalParticleContainer::Evolve (int lev, // Safeguard for GPU exeli = filtered_Ex.elixir(); // Apply filter on Ex, result stored in filtered_Ex - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], tbox); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ex, Ex[pti], filtered_Ex.box()); // Update exfab reference exfab = &filtered_Ex; // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], tbox); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Ez, Ez[pti], filtered_Ez.box()); ezfab = &filtered_Ez; // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], tbox); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_By, By[pti], filtered_By.box()); byfab = &filtered_By; #if (AMREX_SPACEDIM == 3) // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], tbox); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Ey, Ey[pti], filtered_Ey.box()); eyfab = &filtered_Ey; // Filter Bx filtered_Bx.resize(amrex::convert(tbox,WarpX::Bx_nodal_flag)); bxeli = filtered_Bx.elixir(); - nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], tbox); + nci_godfrey_filter_bxbyez[lev]->ApplyStencil(filtered_Bx, Bx[pti], filtered_Bx.box()); bxfab = &filtered_Bx; // Filter Bz filtered_Bz.resize(amrex::convert(tbox,WarpX::Bz_nodal_flag)); bzeli = filtered_Bz.elixir(); - nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], tbox); + nci_godfrey_filter_exeybz[lev]->ApplyStencil(filtered_Bz, Bz[pti], filtered_Bz.box()); bzfab = &filtered_Bz; #endif } @@ -1388,44 +1385,44 @@ PhysicalParticleContainer::Evolve (int lev, static_cast(WarpX::noy), static_cast(WarpX::noz)}); #endif - + // Filter Ex (both 2D and 3D) filtered_Ex.resize(amrex::convert(tbox,WarpX::Ex_nodal_flag)); // Safeguard for GPU exeli = filtered_Ex.elixir(); // Apply filter on Ex, result stored in filtered_Ex - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], tbox); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ex, (*cEx)[pti], filtered_Ex.box()); // Update exfab reference cexfab = &filtered_Ex; // Filter Ez filtered_Ez.resize(amrex::convert(tbox,WarpX::Ez_nodal_flag)); ezeli = filtered_Ez.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], tbox); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_Ez, (*cEz)[pti], filtered_Ez.box()); cezfab = &filtered_Ez; // Filter By filtered_By.resize(amrex::convert(tbox,WarpX::By_nodal_flag)); byeli = filtered_By.elixir(); - nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], tbox); + nci_godfrey_filter_bxbyez[lev-1]->ApplyStencil(filtered_By, (*cBy)[pti], filtered_By.box()); cbyfab = &filtered_By; #if (AMREX_SPACEDIM == 3) // Filter Ey filtered_Ey.resize(amrex::convert(tbox,WarpX::Ey_nodal_flag)); eyeli = filtered_Ey.elixir(); - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Ey, (*cEy)[pti], tbox); + 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], tbox); + 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(); - nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], tbox); + nci_godfrey_filter_exeybz[lev-1]->ApplyStencil(filtered_Bz, (*cBz)[pti], filtered_Bz.box()); cbzfab = &filtered_Bz; #endif } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 5f63e2ad6..b24058a0e 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -651,7 +651,7 @@ WarpX::AllocLevelData (int lev, const BoxArray& ba, const DistributionMapping& d int ngz_nonci = (ngz_tmp % 2) ? ngz_tmp+1 : ngz_tmp; // Always even number int ngz; if (WarpX::use_fdtd_nci_corr) { - int ng = ngz_tmp + (mypc->nstencilz_fdtd_nci_corr-1); + int ng = ngz_tmp + 4; ngz = (ng % 2) ? ng+1 : ng; } else { ngz = ngz_nonci; -- cgit v1.2.3 From 97e66052cc96ed40bef74ed3d62bed3fdb26427a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Wed, 8 May 2019 18:11:33 -0700 Subject: remove redundant Elixir declaration --- Source/Particles/PhysicalParticleContainer.cpp | 1 - 1 file changed, 1 deletion(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d98957c28..680086786 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1374,7 +1374,6 @@ PhysicalParticleContainer::Evolve (int lev, const FArrayBox* cbyfab = &(*cBy)[pti]; const FArrayBox* cbzfab = &(*cBz)[pti]; - Elixir exeli, eyeli, ezeli, bxeli, byeli, bzeli; if (WarpX::use_fdtd_nci_corr) { #if (AMREX_SPACEDIM == 2) -- cgit v1.2.3 From 768093cffd18a2c37d90be1033ae94c3ff38776a Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 9 May 2019 15:53:08 -0700 Subject: move particle plot options to PPC --- Source/Diagnostics/FieldIO.cpp | 29 +++++++++++++++++--------- Source/Diagnostics/ParticleIO.cpp | 5 ++--- Source/Diagnostics/WarpXIO.cpp | 2 +- Source/Particles/MultiParticleContainer.H | 1 - Source/Particles/PhysicalParticleContainer.H | 3 +++ Source/Particles/PhysicalParticleContainer.cpp | 20 ++++++++++++++++++ Source/Particles/WarpXParticleContainer.H | 3 +++ Source/WarpX.H | 6 +++--- Source/WarpX.cpp | 5 +++++ 9 files changed, 56 insertions(+), 18 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index 209d8e9b4..ced1f8ea3 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -299,7 +299,10 @@ WarpX::AverageAndPackFields ( Vector& varnames, amrex::Vector& mf_avg, const int ngrow) const { // Count how many different fields should be written (ncomp) - const int ncomp = 3*3 + const int ncomp = 0 + + static_cast(plot_E_field)*3 + + static_cast(plot_B_field)*3 + + static_cast(plot_J_field)*3 + static_cast(plot_part_per_cell) + static_cast(plot_part_per_grid) + static_cast(plot_part_per_proc) @@ -321,15 +324,21 @@ WarpX::AverageAndPackFields ( Vector& varnames, // Go through the different fields, pack them into mf_avg[lev], // add the corresponding names to `varnames` and increment dcomp int dcomp = 0; - AverageAndPackVectorField(mf_avg[lev], current_fp[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"jx","jy","jz"}) varnames.push_back(name); - dcomp += 3; - AverageAndPackVectorField(mf_avg[lev], Efield_aux[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"Ex","Ey","Ez"}) varnames.push_back(name); - dcomp += 3; - AverageAndPackVectorField(mf_avg[lev], Bfield_aux[lev], dcomp, ngrow); - if(lev==0) for(auto name:{"Bx","By","Bz"}) varnames.push_back(name); - dcomp += 3; + if (plot_J_field) { + AverageAndPackVectorField(mf_avg[lev], current_fp[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"jx","jy","jz"}) varnames.push_back(name); + dcomp += 3; + } + if (plot_E_field) { + AverageAndPackVectorField(mf_avg[lev], Efield_aux[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"Ex","Ey","Ez"}) varnames.push_back(name); + dcomp += 3; + } + if (plot_B_field) { + AverageAndPackVectorField(mf_avg[lev], Bfield_aux[lev], dcomp, ngrow); + if(lev==0) for(auto name:{"Bx","By","Bz"}) varnames.push_back(name); + dcomp += 3; + } if (plot_part_per_cell) { diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index f15c084a0..d5de96c96 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -28,7 +28,6 @@ MultiParticleContainer::Checkpoint (const std::string& dir) const void MultiParticleContainer::WritePlotFile (const std::string& dir, - const Vector& real_flags, const Vector& real_names) const { Vector int_names; @@ -36,8 +35,8 @@ MultiParticleContainer::WritePlotFile (const std::string& dir, for (unsigned i = 0, n = species_names.size(); i < n; ++i) { allcontainers[i]->WritePlotFile(dir, species_names[i], - real_flags, int_flags, - real_names, int_names); + allcontainers[i]->plot_flags, int_flags, + real_names, int_names); } } diff --git a/Source/Diagnostics/WarpXIO.cpp b/Source/Diagnostics/WarpXIO.cpp index 186a8d45e..fb81e989f 100644 --- a/Source/Diagnostics/WarpXIO.cpp +++ b/Source/Diagnostics/WarpXIO.cpp @@ -642,7 +642,7 @@ WarpX::WritePlotFile () const particle_varnames.push_back("uzold"); } - mypc->WritePlotFile(plotfilename, particle_plot_flags, particle_varnames); + mypc->WritePlotFile(plotfilename, particle_varnames); WriteJobInfo(plotfilename); diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 0c5e49c04..9291e0358 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -130,7 +130,6 @@ public: void Checkpoint (const std::string& dir) const; void WritePlotFile( const std::string& dir, - const amrex::Vector& real_flags, const amrex::Vector& real_names) const; void Restart (const std::string& dir); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 4f365768b..ac803494d 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -112,6 +112,9 @@ protected: bool boost_adjust_transverse_positions = false; bool do_backward_propagation = false; + amrex::Vector plot_flags; + amrex::Vector plot_vars; + long NumParticlesToAdd (const amrex::Box& overlap_box, const amrex::RealBox& overlap_realbox, const amrex::RealBox& tile_real_box, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 2fa39d87d..c6bc92ff6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -82,6 +82,26 @@ 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); + { + pp.queryarr("plot_vars", plot_vars); + if (plot_vars.size() == 0){ + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + plot_flags.resize(PIdx::nattribs + 6, 1); + } else { + plot_flags.resize(PIdx::nattribs, 1); + } + } else { + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + plot_flags.resize(PIdx::nattribs + 6, 0); + } else { + plot_flags.resize(PIdx::nattribs, 0); + } + + for (const auto& var : plot_vars){ + plot_flags[ParticleStringNames::to_index.at(var)] = 1; + } + } + } } PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core) diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 6fa02b824..d56d14d96 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -272,6 +272,9 @@ protected: amrex::Vector > m_xp, m_yp, m_zp, m_giv; + amrex::Vector plot_flags; + amrex::Vector plot_vars; + private: virtual void particlePostLocate(ParticleType& p, const amrex::ParticleLocData& pld, const int lev) override; diff --git a/Source/WarpX.H b/Source/WarpX.H index 44e84033f..a7a356afc 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -515,6 +515,9 @@ private: bool dump_plotfiles = true; bool dump_openpmd = false; #endif + bool plot_E_field = true; + bool plot_B_field = true; + bool plot_J_field = true; bool plot_part_per_cell = true; bool plot_part_per_grid = false; bool plot_part_per_proc = false; @@ -543,9 +546,6 @@ private: bool is_synchronized = true; - amrex::Vector particle_plot_vars; - amrex::Vector particle_plot_flags; - #ifdef WARPX_USE_PSATD // Store fields in real space on the dual grid (i.e. the grid for the FFT push of the fields) // This includes data for the FFT guard cells (between FFT groups) diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index a3a24897a..f0b990859 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -371,6 +371,9 @@ WarpX::ReadParameters () if (ParallelDescriptor::NProcs() == 1) { plot_proc_number = false; } + pp.query("plot_E_field" , plot_E_field); + pp.query("plot_B_field" , plot_B_field); + pp.query("plot_J_field" , plot_J_field); pp.query("plot_part_per_cell", plot_part_per_cell); pp.query("plot_part_per_grid", plot_part_per_grid); pp.query("plot_part_per_proc", plot_part_per_proc); @@ -429,6 +432,7 @@ WarpX::ReadParameters () } // select which particle comps to write + /* { pp.queryarr("particle_plot_vars", particle_plot_vars); @@ -460,6 +464,7 @@ WarpX::ReadParameters () } } } + */ pp.query("load_balance_int", load_balance_int); pp.query("load_balance_with_sfc", load_balance_with_sfc); -- cgit v1.2.3 From 886d49c55cd179c34197fef2880d9623ce37ce22 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 9 May 2019 19:46:35 -0700 Subject: Can select particle quantities to dump --- Source/Diagnostics/ParticleIO.cpp | 10 +++++-- Source/Particles/PhysicalParticleContainer.H | 2 ++ Source/Particles/PhysicalParticleContainer.cpp | 41 +++++++++++++++++--------- Source/Particles/WarpXParticleContainer.H | 6 ++++ 4 files changed, 42 insertions(+), 17 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index d5de96c96..12e170c24 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -32,11 +32,15 @@ MultiParticleContainer::WritePlotFile (const std::string& dir, { Vector int_names; Vector int_flags; + Vector real_flags; for (unsigned i = 0, n = species_names.size(); i < n; ++i) { - allcontainers[i]->WritePlotFile(dir, species_names[i], - allcontainers[i]->plot_flags, int_flags, - real_names, int_names); + real_flags = allcontainers[i]->plot_flags; + for(int j=0;jWritePlotFile(dir, species_names[i], + real_flags, int_flags, + real_names, int_names); } } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index ac803494d..c9d6b2c2e 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -112,8 +112,10 @@ protected: bool boost_adjust_transverse_positions = false; bool do_backward_propagation = false; + /* amrex::Vector plot_flags; amrex::Vector plot_vars; + */ long NumParticlesToAdd (const amrex::Box& overlap_box, const amrex::RealBox& overlap_realbox, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c6bc92ff6..1bf41816e 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -82,22 +82,35 @@ 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); - { - pp.queryarr("plot_vars", plot_vars); - if (plot_vars.size() == 0){ - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ - plot_flags.resize(PIdx::nattribs + 6, 1); - } else { - plot_flags.resize(PIdx::nattribs, 1); - } - } else { - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ - plot_flags.resize(PIdx::nattribs + 6, 0); - } else { - plot_flags.resize(PIdx::nattribs, 0); - } + 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 + // boosted frame + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + plot_flags.resize(PIdx::nattribs + 6, 1); + } else { + plot_flags.resize(PIdx::nattribs, 1); + } + } else { + // Set plot_flag to 0 for all attribs + if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + plot_flags.resize(PIdx::nattribs + 6, 0); + } else { + // Set plot_flags to 1 for data in plot_vars + plot_flags.resize(PIdx::nattribs, 0); + Print()<<"here2"<.plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } } diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index d56d14d96..51dc5ec05 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -272,7 +272,13 @@ protected: amrex::Vector > m_xp, m_yp, m_zp, m_giv; + // 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 + // to dump to file. amrex::Vector plot_flags; + // list of names of attributes to dump. amrex::Vector plot_vars; private: -- cgit v1.2.3 From 8e5e956f751e47f6eb8435be56dc623ce770c434 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Thu, 9 May 2019 20:08:24 -0700 Subject: some cleaning --- Source/Diagnostics/ParticleIO.cpp | 13 ++++++------- Source/Particles/PhysicalParticleContainer.H | 5 ----- Source/Particles/PhysicalParticleContainer.cpp | 4 +--- 3 files changed, 7 insertions(+), 15 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/ParticleIO.cpp b/Source/Diagnostics/ParticleIO.cpp index 12e170c24..5cc2b9e6e 100644 --- a/Source/Diagnostics/ParticleIO.cpp +++ b/Source/Diagnostics/ParticleIO.cpp @@ -32,15 +32,14 @@ MultiParticleContainer::WritePlotFile (const std::string& dir, { Vector int_names; Vector int_flags; - Vector real_flags; for (unsigned i = 0, n = species_names.size(); i < n; ++i) { - real_flags = allcontainers[i]->plot_flags; - for(int j=0;jWritePlotFile(dir, species_names[i], - real_flags, int_flags, - real_names, int_names); + auto& pc = allcontainers[i]; + if (pc->plot_species) { + pc->WritePlotFile(dir, species_names[i], + pc->plot_flags, int_flags, + real_names, int_names); + } } } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index c9d6b2c2e..4f365768b 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -112,11 +112,6 @@ protected: bool boost_adjust_transverse_positions = false; bool do_backward_propagation = false; - /* - amrex::Vector plot_flags; - amrex::Vector plot_vars; - */ - long NumParticlesToAdd (const amrex::Box& overlap_box, const amrex::RealBox& overlap_realbox, const amrex::RealBox& tile_real_box, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 1bf41816e..b3c598c22 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -100,9 +100,7 @@ PhysicalParticleContainer::PhysicalParticleContainer (AmrCore* amr_core, int isp if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ plot_flags.resize(PIdx::nattribs + 6, 0); } else { - // Set plot_flags to 1 for data in plot_vars plot_flags.resize(PIdx::nattribs, 0); - Print()<<"here2"<.plot_vars argument not in ParticleStringNames"); + "plot_vars argument not in ParticleStringNames"); plot_flags[ParticleStringNames::to_index.at(var)] = 1; } } -- cgit v1.2.3 From 271eb52eafffdefb8dfc3bb43ee521094fbd5a74 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 10 May 2019 17:49:29 -0700 Subject: lab diags is a species argument. add comments --- Source/Diagnostics/BoostedFrameDiagnostic.H | 4 +++ Source/Diagnostics/BoostedFrameDiagnostic.cpp | 25 ++++++++------- Source/Particles/MultiParticleContainer.H | 5 +++ Source/Particles/MultiParticleContainer.cpp | 36 ++++++++++++++++++++-- Source/Particles/PhysicalParticleContainer.cpp | 15 +++++---- .../Particles/RigidInjectedParticleContainer.cpp | 2 +- Source/Particles/WarpXParticleContainer.H | 8 +++++ Source/Particles/WarpXParticleContainer.cpp | 6 ++-- Source/WarpX.cpp | 8 +++-- 9 files changed, 84 insertions(+), 25 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index e35d307a6..ef4bd2ec1 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -75,6 +75,10 @@ class BoostedFrameDiagnostic { int boost_direction_; amrex::Vector > data_buffer_; + // particles_buffer_ is current blind to refinement level. + // particles_buffer_[i][j] is a WarpXParticleContainer::DiagnosticParticleData where + // - i is the back-transformed snapshot number + // - j is the species number amrex::Vector > particles_buffer_; int num_buffer_ = 256; int max_box_size_ = 256; diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 13972075d..5d85fc8f8 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -540,13 +540,14 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) if (WarpX::do_boosted_frame_particles) { for (int j = 0; j < mypc.nSpecies(); ++j) { + std::string species_name = species_names[mypc.map_species_lab_diags[j]]; #ifdef WARPX_USE_HDF5 writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -600,7 +601,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, DistributionMapping buff_dm(buff_ba); data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp, 0) ); } - if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nSpecies()); + if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nspecies_lab_frame_diags); } if (WarpX::do_boosted_frame_fields) { @@ -666,14 +667,15 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) { - for (int j = 0; j < mypc.nSpecies(); ++j) { + for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) { + const std::string species_name = species_names[mypc.map_species_lab_diags[j]]; #ifdef WARPX_USE_HDF5 writeParticleDataHDF5(particles_buffer_[i][j], snapshots_[i].file_name, - species_names[j]); + species_name); #else std::stringstream part_ss; - part_ss << snapshots_[i].file_name + "/" + species_names[j] + "/"; + part_ss << snapshots_[i].file_name + "/" + species_name + "/"; writeParticleData(particles_buffer_[i][j], part_ss.str(), i_lab); #endif } @@ -855,16 +857,16 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, ParallelDescriptor::Barrier(); - if (WarpX::do_boosted_frame_particles) - { + if (WarpX::do_boosted_frame_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); for (int j = 0; j < mypc.nSpecies(); ++j) { - output_create_species_group(file_name, species_names[j]); + std::string species_name = species_names[mypc.map_species_lab_diags[j]]; + output_create_species_group(file_name, species_name); for (int k = 0; k < static_cast(particle_field_names.size()); ++k) { - std::string field_path = species_names[j] + "/" + particle_field_names[k]; + std::string field_path = species_name + "/" + particle_field_names[k]; output_create_particle_field(file_name, field_path); } } @@ -888,7 +890,8 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, const std::string particles_prefix = "particle"; for(int i = 0; i < nspecies; ++i) { - const std::string fullpath = file_name + "/" + species_names[i]; + std::string species_name = species_names[mypc.map_species_lab_diags[i]]; + const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) CreateDirectoryFailed(fullpath); } diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 9291e0358..5a79443d0 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -153,6 +153,8 @@ public: void SetParticleBoxArray (int lev, amrex::BoxArray& new_ba); void SetParticleDistributionMap (int lev, amrex::DistributionMapping& new_dm); + void setSpeciesLabFrameDiags() const; + int nSpecies() const {return nspecies;} int nSpeciesDepositOnMainGrid () const { @@ -184,6 +186,9 @@ public: // Number of coefficients for the stencil of the NCI corrector. // The stencil is applied in the z direction only. static constexpr int nstencilz_fdtd_nci_corr=5; + int nspecies_lab_frame_diags = 0; + std::vector map_species_lab_diags; + int do_boosted_frame_diags = 0; amrex::Vector > fdtd_nci_stencilz_ex; amrex::Vector > fdtd_nci_stencilz_by; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 440906348..bde8d244e 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -31,7 +31,21 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) pc_tmp.reset(new PhysicalParticleContainer(amr_core)); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + // Compute the number of species for which mab frame data is dumped + // nspecies_lab_frame_diags, and map their ID to MultiParticleContainer + // particle IDs in map_species_lab_diags. + map_species_lab_diags.resize(nspecies); + nspecies_lab_frame_diags = 0; + for (int i=0; ido_boosted_frame_diags){ + map_species_lab_diags[nspecies_lab_frame_diags] = i; + nspecies_lab_frame_diags += 1; + do_boosted_frame_diags = 1; + } + } + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { for (int i = 0; i < nspecies + nlasers; ++i) { @@ -376,13 +390,24 @@ MultiParticleContainer BL_PROFILE("MultiParticleContainer::GetLabFrameData"); + // Loop over particle species for (int i = 0; i < nspecies; ++i){ WarpXParticleContainer* pc = allcontainers[i].get(); WarpXParticleContainer::DiagnosticParticles diagnostic_particles; pc->GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); - + // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData + // where "lev" is the AMR level and "index" is a [grid index][tile index] pair. + + // Loop over AMR levels for (int lev = 0; lev <= pc->finestLevel(); ++lev){ + // Loop over [grid index][tile index] pairs + // and Fills parts[species number i] with particle data from all grids and + // tiles in diagnostic_particles. parts contains particles from all + // AMR levels indistinctly. for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){ + // it->first is the [grid index][tile index] key + // it->second is the corresponding + // WarpXParticleContainer::DiagnosticParticleData value parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), it->second.GetRealData(DiagIdx::w ).begin(), it->second.GetRealData(DiagIdx::w ).end()); @@ -459,3 +484,10 @@ MultiParticleContainer::doContinuousInjection() const } return warpx_do_continuous_injection; } + +// Set number of species for which lab frame data is dumped +// and maps their ID to MultiParticleContainer IDs. +//void +//MultiParticleContainer::setSpeciesLabFrameDiags() const +//{ +//} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index b3c598c22..d99cc9c66 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -82,6 +82,9 @@ 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 + // for this species. + pp.query("do_boosted_frame_diags", do_boosted_frame_diags); pp.query("plot_species", plot_species); int do_user_plot_vars; @@ -90,14 +93,14 @@ 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 - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 1); } else { plot_flags.resize(PIdx::nattribs, 1); } } else { // Set plot_flag to 0 for all attribs - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles){ + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags){ plot_flags.resize(PIdx::nattribs + 6, 0); } else { plot_flags.resize(PIdx::nattribs, 0); @@ -216,7 +219,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); @@ -500,7 +503,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -742,7 +745,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -1715,7 +1718,7 @@ PhysicalParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index a5acca281..fd1b2dfb5 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -225,7 +225,7 @@ RigidInjectedParticleContainer::PushPX(WarpXParIter& pti, auto& Bzp = attribs[PIdx::Bz]; const long np = pti.numParticles(); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& xpold = pti.GetAttribs(particle_comps["xold"]); auto& ypold = pti.GetAttribs(particle_comps["yold"]); diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index 51dc5ec05..600061e8d 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -85,7 +85,13 @@ class WarpXParticleContainer public: friend MultiParticleContainer; + // amrex::StructOfArrays with DiagIdx::nattribs amrex::Real components + // for the particle data and 0 int components. 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 + // DiagnosticParticleData (see above) on this tile. using DiagnosticParticles = amrex::Vector, DiagnosticParticleData> >; WarpXParticleContainer (amrex::AmrCore* amr_core, int ispecies); @@ -265,6 +271,8 @@ protected: // support all features allowed by direct injection. int do_continuous_injection = 0; + int do_boosted_frame_diags = 0; + amrex::Vector local_rho; amrex::Vector local_jx; amrex::Vector local_jy; diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 1abbd747d..0cf5c10b4 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -78,7 +78,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -231,7 +231,7 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); @@ -249,7 +249,7 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && WarpX::do_boosted_frame_particles) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index eb84af2c7..c863e9a4c 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -56,9 +56,12 @@ int WarpX::num_mirrors = 0; int WarpX::sort_int = -1; -bool WarpX::do_boosted_frame_diagnostic = false; int WarpX::num_snapshots_lab = std::numeric_limits::lowest(); Real WarpX::dt_snapshots_lab = std::numeric_limits::lowest(); +// bool WarpX::do_boosted_frame_diagnostic = false; +// bool WarpX::do_boosted_frame_fields = true; +// bool WarpX::do_boosted_frame_particles = true; +bool WarpX::do_boosted_frame_diagnostic = false; bool WarpX::do_boosted_frame_fields = true; bool WarpX::do_boosted_frame_particles = true; @@ -117,7 +120,7 @@ WarpX::ResetInstance () { delete m_instance; m_instance = nullptr; -} +} WarpX::WarpX () { @@ -157,6 +160,7 @@ WarpX::WarpX () current_injection_position = geom[0].ProbLo(moving_window_dir); } } + do_boosted_frame_particles = mypc->do_boosted_frame_diags; Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); -- cgit v1.2.3 From 71cdf5ca0bfe5a239382cdfb56df1d2d9e68c65f Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 12 May 2019 06:54:06 -0700 Subject: debugging --- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 17 ++++++++++++----- Source/Evolve/WarpXEvolveEM.cpp | 4 ++++ Source/Particles/MultiParticleContainer.cpp | 1 + Source/Particles/PhysicalParticleContainer.cpp | 2 ++ 4 files changed, 19 insertions(+), 5 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index ad33b68d8..934610b9e 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -498,7 +498,7 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) { BL_PROFILE("BoostedFrameDiagnostic::Flush"); - std::cout<<"in Flush"\n; + std::cout<<"in Flush\n"; VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); @@ -546,6 +546,7 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) int js = mypc.map_species_lab_diags[j]; std::string species_name = species_names[js]; #ifdef WARPX_USE_HDF5 + std::cout<<"particles_buffer_ j "< zhi_boost) or (snapshots_[i].current_z_lab < snapshots_[i].zmin_lab) or (snapshots_[i].current_z_lab > snapshots_[i].zmax_lab) ) continue; + std::cout<<"toto 3\n"; int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 3\n"; @@ -656,9 +661,9 @@ writeLabFrameData(const MultiFab* cell_centered_data, if (WarpX::do_boosted_frame_particles) { std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 6\n"; - mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, - old_z_boost, snapshots_[i].current_z_boost, - t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); + //mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, + // old_z_boost, snapshots_[i].current_z_boost, + // t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 7\n"; } @@ -688,6 +693,7 @@ writeLabFrameData(const MultiFab* cell_centered_data, const std::string species_name = species_names[mypc.map_species_lab_diags[j]]; #ifdef WARPX_USE_HDF5 + std::cout<<"particles_buffer_ j "<writeLabFrameData"<writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); + std::cout<<"after myBFD->writeLabFrameData"<Flush"<Flush(geom[0]); + std::cout<<"after myBFD->Flush"<GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); std::cout<<"GetLabFrameData 3\n"; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index d99cc9c66..c4d4abfd6 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1863,6 +1863,8 @@ void PhysicalParticleContainer::GetParticleSlice(const int direction, const Real // Note the the slice should always move in the negative boost direction. AMREX_ALWAYS_ASSERT(z_new < z_old); + AMREX_ALWAYS_ASSERT(do_boosted_frame_diags == 1); + const int nlevs = std::max(0, finestLevel()+1); // we figure out a box for coarse-grained rejection. If the RealBox corresponding to a -- cgit v1.2.3 From c952dbe2e3d9a2c7bab2774e8036e0e9ff72e0ed Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 12 May 2019 07:58:38 -0700 Subject: only selected species BFD-dumped, but all old attribs initialized --- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 30 ++++---------------------- Source/Particles/MultiParticleContainer.cpp | 19 ++++++++++------ Source/Particles/PhysicalParticleContainer.cpp | 9 +++++--- Source/Particles/WarpXParticleContainer.cpp | 9 +++++--- Source/WarpX.cpp | 2 +- 5 files changed, 29 insertions(+), 40 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 934610b9e..7a44bc66a 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -582,27 +582,19 @@ writeLabFrameData(const MultiFab* cell_centered_data, const std::vector species_names = mypc.GetSpeciesNames(); - Print()<<"in BoostedFrameDiagnostic::writeLabFrameData 1\n"; - for (int i = 0; i < N_snapshots_; ++i) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 2\n"; - std::cout<<"i "< zhi_boost) or (snapshots_[i].current_z_lab < snapshots_[i].zmin_lab) or (snapshots_[i].current_z_lab > snapshots_[i].zmax_lab) ) continue; - std::cout<<"toto 3\n"; int i_lab = (snapshots_[i].current_z_lab - snapshots_[i].zmin_lab) / dz_lab_; - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 3\n"; if (buff_counter_[i] == 0) { if (WarpX::do_boosted_frame_fields) { @@ -617,7 +609,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) particles_buffer_[i].resize(mypc.nspecies_lab_frame_diags); } - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 4\n"; if (WarpX::do_boosted_frame_fields) { const int ncomp = cell_centered_data->nComp(); @@ -657,14 +648,11 @@ writeLabFrameData(const MultiFab* cell_centered_data, &ncomp, &i_boost, &i_lab); } } - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 5\n"; if (WarpX::do_boosted_frame_particles) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 6\n"; - //mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, - // old_z_boost, snapshots_[i].current_z_boost, - // t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 7\n"; + mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, + old_z_boost, snapshots_[i].current_z_boost, + t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); } @@ -686,28 +674,21 @@ writeLabFrameData(const MultiFab* cell_centered_data, } if (WarpX::do_boosted_frame_particles) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 8\n"; for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) { - std::cout<<"in BoostedFrameDiagnostic::writeLabFrameData 9\n"; const std::string species_name = species_names[mypc.map_species_lab_diags[j]]; #ifdef WARPX_USE_HDF5 - std::cout<<"particles_buffer_ j "<(particle_field_names.size()); ++k) @@ -919,10 +899,8 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); -// int nspecies = mypc.nSpecies(); const std::string particles_prefix = "particle"; -// for(int i = 0; i < nspecies; ++i) { for(int i = 0; i < mypc.nspecies_lab_frame_diags; ++i) { int is = mypc.map_species_lab_diags[i]; std::string species_name = species_names[is]; diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 841c41835..337005dc3 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -55,7 +55,19 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[i]->AddRealComp("uxold"); allcontainers[i]->AddRealComp("uyold"); allcontainers[i]->AddRealComp("uzold"); + } + /* + for (int i = 0; i < nspecies_lab_frame_diags; ++i) + { + int is = map_species_lab_diags[i]; + allcontainers[is]->AddRealComp("xold"); + allcontainers[is]->AddRealComp("yold"); + allcontainers[is]->AddRealComp("zold"); + allcontainers[is]->AddRealComp("uxold"); + allcontainers[is]->AddRealComp("uyold"); + allcontainers[is]->AddRealComp("uzold"); } + */ pc_tmp->AddRealComp("xold"); pc_tmp->AddRealComp("yold"); pc_tmp->AddRealComp("zold"); @@ -389,17 +401,13 @@ MultiParticleContainer { BL_PROFILE("MultiParticleContainer::GetLabFrameData"); - std::cout<<"GetLabFrameData 1\n"; // Loop over particle species for (int i = 0; i < nspecies_lab_frame_diags; ++i){ int isp = map_species_lab_diags[i]; -std::cout<<"GetLabFrameData 2\n"; WarpXParticleContainer* pc = allcontainers[isp].get(); - std::cout<<"getparticleslice "<< isp<GetParticleSlice(direction, z_old, z_new, t_boost, t_lab, dt, diagnostic_particles); -std::cout<<"GetLabFrameData 3\n"; // Here, diagnostic_particles[lev][index] is a WarpXParticleContainer::DiagnosticParticleData // where "lev" is the AMR level and "index" is a [grid index][tile index] pair. @@ -409,16 +417,13 @@ std::cout<<"GetLabFrameData 3\n"; // and Fills parts[species number i] with particle data from all grids and // tiles in diagnostic_particles. parts contains particles from all // AMR levels indistinctly. -std::cout<<"GetLabFrameData 4\n"; for (auto it = diagnostic_particles[lev].begin(); it != diagnostic_particles[lev].end(); ++it){ // it->first is the [grid index][tile index] key // it->second is the corresponding // WarpXParticleContainer::DiagnosticParticleData value -std::cout<<"GetLabFrameData 5\n"; parts[i].GetRealData(DiagIdx::w).insert( parts[i].GetRealData(DiagIdx::w ).end(), it->second.GetRealData(DiagIdx::w ).begin(), it->second.GetRealData(DiagIdx::w ).end()); -std::cout<<"GetLabFrameData 6\n"; parts[i].GetRealData(DiagIdx::x).insert( parts[i].GetRealData(DiagIdx::x ).end(), it->second.GetRealData(DiagIdx::x ).begin(), diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c4d4abfd6..a3f555bd0 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -219,7 +219,8 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); @@ -503,7 +504,8 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -745,7 +747,8 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 0cf5c10b4..9bee0031a 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -78,7 +78,8 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -231,7 +232,8 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); @@ -249,7 +251,8 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + if (WarpX::do_boosted_frame_diagnostic) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index c863e9a4c..7c026e105 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -329,7 +329,7 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); pp.query("do_boosted_frame_fields", do_boosted_frame_fields); - pp.query("do_boosted_frame_particles", do_boosted_frame_particles); + // pp.query("do_boosted_frame_particles", do_boosted_frame_particles); AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, -- cgit v1.2.3 From 9a26a71845fde091c7840772bef1b23dbc46d6ac Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 12 May 2019 10:17:28 -0700 Subject: old attribs not allocated if species not BFD --- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 22 ++++++++++------------ Source/Evolve/WarpXEvolveEM.cpp | 2 -- Source/Laser/LaserParticleContainer.cpp | 5 +++-- Source/Particles/MultiParticleContainer.H | 11 ++++++++--- Source/Particles/MultiParticleContainer.cpp | 7 +++++-- Source/Particles/PhysicalParticleContainer.cpp | 12 +++++------- Source/Particles/WarpXParticleContainer.H | 1 - Source/Particles/WarpXParticleContainer.cpp | 19 +++++++++++++------ Source/WarpX.cpp | 2 +- 9 files changed, 45 insertions(+), 36 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 7a44bc66a..413730f91 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -497,8 +497,6 @@ BoostedFrameDiagnostic(Real zmin_lab, Real zmax_lab, Real v_window_lab, void BoostedFrameDiagnostic::Flush(const Geometry& geom) { BL_PROFILE("BoostedFrameDiagnostic::Flush"); - - std::cout<<"in Flush\n"; VisMF::Header::Version current_version = VisMF::GetHeaderVersion(); VisMF::SetHeaderVersion(amrex::VisMF::Header::NoFabHeader_v1); @@ -542,11 +540,10 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) if (WarpX::do_boosted_frame_particles) { // for (int j = 0; j < mypc.nSpecies(); ++j) { - for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) { - int js = mypc.map_species_lab_diags[j]; + for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { + int js = mypc.mapSpeciesLabDiags(j); std::string species_name = species_names[js]; #ifdef WARPX_USE_HDF5 - std::cout<<"particles_buffer_ j "<(particle_field_names.size()); ++k) { std::string field_path = species_name + "/" + particle_field_names[k]; @@ -872,9 +870,9 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); // for (int j = 0; j < mypc.nSpecies(); ++j) - for (int j = 0; j < mypc.nspecies_lab_frame_diags; ++j) + for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { - int js = mypc.map_species_lab_diags[j]; + int js = mypc.mapSpeciesLabDiags(j); std::string species_name = species_names[js]; output_create_species_group(file_name, species_name); for (int k = 0; k < static_cast(particle_field_names.size()); ++k) @@ -901,8 +899,8 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, const std::vector species_names = mypc.GetSpeciesNames(); const std::string particles_prefix = "particle"; - for(int i = 0; i < mypc.nspecies_lab_frame_diags; ++i) { - int is = mypc.map_species_lab_diags[i]; + for(int i = 0; i < mypc.nSpeciesLabFrameDiags(); ++i) { + int is = mypc.mapSpeciesLabDiags(i); std::string species_name = species_names[is]; const std::string fullpath = file_name + "/" + species_name; if (!UtilCreateDirectory(fullpath, 0755)) diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index 5f6760f22..d88f178df 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -173,9 +173,7 @@ WarpX::EvolveEM (int numsteps) if (WarpX::do_boosted_frame_fields) { cell_centered_data = GetCellCenteredData(); } - std::cout<<"before myBFD->writeLabFrameData"<writeLabFrameData(cell_centered_data.get(), *mypc, geom[0], cur_time, dt[0]); - std::cout<<"after myBFD->writeLabFrameData"<::max(); - - ParmParse pp(laser_name); + do_boosted_frame_diags = 0; + + ParmParse pp(laser_name); // Parse the type of laser profile and set the corresponding flag `profile` std::string laser_type_s; diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 5a79443d0..217f727b0 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -157,6 +157,10 @@ public: int nSpecies() const {return nspecies;} + int nSpeciesLabFrameDiags() const {return nspecies_lab_frame_diags;} + int mapSpeciesLabDiags(int i) const {return map_species_lab_diags[i];} + int doBoostedFrameDiags() const {return do_boosted_frame_diags;} + int nSpeciesDepositOnMainGrid () const { int r = 0; for (int i : deposit_on_main_grid) { @@ -186,9 +190,6 @@ public: // Number of coefficients for the stencil of the NCI corrector. // The stencil is applied in the z direction only. static constexpr int nstencilz_fdtd_nci_corr=5; - int nspecies_lab_frame_diags = 0; - std::vector map_species_lab_diags; - int do_boosted_frame_diags = 0; amrex::Vector > fdtd_nci_stencilz_ex; amrex::Vector > fdtd_nci_stencilz_by; @@ -219,6 +220,10 @@ private: void ReadParameters (); + int nspecies_lab_frame_diags = 0; + std::vector map_species_lab_diags; + int do_boosted_frame_diags = 0; + // runtime parameters int nlasers = 0; int nspecies = 1; // physical particles only. nspecies+nlasers == allcontainers.size(). diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index 337005dc3..508f3f606 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -47,6 +47,8 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { + //maxoldattribs + /* for (int i = 0; i < nspecies + nlasers; ++i) { allcontainers[i]->AddRealComp("xold"); @@ -56,7 +58,7 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[i]->AddRealComp("uyold"); allcontainers[i]->AddRealComp("uzold"); } - /* + */ for (int i = 0; i < nspecies_lab_frame_diags; ++i) { int is = map_species_lab_diags[i]; @@ -67,13 +69,14 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[is]->AddRealComp("uyold"); allcontainers[is]->AddRealComp("uzold"); } - */ + /* pc_tmp->AddRealComp("xold"); pc_tmp->AddRealComp("yold"); pc_tmp->AddRealComp("zold"); pc_tmp->AddRealComp("uxold"); pc_tmp->AddRealComp("uyold"); pc_tmp->AddRealComp("uzold"); + */ } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index a3f555bd0..dd167ba41 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -219,8 +219,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, attribs[PIdx::uz] = u[2]; attribs[PIdx::w ] = weight; - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); @@ -504,8 +503,7 @@ PhysicalParticleContainer::AddPlasmaCPU (int lev, RealBox part_realbox) attribs[PIdx::uy] = u[1]; attribs[PIdx::uz] = u[2]; - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -747,8 +745,7 @@ PhysicalParticleContainer::AddPlasmaGPU (int lev, RealBox part_realbox) attribs[PIdx::uz] = u[2]; // note - this will be slow on the GPU, need to revisit - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(lev, grid_id, tile_id); particle_tile.push_back_real(particle_comps["xold"], x); @@ -838,8 +835,9 @@ FieldGatherES (const amrex::Vector, const auto& particles = pti.GetArrayOfStructs(); int nstride = particles.dataShape().first; const long np = pti.numParticles(); - + std::cout<<"start 1 GetAttribs\n"; auto& attribs = pti.GetAttribs(); + std::cout<<"end 2 GetAttribs\n"; auto& Exp = attribs[PIdx::Ex]; auto& Eyp = attribs[PIdx::Ey]; #if AMREX_SPACEDIM == 3 diff --git a/Source/Particles/WarpXParticleContainer.H b/Source/Particles/WarpXParticleContainer.H index d1e25f3ad..6968f2a61 100644 --- a/Source/Particles/WarpXParticleContainer.H +++ b/Source/Particles/WarpXParticleContainer.H @@ -59,7 +59,6 @@ public: const amrex::Cuda::ManagedDeviceVector& y, const amrex::Cuda::ManagedDeviceVector& z); #endif - const std::array& GetAttribs () const { return GetStructOfArrays().GetRealData(); } diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 9bee0031a..66f0dfb5c 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -6,6 +6,8 @@ #include #include #include +#include +#include using namespace amrex; @@ -22,7 +24,9 @@ WarpXParIter::GetPosition (Cuda::ManagedDeviceVector& x, Cuda::ManagedDevi { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); #ifdef WARPX_RZ + std::cout<<"start 2 GetAttribs()\n"; const auto& attribs = GetAttribs(); + std::cout<<"stop 2 GetAttribs()\n"; const auto& theta = attribs[PIdx::theta]; y.resize(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -39,7 +43,9 @@ void WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector& x, const Cuda::ManagedDeviceVector& y, const Cuda::ManagedDeviceVector& z) { #ifdef WARPX_RZ + std::cout<<"start 3 GetAttribs()\n"; auto& attribs = GetAttribs(); + std::cout<<"stop 3 GetAttribs()\n"; auto& theta = attribs[PIdx::theta]; Cuda::DeviceVector r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -78,8 +84,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) particle_comps["theta"] = PIdx::theta; #endif - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { particle_comps["xold"] = PIdx::nattribs; particle_comps["yold"] = PIdx::nattribs+1; @@ -232,8 +237,7 @@ WarpXParticleContainer::AddNParticles (int lev, p.pos(1) = z[i]; #endif - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x[i]); @@ -251,8 +255,7 @@ WarpXParticleContainer::AddNParticles (int lev, particle_tile.push_back_real(PIdx::uy, vy + ibegin, vy + iend); particle_tile.push_back_real(PIdx::uz, vz + ibegin, vz + iend); - // if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - if (WarpX::do_boosted_frame_diagnostic) + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); @@ -998,7 +1001,9 @@ WarpXParticleContainer::PushXES (Real dt) int nstride = particles.dataShape().first; const long np = pti.numParticles(); + std::cout<<"start 4 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); + std::cout<<"stop 4 GetAttribs()\n"; auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; @@ -1047,7 +1052,9 @@ WarpXParticleContainer::PushX (int lev, Real dt) // - positions are stored as an array of struct, in `ParticleType` ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); // - momenta are stored as a struct of array, in `attribs` + std::cout<<"start 5 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); + std::cout<<"stop 5 GetAttribs()\n"; Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 7c026e105..e4c54cc05 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -160,7 +160,7 @@ WarpX::WarpX () current_injection_position = geom[0].ProbLo(moving_window_dir); } } - do_boosted_frame_particles = mypc->do_boosted_frame_diags; + do_boosted_frame_particles = mypc->doBoostedFrameDiags(); Efield_aux.resize(nlevs_max); Bfield_aux.resize(nlevs_max); -- cgit v1.2.3 From 449c1254dfd126b7b2b6291b05177cb22ed38f5c Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Sun, 12 May 2019 10:35:02 -0700 Subject: cleaning (remove print statements etc.) --- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 5 +---- Source/Evolve/WarpXEvolveEM.cpp | 2 -- Source/Particles/MultiParticleContainer.H | 2 -- Source/Particles/MultiParticleContainer.cpp | 21 --------------------- Source/Particles/PhysicalParticleContainer.cpp | 2 -- Source/Particles/WarpXParticleContainer.cpp | 10 ---------- Source/WarpX.cpp | 7 +------ 7 files changed, 2 insertions(+), 47 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index 413730f91..49038b9e2 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -539,7 +539,6 @@ void BoostedFrameDiagnostic::Flush(const Geometry& geom) } if (WarpX::do_boosted_frame_particles) { - // for (int j = 0; j < mypc.nSpecies(); ++j) { for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { int js = mypc.mapSpeciesLabDiags(j); std::string species_name = species_names[js]; @@ -645,8 +644,8 @@ writeLabFrameData(const MultiFab* cell_centered_data, &ncomp, &i_boost, &i_lab); } } - if (WarpX::do_boosted_frame_particles) { + if (WarpX::do_boosted_frame_particles) { mypc.GetLabFrameData(snapshots_[i].file_name, i_lab, boost_direction_, old_z_boost, snapshots_[i].current_z_boost, t_boost, snapshots_[i].t_lab, dt, particles_buffer_[i]); @@ -732,7 +731,6 @@ writeParticleDataHDF5(const WarpXParticleContainer::DiagnosticParticleData& pdat ParallelDescriptor::ReduceLongMax(old_np); // Write data here - Print()<<"particle_field_names.size()"<(particle_field_names.size()); ++k) { std::string field_path = species_name + "/" + particle_field_names[k]; @@ -869,7 +867,6 @@ LabSnapShot(Real t_lab_in, Real t_boost, Real zmin_lab_in, if (WarpX::do_boosted_frame_particles){ auto & mypc = WarpX::GetInstance().GetPartContainer(); const std::vector species_names = mypc.GetSpeciesNames(); - // for (int j = 0; j < mypc.nSpecies(); ++j) for (int j = 0; j < mypc.nSpeciesLabFrameDiags(); ++j) { int js = mypc.mapSpeciesLabDiags(j); diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index d88f178df..dab58f95b 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -247,9 +247,7 @@ WarpX::EvolveEM (int numsteps) } if (do_boosted_frame_diagnostic) { - std::cout<<"before myBFD->Flush"<Flush(geom[0]); - std::cout<<"after myBFD->Flush"<AddRealComp("xold"); - allcontainers[i]->AddRealComp("yold"); - allcontainers[i]->AddRealComp("zold"); - allcontainers[i]->AddRealComp("uxold"); - allcontainers[i]->AddRealComp("uyold"); - allcontainers[i]->AddRealComp("uzold"); - } - */ for (int i = 0; i < nspecies_lab_frame_diags; ++i) { int is = map_species_lab_diags[i]; @@ -69,14 +57,12 @@ MultiParticleContainer::MultiParticleContainer (AmrCore* amr_core) allcontainers[is]->AddRealComp("uyold"); allcontainers[is]->AddRealComp("uzold"); } - /* pc_tmp->AddRealComp("xold"); pc_tmp->AddRealComp("yold"); pc_tmp->AddRealComp("zold"); pc_tmp->AddRealComp("uxold"); pc_tmp->AddRealComp("uyold"); pc_tmp->AddRealComp("uzold"); - */ } } @@ -500,10 +486,3 @@ MultiParticleContainer::doContinuousInjection() const } return warpx_do_continuous_injection; } - -// Set number of species for which lab frame data is dumped -// and maps their ID to MultiParticleContainer IDs. -//void -//MultiParticleContainer::setSpeciesLabFrameDiags() const -//{ -//} diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index dd167ba41..37c136a3d 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -835,9 +835,7 @@ FieldGatherES (const amrex::Vector, const auto& particles = pti.GetArrayOfStructs(); int nstride = particles.dataShape().first; const long np = pti.numParticles(); - std::cout<<"start 1 GetAttribs\n"; auto& attribs = pti.GetAttribs(); - std::cout<<"end 2 GetAttribs\n"; auto& Exp = attribs[PIdx::Ex]; auto& Eyp = attribs[PIdx::Ey]; #if AMREX_SPACEDIM == 3 diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 66f0dfb5c..0cf5c10b4 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -6,8 +6,6 @@ #include #include #include -#include -#include using namespace amrex; @@ -24,9 +22,7 @@ WarpXParIter::GetPosition (Cuda::ManagedDeviceVector& x, Cuda::ManagedDevi { amrex::ParIter<0,0,PIdx::nattribs>::GetPosition(x, z); #ifdef WARPX_RZ - std::cout<<"start 2 GetAttribs()\n"; const auto& attribs = GetAttribs(); - std::cout<<"stop 2 GetAttribs()\n"; const auto& theta = attribs[PIdx::theta]; y.resize(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -43,9 +39,7 @@ void WarpXParIter::SetPosition (const Cuda::ManagedDeviceVector& x, const Cuda::ManagedDeviceVector& y, const Cuda::ManagedDeviceVector& z) { #ifdef WARPX_RZ - std::cout<<"start 3 GetAttribs()\n"; auto& attribs = GetAttribs(); - std::cout<<"stop 3 GetAttribs()\n"; auto& theta = attribs[PIdx::theta]; Cuda::DeviceVector r(x.size()); for (unsigned int i=0 ; i < x.size() ; i++) { @@ -1001,9 +995,7 @@ WarpXParticleContainer::PushXES (Real dt) int nstride = particles.dataShape().first; const long np = pti.numParticles(); - std::cout<<"start 4 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); - std::cout<<"stop 4 GetAttribs()\n"; auto& uxp = attribs[PIdx::ux]; auto& uyp = attribs[PIdx::uy]; auto& uzp = attribs[PIdx::uz]; @@ -1052,9 +1044,7 @@ WarpXParticleContainer::PushX (int lev, Real dt) // - positions are stored as an array of struct, in `ParticleType` ParticleType * AMREX_RESTRICT pstructs = &(pti.GetArrayOfStructs()[0]); // - momenta are stored as a struct of array, in `attribs` - std::cout<<"start 5 GetAttribs()\n"; auto& attribs = pti.GetAttribs(); - std::cout<<"stop 5 GetAttribs()\n"; Real* AMREX_RESTRICT ux = attribs[PIdx::ux].dataPtr(); Real* AMREX_RESTRICT uy = attribs[PIdx::uy].dataPtr(); Real* AMREX_RESTRICT uz = attribs[PIdx::uz].dataPtr(); diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index e4c54cc05..6b1edaccd 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -56,12 +56,9 @@ int WarpX::num_mirrors = 0; int WarpX::sort_int = -1; +bool WarpX::do_boosted_frame_diagnostic = false; int WarpX::num_snapshots_lab = std::numeric_limits::lowest(); Real WarpX::dt_snapshots_lab = std::numeric_limits::lowest(); -// bool WarpX::do_boosted_frame_diagnostic = false; -// bool WarpX::do_boosted_frame_fields = true; -// bool WarpX::do_boosted_frame_particles = true; -bool WarpX::do_boosted_frame_diagnostic = false; bool WarpX::do_boosted_frame_fields = true; bool WarpX::do_boosted_frame_particles = true; @@ -329,8 +326,6 @@ WarpX::ReadParameters () pp.get("gamma_boost", gamma_boost); pp.query("do_boosted_frame_fields", do_boosted_frame_fields); - // pp.query("do_boosted_frame_particles", do_boosted_frame_particles); - AMREX_ALWAYS_ASSERT_WITH_MESSAGE(do_moving_window, "The moving window should be on if using the boosted frame diagnostic."); -- cgit v1.2.3 From 5308cbefa8fbde30bed8a88943d58ec646731d9e Mon Sep 17 00:00:00 2001 From: Andrew Myers Date: Thu, 23 May 2019 15:27:17 -0700 Subject: fix a bunch of unused variable / parameter shadowing warnings --- Source/Diagnostics/BoostedFrameDiagnostic.H | 2 +- Source/Diagnostics/BoostedFrameDiagnostic.cpp | 5 --- Source/Diagnostics/FieldIO.cpp | 1 - Source/Evolve/WarpXEvolveEM.cpp | 28 +++++++------- Source/FieldSolver/WarpXPushFieldsEM.cpp | 48 ++++++++++++------------ Source/FortranInterface/WarpX_picsar.F90 | 7 +++- Source/Laser/LaserParticleContainer.cpp | 6 +-- Source/Parallelization/WarpXComm.cpp | 52 +++++++++++++------------- Source/Particles/PhysicalParticleContainer.cpp | 3 -- Source/Particles/WarpXParticleContainer.cpp | 44 +++++++++++----------- Source/Utils/WarpXMovingWindow.cpp | 6 +-- Source/WarpX.cpp | 1 - 12 files changed, 96 insertions(+), 107 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.H b/Source/Diagnostics/BoostedFrameDiagnostic.H index 55a124c51..3a09259b0 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.H +++ b/Source/Diagnostics/BoostedFrameDiagnostic.H @@ -40,8 +40,8 @@ class BoostedFrameDiagnostic { amrex::Real current_z_lab; amrex::Real current_z_boost; - std::vector mesh_field_names_; int ncomp_to_dump_; + std::vector mesh_field_names_; int file_num; int initial_i; diff --git a/Source/Diagnostics/BoostedFrameDiagnostic.cpp b/Source/Diagnostics/BoostedFrameDiagnostic.cpp index f62369091..57872790f 100644 --- a/Source/Diagnostics/BoostedFrameDiagnostic.cpp +++ b/Source/Diagnostics/BoostedFrameDiagnostic.cpp @@ -497,8 +497,6 @@ LorentzTransformZ(MultiFab& data, Real gamma_boost, Real beta_boost, int ncomp) // Transform the transverse E and B fields. Note that ez and bz are not // changed by the tranform. Real e_lab, b_lab, j_lab, r_lab; - int i0 = 0; - int i4 = 4; e_lab = gamma_boost * (arr(i, j, k, 0) + beta_boost*clight*arr(i, j, k, 4)); b_lab = gamma_boost * (arr(i, j, k, 4) + beta_boost*arr(i, j, k, 0)/clight); @@ -718,7 +716,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, if (buff_counter_[i] == 0) { // ... reset fields buffer data_buffer_[i] if (WarpX::do_boosted_frame_fields) { - const int ncomp = cell_centered_data->nComp(); Box buff_box = geom.Domain(); buff_box.setSmall(boost_direction_, i_lab - num_buffer_ + 1); buff_box.setBig(boost_direction_, i_lab); @@ -726,7 +723,6 @@ writeLabFrameData(const MultiFab* cell_centered_data, buff_ba.maxSize(max_box_size_); DistributionMapping buff_dm(buff_ba); data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp_to_dump, 0) ); - // data_buffer_[i].reset( new MultiFab(buff_ba, buff_dm, ncomp, 0) ); } // ... reset particle buffer particles_buffer_[i] if (WarpX::do_boosted_frame_particles) @@ -961,7 +957,6 @@ LabSnapShot(Real t_lab_in, Real t_boost, RealBox prob_domain_lab, my_bfd(bfd) { Real zmin_lab = prob_domain_lab_.lo(AMREX_SPACEDIM-1); - Real zmax_lab = prob_domain_lab_.hi(AMREX_SPACEDIM-1); current_z_lab = 0.0; current_z_boost = 0.0; updateCurrentZPositions(t_boost, my_bfd.inv_gamma_boost_, my_bfd.inv_beta_boost_); diff --git a/Source/Diagnostics/FieldIO.cpp b/Source/Diagnostics/FieldIO.cpp index ced1f8ea3..e3d44d1fc 100644 --- a/Source/Diagnostics/FieldIO.cpp +++ b/Source/Diagnostics/FieldIO.cpp @@ -660,7 +660,6 @@ getInterpolatedScalar( interpolated_F->setVal(0.); // Loop through the boxes and interpolate the values from the _cp data - const int use_limiter = 0; #ifdef _OPEMP #pragma omp parallel #endif diff --git a/Source/Evolve/WarpXEvolveEM.cpp b/Source/Evolve/WarpXEvolveEM.cpp index dab58f95b..ad7c7d840 100644 --- a/Source/Evolve/WarpXEvolveEM.cpp +++ b/Source/Evolve/WarpXEvolveEM.cpp @@ -514,7 +514,7 @@ WarpX::ComputeDt () * simulation box passes input parameter zmax_plasma_to_compute_max_step. */ void -WarpX::computeMaxStepBoostAccelerator(amrex::Geometry geom){ +WarpX::computeMaxStepBoostAccelerator(amrex::Geometry a_geom){ // Sanity checks: can use zmax_plasma_to_compute_max_step only if // the moving window and the boost are all in z direction. AMREX_ALWAYS_ASSERT_WITH_MESSAGE( @@ -534,7 +534,7 @@ WarpX::computeMaxStepBoostAccelerator(amrex::Geometry geom){ // Lower end of the simulation domain. All quantities are given in boosted // frame except zmax_plasma_to_compute_max_step. - const Real zmin_domain_boost = geom.ProbLo(AMREX_SPACEDIM-1); + const Real zmin_domain_boost = a_geom.ProbLo(AMREX_SPACEDIM-1); // End of the plasma: Transform input argument // zmax_plasma_to_compute_max_step to boosted frame. const Real len_plasma_boost = zmax_plasma_to_compute_max_step/gamma_boost; @@ -591,19 +591,19 @@ WarpX::applyMirrors(Real time){ NullifyMF(Bz, lev, z_min, z_max); if (lev>0){ // Get coarse patch field MultiFabs - MultiFab& Ex = *Efield_cp[lev][0].get(); - MultiFab& Ey = *Efield_cp[lev][1].get(); - MultiFab& Ez = *Efield_cp[lev][2].get(); - MultiFab& Bx = *Bfield_cp[lev][0].get(); - MultiFab& By = *Bfield_cp[lev][1].get(); - MultiFab& Bz = *Bfield_cp[lev][2].get(); + MultiFab& cEx = *Efield_cp[lev][0].get(); + MultiFab& cEy = *Efield_cp[lev][1].get(); + MultiFab& cEz = *Efield_cp[lev][2].get(); + MultiFab& cBx = *Bfield_cp[lev][0].get(); + MultiFab& cBy = *Bfield_cp[lev][1].get(); + MultiFab& cBz = *Bfield_cp[lev][2].get(); // Set each field to zero between z_min and z_max - NullifyMF(Ex, lev, z_min, z_max); - NullifyMF(Ey, lev, z_min, z_max); - NullifyMF(Ez, lev, z_min, z_max); - NullifyMF(Bx, lev, z_min, z_max); - NullifyMF(By, lev, z_min, z_max); - NullifyMF(Bz, lev, z_min, z_max); + NullifyMF(cEx, lev, z_min, z_max); + NullifyMF(cEy, lev, z_min, z_max); + NullifyMF(cEz, lev, z_min, z_max); + NullifyMF(cBx, lev, z_min, z_max); + NullifyMF(cBy, lev, z_min, z_max); + NullifyMF(cBz, lev, z_min, z_max); } } } diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp index a66d14980..c53e13f8f 100644 --- a/Source/FieldSolver/WarpXPushFieldsEM.cpp +++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp @@ -17,30 +17,30 @@ using namespace amrex; void -WarpX::EvolveB (Real dt) +WarpX::EvolveB (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveB(lev, dt); + EvolveB(lev, a_dt); } } void -WarpX::EvolveB (int lev, Real dt) +WarpX::EvolveB (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveB()"); - EvolveB(lev, PatchType::fine, dt); + EvolveB(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveB(lev, PatchType::coarse, dt); + EvolveB(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt) { const int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); - Real dtsdx = dt/dx[0], dtsdy = dt/dx[1], dtsdz = dt/dx[2]; + Real dtsdx = a_dt/dx[0], dtsdy = a_dt/dx[1], dtsdz = a_dt/dx[2]; MultiFab *Ex, *Ey, *Ez, *Bx, *By, *Bz; if (patch_type == PatchType::fine) @@ -164,30 +164,30 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveE (Real dt) +WarpX::EvolveE (Real a_dt) { for (int lev = 0; lev <= finest_level; ++lev) { - EvolveE(lev, dt); + EvolveE(lev, a_dt); } } void -WarpX::EvolveE (int lev, Real dt) +WarpX::EvolveE (int lev, Real a_dt) { BL_PROFILE("WarpX::EvolveE()"); - EvolveE(lev, PatchType::fine, dt); + EvolveE(lev, PatchType::fine, a_dt); if (lev > 0) { - EvolveE(lev, PatchType::coarse, dt); + EvolveE(lev, PatchType::coarse, a_dt); } } void -WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) +WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt) { - const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * dt; - const Real c2dt = (PhysConst::c*PhysConst::c) * dt; + const Real mu_c2_dt = (PhysConst::mu0*PhysConst::c*PhysConst::c) * a_dt; + const Real c2dt = (PhysConst::c*PhysConst::c) * a_dt; int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const std::array& dx = WarpX::CellSize(patch_level); @@ -358,27 +358,27 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real dt) } void -WarpX::EvolveF (Real dt, DtType dt_type) +WarpX::EvolveF (Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; for (int lev = 0; lev <= finest_level; ++lev) { - EvolveF(lev, dt, dt_type); + EvolveF(lev, a_dt, a_dt_type); } } void -WarpX::EvolveF (int lev, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; - EvolveF(lev, PatchType::fine, dt, dt_type); - if (lev > 0) EvolveF(lev, PatchType::coarse, dt, dt_type); + EvolveF(lev, PatchType::fine, a_dt, a_dt_type); + if (lev > 0) EvolveF(lev, PatchType::coarse, a_dt, a_dt_type); } void -WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) +WarpX::EvolveF (int lev, PatchType patch_type, Real a_dt, DtType a_dt_type) { if (!do_dive_cleaning) return; @@ -388,7 +388,7 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) int patch_level = (patch_type == PatchType::fine) ? lev : lev-1; const auto& dx = WarpX::CellSize(patch_level); - const std::array dtsdx {dt/dx[0], dt/dx[1], dt/dx[2]}; + const std::array dtsdx {a_dt/dx[0], a_dt/dx[1], a_dt/dx[2]}; MultiFab *Ex, *Ey, *Ez, *rho, *F; if (patch_type == PatchType::fine) @@ -408,12 +408,12 @@ WarpX::EvolveF (int lev, PatchType patch_type, Real dt, DtType dt_type) F = F_cp[lev].get(); } - const int rhocomp = (dt_type == DtType::FirstHalf) ? 0 : 1; + const int rhocomp = (a_dt_type == DtType::FirstHalf) ? 0 : 1; MultiFab src(rho->boxArray(), rho->DistributionMap(), 1, 0); ComputeDivE(src, 0, {Ex,Ey,Ez}, dx); MultiFab::Saxpy(src, -mu_c2, *rho, rhocomp, 0, 1, 0); - MultiFab::Saxpy(*F, dt, src, 0, 0, 1, 0); + MultiFab::Saxpy(*F, a_dt, src, 0, 0, 1, 0); if (do_pml && pml[lev]->ok()) { diff --git a/Source/FortranInterface/WarpX_picsar.F90 b/Source/FortranInterface/WarpX_picsar.F90 index a4cc99926..c17e8861b 100644 --- a/Source/FortranInterface/WarpX_picsar.F90 +++ b/Source/FortranInterface/WarpX_picsar.F90 @@ -265,8 +265,10 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: rho(*) real(amrex_real), intent(IN) :: rmin, dr +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 - +#endif + ! Compute the number of valid cells and guard cells integer(c_long) :: rho_nvalid(AMREX_SPACEDIM), rho_nguards(AMREX_SPACEDIM) rho_nvalid = rho_ntot - 2*rho_ng @@ -385,8 +387,9 @@ subroutine warpx_charge_deposition(rho,np,xp,yp,zp,w,q,xmin,ymin,zmin,dx,dy,dz,n real(amrex_real), intent(IN OUT):: jx(*), jy(*), jz(*) real(amrex_real), intent(IN) :: rmin, dr +#ifdef WARPX_RZ integer(c_long) :: type_rz_depose = 1 - +#endif ! Compute the number of valid cells and guard cells integer(c_long) :: jx_nvalid(AMREX_SPACEDIM), jy_nvalid(AMREX_SPACEDIM), jz_nvalid(AMREX_SPACEDIM), & jx_nguards(AMREX_SPACEDIM), jy_nguards(AMREX_SPACEDIM), jz_nguards(AMREX_SPACEDIM) diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp index 52d506bb8..2f964b6f3 100644 --- a/Source/Laser/LaserParticleContainer.cpp +++ b/Source/Laser/LaserParticleContainer.cpp @@ -78,14 +78,14 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies, parser.define(field_function); parser.registerVariables({"X","Y","t"}); - ParmParse pp("my_constants"); + ParmParse ppc("my_constants"); std::set symbols = parser.symbols(); symbols.erase("X"); symbols.erase("Y"); symbols.erase("t"); // after removing variables, we are left with constants for (auto it = symbols.begin(); it != symbols.end(); ) { Real v; - if (pp.query(it->c_str(), v)) { + if (ppc.query(it->c_str(), v)) { parser.setConstant(*it, v); it = symbols.erase(it); } else { @@ -429,8 +429,6 @@ LaserParticleContainer::Evolve (int lev, { Real wt = amrex::second(); - const Box& box = pti.validbox(); - auto& attribs = pti.GetAttribs(); auto& wp = attribs[PIdx::w ]; diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp index 348b31c8b..1cf347420 100644 --- a/Source/Parallelization/WarpXComm.cpp +++ b/Source/Parallelization/WarpXComm.cpp @@ -79,7 +79,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dBz, *Bfield_cp[lev][2], 0, 0, 1, ng); const Real* dx = Geom(lev-1).CellSize(); - const int ref_ratio = refRatio(lev-1)[0]; + const int rr = refRatio(lev-1)[0]; #ifdef _OPENMP #pragma omp parallel #endif @@ -89,7 +89,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(rr).refine(rr); // so that ccbx is coarsenable const FArrayBox& cxfab = dBx[mfi]; const FArrayBox& cyfab = dBy[mfi]; @@ -106,18 +106,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &rr,&use_limiter); #else amrex_interp_div_free_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[0]), BL_TO_FORTRAN_ANYD(bfab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - dx, &ref_ratio,&use_limiter); + dx, &rr,&use_limiter); amrex_interp_cc_bfield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(bfab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); #endif for (int idim = 0; idim < 3; ++idim) @@ -153,7 +153,7 @@ WarpX::UpdateAuxilaryData () MultiFab::Subtract(dEy, *Efield_cp[lev][1], 0, 0, 1, ng); MultiFab::Subtract(dEz, *Efield_cp[lev][2], 0, 0, 1, ng); - const int ref_ratio = refRatio(lev-1)[0]; + const int rr = refRatio(lev-1)[0]; #ifdef _OPEMP #pragma omp parallel #endif @@ -163,7 +163,7 @@ WarpX::UpdateAuxilaryData () { Box ccbx = mfi.fabbox(); ccbx.enclosedCells(); - ccbx.coarsen(ref_ratio).refine(ref_ratio); // so that ccbx is coarsenable + ccbx.coarsen(rr).refine(rr); // so that ccbx is coarsenable const FArrayBox& cxfab = dEx[mfi]; const FArrayBox& cyfab = dEy[mfi]; @@ -180,18 +180,18 @@ WarpX::UpdateAuxilaryData () BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(cyfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); #else amrex_interp_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[0]), BL_TO_FORTRAN_ANYD(efab[2]), BL_TO_FORTRAN_ANYD(cxfab), BL_TO_FORTRAN_ANYD(czfab), - &ref_ratio,&use_limiter); + &rr,&use_limiter); amrex_interp_nd_efield(ccbx.loVect(), ccbx.hiVect(), BL_TO_FORTRAN_ANYD(efab[1]), BL_TO_FORTRAN_ANYD(cyfab), - &ref_ratio); + &rr); #endif for (int idim = 0; idim < 3; ++idim) @@ -364,7 +364,7 @@ WarpX::SyncCurrent () current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& rr = refRatio(lev-1); std::array fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -372,7 +372,7 @@ WarpX::SyncCurrent () std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, rr[0]); } Vector,3> > j_fp(finest_level+1); @@ -524,10 +524,10 @@ WarpX::SyncCurrent () void WarpX::SyncCurrent (const std::array& fine, const std::array< amrex::MultiFab*,3>& crse, - int ref_ratio) + int rr) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine[0]->nGrowVect() + 1) /ref_ratio; + BL_ASSERT(rr == 2); + const IntVect& ng = (fine[0]->nGrowVect() + 1) /rr; #ifdef _OPEMP #pragma omp parallel @@ -539,7 +539,7 @@ WarpX::SyncCurrent (const std::array& fine, for (MFIter mfi(*crse[idim],true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,rr),1); ffab.resize(fbx); fbx &= (*fine[idim])[mfi].box(); ffab.setVal(0.0); @@ -564,8 +564,8 @@ WarpX::SyncRho (amrex::Vector >& rhof, for (int lev = 1; lev <= finest_level; ++lev) { rhoc[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rhof[lev], *rhoc[lev], ref_ratio[0]); + const IntVect& rr = refRatio(lev-1); + SyncRho(*rhof[lev], *rhoc[lev], rr[0]); } Vector > rho_f_g(finest_level+1); @@ -673,10 +673,10 @@ WarpX::SyncRho (amrex::Vector >& rhof, * averaging the values of the charge density of the fine patch (on the same level). */ void -WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) +WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int rr) { - BL_ASSERT(ref_ratio == 2); - const IntVect& ng = (fine.nGrowVect()+1)/ref_ratio; + BL_ASSERT(rr == 2); + const IntVect& ng = (fine.nGrowVect()+1)/rr; const int nc = fine.nComp(); #ifdef _OPEMP @@ -687,7 +687,7 @@ WarpX::SyncRho (const MultiFab& fine, MultiFab& crse, int ref_ratio) for (MFIter mfi(crse,true); mfi.isValid(); ++mfi) { const Box& bx = mfi.growntilebox(ng); - Box fbx = amrex::grow(amrex::refine(bx,ref_ratio),1); + Box fbx = amrex::grow(amrex::refine(bx,rr),1); ffab.resize(fbx, nc); fbx &= fine[mfi].box(); ffab.setVal(0.0); @@ -710,7 +710,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) current_cp[lev][1]->setVal(0.0); current_cp[lev][2]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); + const IntVect& rr = refRatio(lev-1); std::array fine { current_fp[lev][0].get(), current_fp[lev][1].get(), @@ -718,7 +718,7 @@ WarpX::RestrictCurrentFromFineToCoarsePatch (int lev) std::array< MultiFab*,3> crse { current_cp[lev][0].get(), current_cp[lev][1].get(), current_cp[lev][2].get() }; - SyncCurrent(fine, crse, ref_ratio[0]); + SyncCurrent(fine, crse, rr[0]); } void @@ -824,8 +824,8 @@ WarpX::RestrictRhoFromFineToCoarsePatch (int lev) { if (rho_fp[lev]) { rho_cp[lev]->setVal(0.0); - const IntVect& ref_ratio = refRatio(lev-1); - SyncRho(*rho_fp[lev], *rho_cp[lev], ref_ratio[0]); + const IntVect& rr = refRatio(lev-1); + SyncRho(*rho_fp[lev], *rho_cp[lev], rr[0]); } } diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 37c136a3d..212084e64 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -1024,9 +1024,6 @@ PhysicalParticleContainer::FieldGather (int lev, { const std::array& dx = WarpX::CellSize(lev); - // WarpX assumes the same number of guard cells for Ex, Ey, Ez, Bx, By, Bz - long ng = Ex.nGrow(); - BL_ASSERT(OnSameGrids(lev,Ex)); MultiFab* cost = WarpX::getCosts(lev); diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index c7ffe956b..9791eee80 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -237,10 +237,10 @@ WarpXParticleContainer::AddNParticles (int lev, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x[i]); - particle_tile.push_back_real(particle_comps["yold"], y[i]); - particle_tile.push_back_real(particle_comps["zold"], z[i]); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + ptile.push_back_real(particle_comps["xold"], x[i]); + ptile.push_back_real(particle_comps["yold"], y[i]); + ptile.push_back_real(particle_comps["zold"], z[i]); } particle_tile.push_back(p); @@ -255,10 +255,10 @@ WarpXParticleContainer::AddNParticles (int lev, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); - particle_tile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); - particle_tile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); + auto& ptile = DefineAndReturnParticleTile(0, 0, 0); + ptile.push_back_real(particle_comps["uxold"], vx + ibegin, vx + iend); + ptile.push_back_real(particle_comps["uyold"], vy + ibegin, vy + iend); + ptile.push_back_real(particle_comps["uzold"], vz + ibegin, vz + iend); } for (int comp = PIdx::uz+1; comp < PIdx::nattribs; ++comp) @@ -737,7 +737,6 @@ WarpXParticleContainer::DepositCharge (Vector >& rho, const auto& gm = m_gdb->Geom(lev); const auto& ba = m_gdb->ParticleBoxArray(lev); - const auto& dm = m_gdb->DistributionMap(lev); const Real* dx = gm.CellSize(); const Real* plo = gm.ProbLo(); @@ -807,36 +806,36 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #ifdef _OPENMP #pragma omp parallel -#endif { +#endif Cuda::ManagedDeviceVector xp, yp, zp; - FArrayBox local_rho; +#ifdef _OPENMP + FArrayBox rho_loc; +#endif for (WarpXParIter pti(*this, lev); pti.isValid(); ++pti) { - const Box& box = pti.validbox(); - auto& wp = pti.GetAttribs(PIdx::w); const long np = pti.numParticles(); pti.GetPosition(xp, yp, zp); - const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); - const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); - // Data on the grid Real* data_ptr; FArrayBox& rhofab = (*rho)[pti]; #ifdef _OPENMP + const std::array& xyzmin_tile = WarpX::LowerCorner(pti.tilebox(), lev); Box tile_box = convert(pti.tilebox(), IntVect::TheUnitVector()); const std::array& xyzmin = xyzmin_tile; tile_box.grow(ng); - local_rho.resize(tile_box); - local_rho = 0.0; - data_ptr = local_rho.dataPtr(); - auto rholen = local_rho.length(); + rho_loc.resize(tile_box); + rho_loc = 0.0; + data_ptr = rho_loc.dataPtr(); + auto rholen = rho_loc.length(); #else + const Box& box = pti.validbox(); + const std::array& xyzmin_grid = WarpX::LowerCorner(box, lev); const std::array& xyzmin = xyzmin_grid; data_ptr = rhofab.dataPtr(); auto rholen = rhofab.length(); @@ -874,10 +873,9 @@ WarpXParticleContainer::GetChargeDensity (int lev, bool local) #endif #ifdef _OPENMP - rhofab.atomicAdd(local_rho); -#endif + rhofab.atomicAdd(rho_loc); } - +#endif } if (!local) rho->SumBoundary(gm.periodicity()); diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 5f9e36bf6..a0ab1f26f 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -5,7 +5,7 @@ using namespace amrex; void -WarpX::UpdatePlasmaInjectionPosition (Real dt) +WarpX::UpdatePlasmaInjectionPosition (Real a_dt) { int dir = moving_window_dir; // Continuously inject plasma in new cells (by default only on level 0) @@ -14,12 +14,12 @@ WarpX::UpdatePlasmaInjectionPosition (Real dt) // call to this function, and injection position needs to be updated current_injection_position -= WarpX::beta_boost * #if ( AMREX_SPACEDIM == 3 ) - WarpX::boost_direction[dir] * PhysConst::c * dt; + WarpX::boost_direction[dir] * PhysConst::c * a_dt; #elif ( AMREX_SPACEDIM == 2 ) // In 2D, dir=0 corresponds to x and dir=1 corresponds to z // This needs to be converted in order to index `boost_direction` // which has 3 components, for both 2D and 3D simulations. - WarpX::boost_direction[2*dir] * PhysConst::c * dt; + WarpX::boost_direction[2*dir] * PhysConst::c * a_dt; #endif } } diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 6b6752bf1..3d7f7dcc5 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -1001,7 +1001,6 @@ WarpX::ComputeDivE (MultiFab& divE, int dcomp, void WarpX::BuildBufferMasks () { - int ngbuffer = std::max(n_field_gather_buffer, n_current_deposition_buffer); for (int lev = 1; lev <= maxLevel(); ++lev) { for (int ipass = 0; ipass < 2; ++ipass) -- cgit v1.2.3 From 46b9fde64b935703e9bda18f6150a165c14d4032 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 24 May 2019 16:19:14 -0700 Subject: option do symmetrize gaussian particle beam --- Source/Initialization/PlasmaInjector.H | 1 + Source/Initialization/PlasmaInjector.cpp | 1 + Source/Particles/PhysicalParticleContainer.H | 5 +- Source/Particles/PhysicalParticleContainer.cpp | 91 +++++++++++++++++--------- 4 files changed, 66 insertions(+), 32 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index cd5802a55..f998e217e 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -293,6 +293,7 @@ public: amrex::Real z_rms; amrex::Real q_tot; long npart; + int do_symmetrize = 0; bool radially_weighted = true; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 52b5d8b61..f9642d1b6 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -286,6 +286,7 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.get("z_rms", z_rms); pp.get("q_tot", q_tot); pp.get("npart", npart); + pp.query("do_symmetrize", do_symmetrize); gaussian_beam = true; parseMomentum(pp); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index df1e0296a..bd8cfb47e 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -94,7 +94,10 @@ public: void AddGaussianBeam(amrex::Real x_m, amrex::Real y_m, amrex::Real z_m, amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, - amrex::Real q_tot, long npart); + amrex::Real q_tot, long npart, int do_symmetrize); + + void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z, + 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, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 212084e64..38ea55d6e 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -181,7 +181,8 @@ 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); RealBox containing_bx = geom.ProbDomain(); @@ -191,13 +192,15 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution disty(y_m, y_rms); std::normal_distribution distz(z_m, z_rms); - std::array attribs; - attribs.fill(0.0); - if (ParallelDescriptor::IOProcessor()) { - std::array u; - Real weight; - for (long i = 0; i < npart; ++i) { + std::array u; + Real weight; + // 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; + } + for (long i = 0; i < npart; ++i) { #if ( AMREX_SPACEDIM == 3 | WARPX_RZ) weight = q_tot/npart/charge; Real x = distx(mt); @@ -209,35 +212,60 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real y = 0.; Real z = distz(mt); #endif - if (plasma_injector->insideBounds(x, y, z)) { - plasma_injector->getMomentum(u, x, y, z); - if (WarpX::gamma_boost > 1.) { - MapParticletoBoostedFrame(x, y, z, u); - } - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - attribs[PIdx::w ] = weight; - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); - } - - AddOneParticle(0, 0, 0, x, y, z, attribs); + if (plasma_injector->insideBounds(x, y, z)) { + plasma_injector->getMomentum(u, x, y, z); + if (do_symmetrize){ + // Add four particles to the beam: + // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) + std::array u_tmp = u; + for (int ix=0; ix<2; ix++){ + for (int iy=0; iy<2; iy++){ + x *= std::pow(-1,ix); + u_tmp[0] *= std::pow(-1,ix); + y *= std::pow(-1,iy); + u_tmp[1] *= std::pow(-1,iy); + CheckAndAddParticle(x, y, z, u_tmp, weight); + } + } + } else { + CheckAndAddParticle(x, y, z, u, weight); + } } } } Redistribute(); } +void +PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, + std::array u, + Real weight) +{ + std::array attribs; + attribs.fill(0.0); + + if (WarpX::gamma_boost > 1.) { + MapParticletoBoostedFrame(x, y, z, u); + } + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + attribs[PIdx::w ] = weight; + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } + AddOneParticle(0, 0, 0, x, y, z, attribs); +} + void PhysicalParticleContainer::AddParticles (int lev) { @@ -263,7 +291,8 @@ PhysicalParticleContainer::AddParticles (int lev) plasma_injector->y_rms, plasma_injector->z_rms, plasma_injector->q_tot, - plasma_injector->npart); + plasma_injector->npart, + plasma_injector->do_symmetrize); return; -- cgit v1.2.3 From 5473027f14805cae5f00f21a8eef2ca988d9f8ab Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Fri, 24 May 2019 16:27:09 -0700 Subject: add doc --- Docs/source/running_cpp/parameters.rst | 9 +++++++++ Source/Particles/PhysicalParticleContainer.cpp | 3 +++ 2 files changed, 12 insertions(+) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index ac378dc03..2ab071673 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -187,6 +187,15 @@ Particle initialization * ``NRandomPerCell``: injection with a fixed number of randomly-distributed particles per cell. This requires the additional parameter ``.num_particles_per_cell``. + * ``gaussian_beam``: Inject particle beam with gaussian distribution in + space in all directions. This requires additional parameters: + ``.q_tot`` (beam charge), + ``.npart`` (number of particles in the beam), + ``.x/y/z_m`` (average position in `x/y/z`), + ``.x/y/z_rms`` (standard deviation in `x/y/z`), + and optional argument ``.do_symmetrize`` (whether to + symmetrize the beam in the x and y directions). + * ``.do_continuous_injection`` (`0` or `1`) Whether to inject particles during the simulation, and not only at initialization. This can be required whith a moving window and/or when diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 38ea55d6e..bf93375fe 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -244,6 +244,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, std::array attribs; attribs.fill(0.0); + // update attribs with input arguments if (WarpX::gamma_boost > 1.) { MapParticletoBoostedFrame(x, y, z, u); } @@ -254,6 +255,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) { + // need to create old values auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); particle_tile.push_back_real(particle_comps["xold"], x); particle_tile.push_back_real(particle_comps["yold"], y); @@ -263,6 +265,7 @@ PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, particle_tile.push_back_real(particle_comps["uyold"], u[1]); particle_tile.push_back_real(particle_comps["uzold"], u[2]); } + // add particle AddOneParticle(0, 0, 0, x, y, z, attribs); } -- cgit v1.2.3 From b7dbfc6b22fb3d94310cc5f5e4b2f34fccfe49ce Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Mon, 27 May 2019 07:33:26 -0700 Subject: typo --- 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 bf93375fe..c762bdbb3 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -224,7 +224,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, u_tmp[0] *= std::pow(-1,ix); y *= std::pow(-1,iy); u_tmp[1] *= std::pow(-1,iy); - CheckAndAddParticle(x, y, z, u_tmp, weight); + CheckAndAddParticle(x, y, z, u_tmp, weight/4); } } } else { -- cgit v1.2.3 From 2ddfcb10aeb31b496153658269147e203a362657 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 28 May 2019 17:02:14 -0700 Subject: fix issue for symmetric beam --- Source/Particles/PhysicalParticleContainer.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index c762bdbb3..2c501985f 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -215,14 +215,16 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, if (plasma_injector->insideBounds(x, y, z)) { plasma_injector->getMomentum(u, x, y, z); if (do_symmetrize){ + std::array u_tmp; + Real x_tmp, y_tmp, z_tmp; // Add four particles to the beam: // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) - std::array u_tmp = u; for (int ix=0; ix<2; ix++){ for (int iy=0; iy<2; iy++){ - x *= std::pow(-1,ix); + u_tmp = u; + x_tmp = x*std::pow(-1,ix); u_tmp[0] *= std::pow(-1,ix); - y *= std::pow(-1,iy); + y_tmp = y*std::pow(-1,iy); u_tmp[1] *= std::pow(-1,iy); CheckAndAddParticle(x, y, z, u_tmp, weight/4); } -- cgit v1.2.3 From 6c20c43e0ecc9a03b21ef6c57adb2f13911004e2 Mon Sep 17 00:00:00 2001 From: MaxThevenet Date: Tue, 28 May 2019 17:19:53 -0700 Subject: typos --- Source/Particles/PhysicalParticleContainer.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'Source/Particles/PhysicalParticleContainer.cpp') diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 2c501985f..6faf7a127 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -216,7 +216,7 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, plasma_injector->getMomentum(u, x, y, z); if (do_symmetrize){ std::array u_tmp; - Real x_tmp, y_tmp, z_tmp; + Real x_tmp, y_tmp; // Add four particles to the beam: // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) for (int ix=0; ix<2; ix++){ @@ -226,7 +226,8 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, u_tmp[0] *= std::pow(-1,ix); y_tmp = y*std::pow(-1,iy); u_tmp[1] *= std::pow(-1,iy); - CheckAndAddParticle(x, y, z, u_tmp, weight/4); + CheckAndAddParticle(x_tmp, y_tmp, z, + u_tmp, weight/4); } } } else { -- cgit v1.2.3