diff options
22 files changed, 293 insertions, 66 deletions
diff --git a/Docs/source/running_cpp/parameters.rst b/Docs/source/running_cpp/parameters.rst index b2810329a..3cd7af96c 100644 --- a/Docs/source/running_cpp/parameters.rst +++ b/Docs/source/running_cpp/parameters.rst @@ -454,10 +454,6 @@ Particle initialization this process (see "Lookup tables for QED modules" section below). **Implementation of this feature is in progress. It requires to compile with QED=TRUE** -* ``warpx.E_external_particle`` & ``warpx.B_external_particle`` (list of `float`) optional (default `0. 0. 0.`) - Two separate parameters which add a uniform E-field or B-field to each particle - which is then added to the field values gathered from the grid in the - PIC cycle. Laser initialization -------------------- @@ -732,6 +728,53 @@ Laser initialization the field solver. In particular, do not use any other boundary condition than periodic. +* ``particles.B_ext_particle_init_style`` (string) optional (default is "default") + This parameter determines the type of initialization for the external + magnetic field that is applied directly to the particles at every timestep. + The "default" style sets the external B-field (Bx,By,Bz) to (0.0,0.0,0.0). + The string can be set to "constant" if a constant external B-field is applied + every timestep. If this parameter is set to "constant", then an additional + parameter, namely, ``particles.B_external_particle`` must be specified in + the input file. + To parse a mathematical function for the external B-field, use the option + ``parse_B_ext_particle_function``. This option requires additional parameters + in the input file, namely, + ``particles.Bx_external_particle_function(x,y,z,t)``, + ``particles.By_external_particle_function(x,y,z,t)``, + ``particles.Bz_external_particle_function(x,y,z,t)`` to apply the external B-field + on the particles. Constants required in the mathematical expression can be set + using ``my_constants``. For a two-dimensional simulation, it is assumed that + the first and second dimensions are `x` and `z`, respectively, and the + value of the `By` component is set to zero. + Note that the current implementation of the parser for B-field on particles + does not work with RZ and the code will abort with an error message. + +* ``particles.E_ext_particle_init_style`` (string) optional (default is "default") + This parameter determines the type of initialization for the external + electric field that is applied directly to the particles at every timestep. + The "default" style set the external E-field (Ex,Ey,Ez) to (0.0,0.0,0.0). + The string can be set to "constant" if a cosntant external E-field is to be + used in the simulation at every timestep. If this parameter is set to "constant", + then an additional parameter, namely, ``particles.E_external_particle`` must be + specified in the input file. + To parse a mathematical function for the external E-field, use the option + ``parse_E_ext_particle_function``. This option requires additional + parameters in the input file, namely, + ``particles.Ex_external_particle_function(x,y,z,t)``, + ``particles.Ey_external_particle_function(x,y,z,t)``, + ``particles.Ez_external_particle_function(x,y,z,t)`` to apply the external E-field + on the particles. Constants required in the mathematical expression can be set + using ``my_constants``. For a two-dimensional simulation, similar to the B-field, + it is assumed that the first and second dimensions are `x` and `z`, respectively, + and the value of the `Ey` component is set to zero. + The current implementation of the parser for the E-field on particles does not work + with RZ and the code will abort with an error message. + +* ``particles.E_external_particle`` & ``particles.B_external_particle`` (list of `float`) optional (default `0. 0. 0.`) + Two separate parameters which add an externally applied uniform E-field or + B-field to each particle which is then added to the field values gathered + from the grid in the PIC cycle. + Collision initialization ------------------------ diff --git a/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution b/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution index 0d4cad2b1..0e1b37a00 100644 --- a/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution +++ b/Examples/Modules/qed/breit_wheeler/inputs_3d_optical_depth_evolution @@ -158,9 +158,11 @@ qed_bw.lookup_table_mode = "dummy_builtin" ################################# ### EXTERNAL E FIELD ### (3e15 * [-0.811107105653813 0.324442842261525 0.486664263392288] ) -warpx.E_external_particle = -2433321316961438 973328526784575 1459992790176863 +particles.E_ext_particle_init_style = "constant" +particles.E_external_particle = -2433321316961438 973328526784575 1459992790176863 #### ### EXTERNAL B FIELD ### (1e7 * [0.28571429 0.42857143 0.85714286] ) -warpx.B_external_particle = 2857142.85714286 4285714.28571428 8571428.57142857 +particles.B_ext_particle_init_style = "constant" +particles.B_external_particle = 2857142.85714286 4285714.28571428 8571428.57142857 #### diff --git a/Examples/Tests/Larmor/inputs_2d_mr b/Examples/Tests/Larmor/inputs_2d_mr index f933aa77a..8b197543c 100644 --- a/Examples/Tests/Larmor/inputs_2d_mr +++ b/Examples/Tests/Larmor/inputs_2d_mr @@ -34,7 +34,8 @@ geometry.prob_hi = 2.0 2.0 warpx.do_pml = 1 warpx.pml_ncell = 10 -warpx.B_external_particle = 0.0 0.00078110417851950768 0.0 +particles.B_ext_particle_init_style = "constant" +particles.B_external_particle = 0.0 0.00078110417851950768 0.0 # Verbosity warpx.verbose = 1 diff --git a/Examples/Tests/particle_pusher/inputs_3d b/Examples/Tests/particle_pusher/inputs_3d index 45ba7fa70..c61ac72c2 100644 --- a/Examples/Tests/particle_pusher/inputs_3d +++ b/Examples/Tests/particle_pusher/inputs_3d @@ -40,5 +40,7 @@ warpx.plot_raw_fields = 0 # External fields # Ex is set to be Ex = -vy*Bz -warpx.B_external_particle = 0.0 0.0 1.0 -warpx.E_external_particle = -2.994174829214179e+08 0.0 0.0 +particles.B_ext_particle_init_style = "constant" +particles.B_external_particle = 0.0 0.0 1.0 +particles.E_ext_particle_init_style = "constant" +particles.E_external_particle = -2.994174829214179e+08 0.0 0.0 diff --git a/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d b/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d index e8fdbe984..07ae2f17f 100644 --- a/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d +++ b/Examples/Tests/radiation_reaction/test_const_B_analytical/inputs_3d @@ -64,4 +64,5 @@ pos_perp2.single_particle_vel = -958.3148474999095 127.77531299998823 255.550625 pos_perp2.single_particle_weight = 1e-08 pos_perp2.do_classical_radiation_reaction = 1 -warpx.B_external_particle = 917978.2333474257 1376967.350021139 2753934.700042278 +particles.B_ext_particle_init_style = "constant" +particles.B_external_particle = 917978.2333474257 1376967.350021139 2753934.700042278 diff --git a/Source/Initialization/InjectorDensity.cpp b/Source/Initialization/InjectorDensity.cpp index 9f711a7af..fa54b342c 100644 --- a/Source/Initialization/InjectorDensity.cpp +++ b/Source/Initialization/InjectorDensity.cpp @@ -36,8 +36,8 @@ InjectorDensity::sharedMemoryNeeded () const noexcept case Type::parser: { // For parser injector, the 3D position of each particle - // is stored in shared memory. - return amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3; + // and time, t, is stored in shared memory. + return amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4; } default: return 0; diff --git a/Source/Initialization/InjectorMomentum.cpp b/Source/Initialization/InjectorMomentum.cpp index 8fadf0c4b..255883a34 100644 --- a/Source/Initialization/InjectorMomentum.cpp +++ b/Source/Initialization/InjectorMomentum.cpp @@ -30,9 +30,9 @@ InjectorMomentum::sharedMemoryNeeded () const noexcept { case Type::parser: { - // For parser injector, the 3D position of each particle + // For parser injector, the 3D position of each particle and time, t, // is stored in shared memory. - return amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3; + return amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4; } default: return 0; diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp index be29a1cbc..48c30ae93 100644 --- a/Source/Initialization/WarpXInitData.cpp +++ b/Source/Initialization/WarpXInitData.cpp @@ -535,8 +535,8 @@ WarpX::InitializeExternalFieldsOnGridUsingParser ( mfzfab(i,j,k) = zfield_parser->getField(x,y,z); }, /* To allocate shared memory for the GPU threads. */ - /* But, for now only 3 doubles (x,y,z) are allocated. */ - amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 3 + /* But, for now only 4 doubles (x,y,z,t) are allocated. */ + amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4 ); } diff --git a/Source/Parser/GpuParser.H b/Source/Parser/GpuParser.H index c158ee314..ff855d275 100644 --- a/Source/Parser/GpuParser.H +++ b/Source/Parser/GpuParser.H @@ -17,7 +17,7 @@ public: AMREX_GPU_HOST_DEVICE amrex::Real - operator() (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + operator() (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t=0.0) const noexcept { #ifdef AMREX_USE_GPU @@ -27,15 +27,17 @@ public: amrex::Gpu::SharedMemory<amrex::Real> gsm; amrex::Real* p = gsm.dataPtr(); int tid = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*(blockDim.x*blockDim.y); - p[tid*3] = x; - p[tid*3+1] = y; - p[tid*3+2] = z; + p[tid*4] = x; + p[tid*4+1] = y; + p[tid*4+2] = z; + p[tid*4+3] = t; return wp_ast_eval(m_gpu_parser.ast); #else // WarpX compiled for GPU, function compiled for __host__ m_var.x = x; m_var.y = y; m_var.z = z; + m_t = t; return wp_ast_eval(m_cpu_parser.ast); #endif @@ -49,10 +51,12 @@ public: m_var[tid].x = x; m_var[tid].y = y; m_var[tid].z = z; + m_t[tid] = t; return wp_ast_eval(m_parser[tid]->ast); #endif } + private: #ifdef AMREX_USE_GPU @@ -61,10 +65,12 @@ private: // Copy of the parser running on __host__ struct wp_parser m_cpu_parser; mutable amrex::XDim3 m_var; + mutable amrex::Real m_t; #else // Only one parser struct wp_parser** m_parser; mutable amrex::XDim3* m_var; + mutable amrex::Real* m_t; int nthreads; #endif }; diff --git a/Source/Parser/GpuParser.cpp b/Source/Parser/GpuParser.cpp index 5078b498b..ba904666b 100644 --- a/Source/Parser/GpuParser.cpp +++ b/Source/Parser/GpuParser.cpp @@ -16,6 +16,7 @@ GpuParser::GpuParser (WarpXParser const& wp) wp_parser_regvar_gpu(&m_gpu_parser, "x", 0); wp_parser_regvar_gpu(&m_gpu_parser, "y", 1); wp_parser_regvar_gpu(&m_gpu_parser, "z", 2); + wp_parser_regvar_gpu(&m_gpu_parser, "t", 3); // Initialize CPU parser: allocate memory in CUDA managed memory, // copy all data needed on CPU to m_cpu_parser @@ -28,6 +29,7 @@ GpuParser::GpuParser (WarpXParser const& wp) wp_parser_regvar(&m_cpu_parser, "x", &(m_var.x)); wp_parser_regvar(&m_cpu_parser, "y", &(m_var.y)); wp_parser_regvar(&m_cpu_parser, "z", &(m_var.z)); + wp_parser_regvar(&m_cpu_parser, "t", &(m_t)); #else // not defined AMREX_USE_GPU @@ -39,6 +41,7 @@ GpuParser::GpuParser (WarpXParser const& wp) m_parser = ::new struct wp_parser*[nthreads]; m_var = ::new amrex::XDim3[nthreads]; + m_t = ::new amrex::Real[nthreads]; for (int tid = 0; tid < nthreads; ++tid) { @@ -50,6 +53,7 @@ GpuParser::GpuParser (WarpXParser const& wp) wp_parser_regvar(m_parser[tid], "x", &(m_var[tid].x)); wp_parser_regvar(m_parser[tid], "y", &(m_var[tid].y)); wp_parser_regvar(m_parser[tid], "z", &(m_var[tid].z)); + wp_parser_regvar(m_parser[tid], "t", &(m_t[tid])); } #endif // AMREX_USE_GPU diff --git a/Source/Parser/WarpXParserWrapper.H b/Source/Parser/WarpXParserWrapper.H index 2dd7f72c7..2a4ff6fd2 100644 --- a/Source/Parser/WarpXParserWrapper.H +++ b/Source/Parser/WarpXParserWrapper.H @@ -24,7 +24,7 @@ struct ParserWrapper AMREX_GPU_HOST_DEVICE amrex::Real - getField (amrex::Real x, amrex::Real y, amrex::Real z) const noexcept + getField (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real t=0.0) const noexcept { return m_parser(x,y,z); } diff --git a/Source/Parser/wp_parser_c.h b/Source/Parser/wp_parser_c.h index 970d6b355..2cf0e2c00 100644 --- a/Source/Parser/wp_parser_c.h +++ b/Source/Parser/wp_parser_c.h @@ -30,7 +30,7 @@ wp_ast_eval (struct wp_node* node) #ifdef AMREX_DEVICE_COMPILE extern __shared__ amrex_real extern_xyz[]; int tid = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*(blockDim.x*blockDim.y); - amrex_real* x = extern_xyz + tid*3; + amrex_real* x = extern_xyz + tid*4; // parser assumes 4 independent variables (x,y,z,t) #endif switch (node->type) diff --git a/Source/Particles/MultiParticleContainer.H b/Source/Particles/MultiParticleContainer.H index 9db129b05..ed1c2f371 100644 --- a/Source/Particles/MultiParticleContainer.H +++ b/Source/Particles/MultiParticleContainer.H @@ -8,6 +8,7 @@ #include <RigidInjectedParticleContainer.H> #include <PhotonParticleContainer.H> #include <LaserParticleContainer.H> +#include <WarpXParserWrapper.H> #include <AMReX_Particles.H> #ifdef WARPX_QED @@ -215,6 +216,21 @@ public: IonizationProcess ionization_process; + std::string m_B_ext_particle_s = "default"; + std::string m_E_ext_particle_s = "default"; + // External fields added to particle fields. + amrex::Vector<amrex::Real> m_B_external_particle; + amrex::Vector<amrex::Real> m_E_external_particle; + // ParserWrapper for B_external on the particle + std::unique_ptr<ParserWrapper> m_Bx_particle_parser; + std::unique_ptr<ParserWrapper> m_By_particle_parser; + std::unique_ptr<ParserWrapper> m_Bz_particle_parser; + // ParserWrapper for E_external on the particle + std::unique_ptr<ParserWrapper> m_Ex_particle_parser; + std::unique_ptr<ParserWrapper> m_Ey_particle_parser; + std::unique_ptr<ParserWrapper> m_Ez_particle_parser; + + protected: // Particle container types diff --git a/Source/Particles/MultiParticleContainer.cpp b/Source/Particles/MultiParticleContainer.cpp index d84bc1afa..ab836ce9d 100644 --- a/Source/Particles/MultiParticleContainer.cpp +++ b/Source/Particles/MultiParticleContainer.cpp @@ -72,6 +72,93 @@ MultiParticleContainer::ReadParameters () { ParmParse pp("particles"); + // allocating and initializing default values of external fields for particles + m_E_external_particle.resize(3); + m_B_external_particle.resize(3); + // initialize E and B fields to 0.0 + for (int idim = 0; idim < 3; ++idim) { + m_E_external_particle[idim] = 0.0; + m_B_external_particle[idim] = 0.0; + } + // default values of E_external_particle and B_external_particle + // are used to set the E and B field when "constant" or "parser" + // is not explicitly used in the input + pp.query("B_ext_particle_init_style", m_B_ext_particle_s); + std::transform(m_B_ext_particle_s.begin(), + m_B_ext_particle_s.end(), + m_B_ext_particle_s.begin(), + ::tolower); + pp.query("E_ext_particle_init_style", m_E_ext_particle_s); + std::transform(m_E_ext_particle_s.begin(), + m_E_ext_particle_s.end(), + m_E_ext_particle_s.begin(), + ::tolower); + // if the input string for B_external on particles is "constant" + // then the values for the external B on particles must + // be provided in the input file. + if (m_B_ext_particle_s == "constant") + pp.getarr("B_external_particle", m_B_external_particle); + + // if the input string for E_external on particles is "constant" + // then the values for the external E on particles must + // be provided in the input file. + if (m_E_ext_particle_s == "constant") + pp.getarr("E_external_particle", m_E_external_particle); + + // if the input string for B_ext_particle_s is + // "parse_b_ext_particle_function" then the mathematical expression + // for the Bx_, By_, Bz_external_particle_function(x,y,z) + // must be provided in the input file. + if (m_B_ext_particle_s == "parse_b_ext_particle_function") { + // store the mathematical expression as string + std::string str_Bx_ext_particle_function; + std::string str_By_ext_particle_function; + std::string str_Bz_ext_particle_function; + Store_parserString(pp, "Bx_external_particle_function(x,y,z,t)", + str_Bx_ext_particle_function); + Store_parserString(pp, "By_external_particle_function(x,y,z,t)", + str_By_ext_particle_function); + Store_parserString(pp, "Bz_external_particle_function(x,y,z,t)", + str_Bz_ext_particle_function); + + // Parser for B_external on the particle + m_Bx_particle_parser.reset(new ParserWrapper( + makeParser(str_Bx_ext_particle_function))); + m_By_particle_parser.reset(new ParserWrapper( + makeParser(str_By_ext_particle_function))); + m_Bz_particle_parser.reset(new ParserWrapper( + makeParser(str_Bz_ext_particle_function))); + + } + + // if the input string for E_ext_particle_s is + // "parse_e_ext_particle_function" then the mathematical expression + // for the Ex_, Ey_, Ez_external_particle_function(x,y,z) + // must be provided in the input file. + if (m_E_ext_particle_s == "parse_e_ext_particle_function") { + // store the mathematical expression as string + std::string str_Ex_ext_particle_function; + std::string str_Ey_ext_particle_function; + std::string str_Ez_ext_particle_function; + Store_parserString(pp, "Ex_external_particle_function(x,y,z,t)", + str_Ex_ext_particle_function); + Store_parserString(pp, "Ey_external_particle_function(x,y,z,t)", + str_Ey_ext_particle_function); + Store_parserString(pp, "Ez_external_particle_function(x,y,z,t)", + str_Ez_ext_particle_function); + // Parser for E_external on the particle + m_Ex_particle_parser.reset(new ParserWrapper( + makeParser(str_Ex_ext_particle_function))); + m_Ey_particle_parser.reset(new ParserWrapper( + makeParser(str_Ey_ext_particle_function))); + m_Ez_particle_parser.reset(new ParserWrapper( + makeParser(str_Ez_ext_particle_function))); + + } + + + + pp.query("nspecies", nspecies); BL_ASSERT(nspecies >= 0); diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index 0192ffb37..74d1a0f62 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -47,6 +47,26 @@ public: amrex::Real t, amrex::Real dt) override; #endif // WARPX_DO_ELECTROSTATIC + /** + * \brief Apply external E and B fields on the particles. The E and B + * fields could be defined as a constant or using a parser for reading + * in a mathematical expression. The default value for the E- and B-fields + * is (0.0,0.0,0.0). + * + * \param[in,out] Exp-Bzp pointer to fields on particles modified based + * on external E and B + * \param[in] xp,yp,zp arrays of particle positions required to compute + * mathematical expression for the external fields + * using parser. + */ + void AssignExternalFieldOnParticles ( WarpXParIter& pti, + RealVector& Exp, RealVector& Eyp, + RealVector& Ezp, RealVector& Bxp, + RealVector& Byp, RealVector& Bzp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal> xp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal> yp, + amrex::Gpu::ManagedDeviceVector<amrex::ParticleReal> zp, int lev); + virtual void FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -318,6 +338,7 @@ protected: //radiation reaction bool do_classical_radiation_reaction = false; + #ifdef WARPX_QED // A flag to enable quantum_synchrotron process for leptons bool m_do_qed_quantum_sync = false; diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index 936a64c9c..21c8ddd31 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -967,6 +967,80 @@ PhysicalParticleContainer::EvolveES (const Vector<std::array<std::unique_ptr<Mul #endif // WARPX_DO_ELECTROSTATIC void +PhysicalParticleContainer::AssignExternalFieldOnParticles(WarpXParIter& pti, + RealVector& Exp, RealVector& Eyp, RealVector& Ezp, + RealVector& Bxp, RealVector& Byp, RealVector& Bzp, + Gpu::ManagedDeviceVector<ParticleReal> xp, + Gpu::ManagedDeviceVector<ParticleReal> yp, + Gpu::ManagedDeviceVector<ParticleReal> zp, int lev) +{ + const long np = pti.numParticles(); + /// get WarpX class object + auto & warpx = WarpX::GetInstance(); + /// get MultiParticleContainer class object + auto & mypc = warpx.GetPartContainer(); + if (mypc.m_E_ext_particle_s=="constant" || + mypc.m_E_ext_particle_s=="default") { + Exp.assign(np,mypc.m_E_external_particle[0]); + Eyp.assign(np,mypc.m_E_external_particle[1]); + Ezp.assign(np,mypc.m_E_external_particle[2]); + } + if (mypc.m_B_ext_particle_s=="constant" || + mypc.m_B_ext_particle_s=="default") { + Bxp.assign(np,mypc.m_B_external_particle[0]); + Byp.assign(np,mypc.m_B_external_particle[1]); + Bzp.assign(np,mypc.m_B_external_particle[2]); + } + if (mypc.m_E_ext_particle_s=="parse_e_ext_particle_function") { + Real* const AMREX_RESTRICT xp_data = xp.dataPtr(); + Real* const AMREX_RESTRICT yp_data = yp.dataPtr(); + Real* const AMREX_RESTRICT zp_data = zp.dataPtr(); + Real* const AMREX_RESTRICT Exp_data = Exp.dataPtr(); + Real* const AMREX_RESTRICT Eyp_data = Eyp.dataPtr(); + Real* const AMREX_RESTRICT Ezp_data = Ezp.dataPtr(); + ParserWrapper *xfield_partparser = mypc.m_Ex_particle_parser.get(); + ParserWrapper *yfield_partparser = mypc.m_Ey_particle_parser.get(); + ParserWrapper *zfield_partparser = mypc.m_Ez_particle_parser.get(); + Real time = warpx.gett_new(lev); + amrex::ParallelFor(pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + Exp_data[i] = xfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); + Eyp_data[i] = yfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); + Ezp_data[i] = zfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); + }, + /* To allocate shared memory for the GPU threads. */ + /* But, for now only 4 doubles (x,y,z,t) are allocated. */ + amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4 + ); + } + if (mypc.m_B_ext_particle_s=="parse_b_ext_particle_function") { + Real* const AMREX_RESTRICT xp_data = xp.dataPtr(); + Real* const AMREX_RESTRICT yp_data = yp.dataPtr(); + Real* const AMREX_RESTRICT zp_data = zp.dataPtr(); + Real* const AMREX_RESTRICT Bxp_data = Bxp.dataPtr(); + Real* const AMREX_RESTRICT Byp_data = Byp.dataPtr(); + Real* const AMREX_RESTRICT Bzp_data = Bzp.dataPtr(); + ParserWrapper *xfield_partparser = mypc.m_Bx_particle_parser.get(); + ParserWrapper *yfield_partparser = mypc.m_By_particle_parser.get(); + ParserWrapper *zfield_partparser = mypc.m_Bz_particle_parser.get(); + Real time = warpx.gett_new(lev); + amrex::ParallelFor(pti.numParticles(), + [=] AMREX_GPU_DEVICE (long i) { + Bxp_data[i] = xfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); + Byp_data[i] = yfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); + Bzp_data[i] = zfield_partparser->getField(xp_data[i],yp_data[i],zp_data[i],time); + }, + /* To allocate shared memory for the GPU threads. */ + /* But, for now only 4 doubles (x,y,z,t) are allocated. */ + amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double) * 4 + ); + } + +} + + + +void PhysicalParticleContainer::FieldGather (int lev, const amrex::MultiFab& Ex, const amrex::MultiFab& Ey, @@ -1015,13 +1089,6 @@ PhysicalParticleContainer::FieldGather (int lev, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,WarpX::E_external_particle[0]); - Eyp.assign(np,WarpX::E_external_particle[1]); - Ezp.assign(np,WarpX::E_external_particle[2]); - Bxp.assign(np,WarpX::B_external_particle[0]); - Byp.assign(np,WarpX::B_external_particle[1]); - Bzp.assign(np,WarpX::B_external_particle[2]); - // // copy data from particle container to temp arrays // @@ -1152,14 +1219,6 @@ PhysicalParticleContainer::Evolve (int lev, exfab, eyfab, ezfab, bxfab, byfab, bzfab); } - Exp.assign(np,WarpX::E_external_particle[0]); - Eyp.assign(np,WarpX::E_external_particle[1]); - Ezp.assign(np,WarpX::E_external_particle[2]); - - Bxp.assign(np,WarpX::B_external_particle[0]); - Byp.assign(np,WarpX::B_external_particle[1]); - Bzp.assign(np,WarpX::B_external_particle[2]); - // Determine which particles deposit/gather in the buffer, and // which particles deposit/gather in the fine patch long nfine_current = np; @@ -1799,14 +1858,6 @@ PhysicalParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,WarpX::E_external_particle[0]); - Eyp.assign(np,WarpX::E_external_particle[1]); - Ezp.assign(np,WarpX::E_external_particle[2]); - - Bxp.assign(np,WarpX::B_external_particle[0]); - Byp.assign(np,WarpX::B_external_particle[1]); - Bzp.assign(np,WarpX::B_external_particle[2]); - // // copy data from particle container to temp arrays // @@ -2176,9 +2227,16 @@ PhysicalParticleContainer::FieldGather (WarpXParIter& pti, AMREX_ALWAYS_ASSERT_WITH_MESSAGE((gather_lev==(lev-1)) || (gather_lev==(lev )), "Gather buffers only work for lev-1"); - // If no particles, do not do anything if (np_to_gather == 0) return; + + // initializing the field value to the externally applied field before + // gathering fields from the grid to the particles. + AssignExternalFieldOnParticles(pti, Exp, Eyp, Ezp, Bxp, Byp, Bzp, + m_xp[thread_num], m_yp[thread_num], + m_zp[thread_num], lev); + + // Get cell size on gather_lev const std::array<Real,3>& dx = WarpX::CellSize(std::max(gather_lev,0)); @@ -2424,4 +2482,5 @@ set_quantum_sync_engine_ptr(std::shared_ptr<QuantumSynchrotronEngine> ptr) { m_shr_p_qs_engine = ptr; } + #endif diff --git a/Source/Particles/RigidInjectedParticleContainer.cpp b/Source/Particles/RigidInjectedParticleContainer.cpp index bee71fba1..f3b502a0a 100644 --- a/Source/Particles/RigidInjectedParticleContainer.cpp +++ b/Source/Particles/RigidInjectedParticleContainer.cpp @@ -392,13 +392,6 @@ RigidInjectedParticleContainer::PushP (int lev, Real dt, const FArrayBox& byfab = By[pti]; const FArrayBox& bzfab = Bz[pti]; - Exp.assign(np,WarpX::E_external_particle[0]); - Eyp.assign(np,WarpX::E_external_particle[1]); - Ezp.assign(np,WarpX::E_external_particle[2]); - Bxp.assign(np,WarpX::B_external_particle[0]); - Byp.assign(np,WarpX::B_external_particle[1]); - Bzp.assign(np,WarpX::B_external_particle[2]); - // // copy data from particle container to temp arrays // diff --git a/Source/Particles/WarpXParticleContainer.cpp b/Source/Particles/WarpXParticleContainer.cpp index 6abc02139..85c4dc0d5 100644 --- a/Source/Particles/WarpXParticleContainer.cpp +++ b/Source/Particles/WarpXParticleContainer.cpp @@ -113,6 +113,7 @@ WarpXParticleContainer::WarpXParticleContainer (AmrCore* amr_core, int ispecies) m_xp.resize(num_threads); m_yp.resize(num_threads); m_zp.resize(num_threads); + } void diff --git a/Source/Utils/WarpXMovingWindow.cpp b/Source/Utils/WarpXMovingWindow.cpp index 0fcbb6a88..3f607615b 100644 --- a/Source/Utils/WarpXMovingWindow.cpp +++ b/Source/Utils/WarpXMovingWindow.cpp @@ -323,7 +323,7 @@ WarpX::shiftMF (MultiFab& mf, const Geometry& geom, int num_shift, int dir, #endif srcfab(i,j,k,n) = field_parser->getField(x,y,z); } - , amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double)*3 + , amrex::Gpu::numThreadsPerBlockParallelFor() * sizeof(double)*4 ); } diff --git a/Source/Utils/WarpXUtil.cpp b/Source/Utils/WarpXUtil.cpp index e9fb958fd..a154e93df 100644 --- a/Source/Utils/WarpXUtil.cpp +++ b/Source/Utils/WarpXUtil.cpp @@ -185,12 +185,13 @@ void Store_parserString(amrex::ParmParse& pp, std::string query_string, WarpXParser makeParser (std::string const& parse_function) { WarpXParser parser(parse_function); - parser.registerVariables({"x","y","z"}); + parser.registerVariables({"x","y","z","t"}); ParmParse pp("my_constants"); std::set<std::string> symbols = parser.symbols(); symbols.erase("x"); symbols.erase("y"); symbols.erase("z"); + symbols.erase("t"); for (auto it = symbols.begin(); it != symbols.end(); ) { Real v; if (pp.query(it->c_str(), v)) { diff --git a/Source/WarpX.H b/Source/WarpX.H index 2c47c193f..1549dded2 100644 --- a/Source/WarpX.H +++ b/Source/WarpX.H @@ -86,15 +86,11 @@ public: //! Author of an input file / simulation setup static std::string authors; - // External fields added to particle fields. - static amrex::Vector<amrex::Real> B_external_particle; - static amrex::Vector<amrex::Real> E_external_particle; - // Initial field on the grid. static amrex::Vector<amrex::Real> E_external_grid; static amrex::Vector<amrex::Real> B_external_grid; - // Initialization Type for External E and B + // Initialization Type for External E and B on grid static std::string B_ext_grid_s; static std::string E_ext_grid_s; diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp index 731e9914b..5f6f2b2e5 100644 --- a/Source/WarpX.cpp +++ b/Source/WarpX.cpp @@ -24,9 +24,6 @@ using namespace amrex; -Vector<Real> WarpX::B_external_particle(3, 0.0); -Vector<Real> WarpX::E_external_particle(3, 0.0); - Vector<Real> WarpX::E_external_grid(3, 0.0); Vector<Real> WarpX::B_external_grid(3, 0.0); @@ -342,9 +339,6 @@ WarpX::ReadParameters () pp.query("zmax_plasma_to_compute_max_step", zmax_plasma_to_compute_max_step); - pp.queryarr("B_external_particle", B_external_particle); - pp.queryarr("E_external_particle", E_external_particle); - pp.query("do_moving_window", do_moving_window); if (do_moving_window) { |