aboutsummaryrefslogtreecommitdiff
path: root/Source/Laser/LaserParticleContainer.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/Laser/LaserParticleContainer.cpp')
-rw-r--r--Source/Laser/LaserParticleContainer.cpp245
1 files changed, 138 insertions, 107 deletions
diff --git a/Source/Laser/LaserParticleContainer.cpp b/Source/Laser/LaserParticleContainer.cpp
index dca73de50..8571c74ad 100644
--- a/Source/Laser/LaserParticleContainer.cpp
+++ b/Source/Laser/LaserParticleContainer.cpp
@@ -40,115 +40,116 @@ LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies,
profile = laser_t::Harris;
} else if(laser_type_s == "parse_field_function") {
profile = laser_t::parse_field_function;
- } else {
- amrex::Abort("Unknown laser type");
- }
+ } else {
+ amrex::Abort("Unknown laser type");
+ }
- // Parse the properties of the antenna
- pp.getarr("position", position);
- pp.getarr("direction", nvec);
- pp.getarr("polarization", p_X);
- pp.query("pusher_algo", pusher_algo);
- pp.get("wavelength", wavelength);
- pp.get("e_max", e_max);
- pp.query("do_continuous_injection", do_continuous_injection);
-
- if ( profile == laser_t::Gaussian ) {
- // Parse the properties of the Gaussian profile
- pp.get("profile_waist", profile_waist);
- pp.get("profile_duration", profile_duration);
- pp.get("profile_t_peak", profile_t_peak);
- pp.get("profile_focal_distance", profile_focal_distance);
- stc_direction = p_X;
- pp.queryarr("stc_direction", stc_direction);
- pp.query("zeta", zeta);
- pp.query("beta", beta);
- pp.query("phi2", phi2);
- }
+ // Parse the properties of the antenna
+ pp.getarr("position", position);
+ pp.getarr("direction", nvec);
+ pp.getarr("polarization", p_X);
+ pp.query("pusher_algo", pusher_algo);
+ pp.get("wavelength", wavelength);
+ pp.get("e_max", e_max);
+ pp.query("do_continuous_injection", do_continuous_injection);
+ pp.query("min_particles_per_mode", min_particles_per_mode);
+
+ if ( profile == laser_t::Gaussian ) {
+ // Parse the properties of the Gaussian profile
+ pp.get("profile_waist", profile_waist);
+ pp.get("profile_duration", profile_duration);
+ pp.get("profile_t_peak", profile_t_peak);
+ pp.get("profile_focal_distance", profile_focal_distance);
+ stc_direction = p_X;
+ pp.queryarr("stc_direction", stc_direction);
+ pp.query("zeta", zeta);
+ pp.query("beta", beta);
+ pp.query("phi2", phi2);
+ }
- if ( profile == laser_t::Harris ) {
- // Parse the properties of the Harris profile
- pp.get("profile_waist", profile_waist);
- pp.get("profile_duration", profile_duration);
- pp.get("profile_focal_distance", profile_focal_distance);
- }
+ if ( profile == laser_t::Harris ) {
+ // Parse the properties of the Harris profile
+ pp.get("profile_waist", profile_waist);
+ pp.get("profile_duration", profile_duration);
+ pp.get("profile_focal_distance", profile_focal_distance);
+ }
- if ( profile == laser_t::parse_field_function ) {
- // Parse the properties of the parse_field_function profile
- pp.get("field_function(X,Y,t)", field_function);
- parser.define(field_function);
- parser.registerVariables({"X","Y","t"});
-
- ParmParse ppc("my_constants");
- std::set<std::string> 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 (ppc.query(it->c_str(), v)) {
- parser.setConstant(*it, v);
- it = symbols.erase(it);
- } else {
- ++it;
- }
- }
- for (auto const& s : symbols) { // make sure there no unknown symbols
- amrex::Abort("Laser Profile: Unknown symbol "+s);
+ if ( profile == laser_t::parse_field_function ) {
+ // Parse the properties of the parse_field_function profile
+ pp.get("field_function(X,Y,t)", field_function);
+ parser.define(field_function);
+ parser.registerVariables({"X","Y","t"});
+
+ ParmParse ppc("my_constants");
+ std::set<std::string> 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 (ppc.query(it->c_str(), v)) {
+ parser.setConstant(*it, v);
+ it = symbols.erase(it);
+ } else {
+ ++it;
}
}
-
- // Plane normal
- Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]);
- nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s };
-
- if (WarpX::gamma_boost > 1.) {
- // Check that the laser direction is equal to the boost direction
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nvec[0]*WarpX::boost_direction[0]
- + nvec[1]*WarpX::boost_direction[1]
- + nvec[2]*WarpX::boost_direction[2] - 1. < 1.e-12,
- "The Lorentz boost should be in the same direction as the laser propagation");
- // Get the position of the plane, along the boost direction, in the lab frame
- // and convert the position of the antenna to the boosted frame
- Z0_lab = nvec[0]*position[0] + nvec[1]*position[1] + nvec[2]*position[2];
- Real Z0_boost = Z0_lab/WarpX::gamma_boost;
- position[0] += (Z0_boost-Z0_lab)*nvec[0];
- position[1] += (Z0_boost-Z0_lab)*nvec[1];
- position[2] += (Z0_boost-Z0_lab)*nvec[2];
+ for (auto const& s : symbols) { // make sure there no unknown symbols
+ amrex::Abort("Laser Profile: Unknown symbol "+s);
}
+ }
- // The first polarization vector
- s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]);
- p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s };
+ // Plane normal
+ Real s = 1.0/std::sqrt(nvec[0]*nvec[0] + nvec[1]*nvec[1] + nvec[2]*nvec[2]);
+ nvec = { nvec[0]*s, nvec[1]*s, nvec[2]*s };
+
+ if (WarpX::gamma_boost > 1.) {
+ // Check that the laser direction is equal to the boost direction
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE( nvec[0]*WarpX::boost_direction[0]
+ + nvec[1]*WarpX::boost_direction[1]
+ + nvec[2]*WarpX::boost_direction[2] - 1. < 1.e-12,
+ "The Lorentz boost should be in the same direction as the laser propagation");
+ // Get the position of the plane, along the boost direction, in the lab frame
+ // and convert the position of the antenna to the boosted frame
+ Z0_lab = nvec[0]*position[0] + nvec[1]*position[1] + nvec[2]*position[2];
+ Real Z0_boost = Z0_lab/WarpX::gamma_boost;
+ position[0] += (Z0_boost-Z0_lab)*nvec[0];
+ position[1] += (Z0_boost-Z0_lab)*nvec[1];
+ position[2] += (Z0_boost-Z0_lab)*nvec[2];
+ }
+
+ // The first polarization vector
+ s = 1.0/std::sqrt(p_X[0]*p_X[0] + p_X[1]*p_X[1] + p_X[2]*p_X[2]);
+ p_X = { p_X[0]*s, p_X[1]*s, p_X[2]*s };
- Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0);
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
- "Laser plane vector is not perpendicular to the main polarization vector");
+ Real dp = std::inner_product(nvec.begin(), nvec.end(), p_X.begin(), 0.0);
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
+ "Laser plane vector is not perpendicular to the main polarization vector");
- p_Y = CrossProduct(nvec, p_X); // The second polarization vector
+ p_Y = CrossProduct(nvec, p_X); // The second polarization vector
- s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]);
- stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s };
- dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0);
- AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
- "stc_direction is not perpendicular to the laser plane vector");
+ s = 1.0/std::sqrt(stc_direction[0]*stc_direction[0] + stc_direction[1]*stc_direction[1] + stc_direction[2]*stc_direction[2]);
+ stc_direction = { stc_direction[0]*s, stc_direction[1]*s, stc_direction[2]*s };
+ dp = std::inner_product(nvec.begin(), nvec.end(), stc_direction.begin(), 0.0);
+ AMREX_ALWAYS_ASSERT_WITH_MESSAGE(std::abs(dp) < 1.0e-14,
+ "stc_direction is not perpendicular to the laser plane vector");
- // Get angle between p_X and stc_direction
- // in 2d, stcs are in the simulation plane
+ // Get angle between p_X and stc_direction
+ // in 2d, stcs are in the simulation plane
#if AMREX_SPACEDIM == 3
- theta_stc = acos(stc_direction[0]*p_X[0] +
+ theta_stc = acos(stc_direction[0]*p_X[0] +
stc_direction[1]*p_X[1] +
stc_direction[2]*p_X[2]);
#else
- theta_stc = 0.;
+ theta_stc = 0.;
#endif
-#if AMREX_SPACEDIM == 3
- u_X = p_X;
- u_Y = p_Y;
+#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
+ u_X = p_X;
+ u_Y = p_Y;
#else
- u_X = CrossProduct({0., 1., 0.}, nvec);
- u_Y = {0., 1., 0.};
+ u_X = CrossProduct({0., 1., 0.}, nvec);
+ u_Y = {0., 1., 0.};
#endif
laser_injection_box= Geom(0).ProbDomain();
@@ -271,9 +272,15 @@ LaserParticleContainer::InitData (int lev)
position[1] + (S_X*(i+0.5))*u_X[1] + (S_Y*(j+0.5))*u_Y[1],
position[2] + (S_X*(i+0.5))*u_X[2] + (S_Y*(j+0.5))*u_Y[2] };
#else
+# if (defined WARPX_DIM_RZ)
+ return { position[0] + (S_X*(i+0.5)),
+ 0.0,
+ position[2]};
+# else
return { position[0] + (S_X*(i+0.5))*u_X[0],
- 0.0,
- position[2] + (S_X*(i+0.5))*u_X[2] };
+ 0.0,
+ position[2] + (S_X*(i+0.5))*u_X[2] };
+# endif
#endif
};
@@ -283,7 +290,11 @@ LaserParticleContainer::InitData (int lev)
return {u_X[0]*(pos[0]-position[0])+u_X[1]*(pos[1]-position[1])+u_X[2]*(pos[2]-position[2]),
u_Y[0]*(pos[0]-position[0])+u_Y[1]*(pos[1]-position[1])+u_Y[2]*(pos[2]-position[2])};
#else
+# if (defined WARPX_DIM_RZ)
+ return {pos[0]-position[0], 0.0};
+# else
return {u_X[0]*(pos[0]-position[0])+u_X[2]*(pos[2]-position[2]), 0.0};
+# endif
#endif
};
@@ -364,6 +375,7 @@ LaserParticleContainer::InitData (int lev)
#endif
if (laser_injection_box.contains(x))
{
+#ifndef WARPX_DIM_RZ
for (int k = 0; k<2; ++k) {
particle_x.push_back(pos[0]);
particle_y.push_back(pos[1]);
@@ -371,6 +383,21 @@ LaserParticleContainer::InitData (int lev)
}
particle_w.push_back( weight);
particle_w.push_back(-weight);
+#else
+ // Particles are laid out in radial spokes
+ const int n_spokes = (WarpX::n_rz_azimuthal_modes - 1)*min_particles_per_mode;
+ for (int spoke = 0 ; spoke < n_spokes ; spoke++) {
+ const Real phase = 2.*MathConst::pi*spoke/n_spokes;
+ for (int k = 0; k<2; ++k) {
+ particle_x.push_back(pos[0]*std::cos(phase));
+ particle_y.push_back(pos[0]*std::sin(phase));
+ particle_z.push_back(pos[2]);
+ }
+ const Real r_weight = weight*2.*MathConst::pi*pos[0]/n_spokes;
+ particle_w.push_back( r_weight);
+ particle_w.push_back(-r_weight);
+ }
+#endif
}
}
}
@@ -564,13 +591,17 @@ LaserParticleContainer::ComputeSpacing (int lev, Real& Sx, Real& Sy) const
#if (AMREX_SPACEDIM == 3)
Sx = std::min(std::min(dx[0]/(std::abs(u_X[0])+eps),
dx[1]/(std::abs(u_X[1])+eps)),
- dx[2]/(std::abs(u_X[2])+eps));
+ dx[2]/(std::abs(u_X[2])+eps));
Sy = std::min(std::min(dx[0]/(std::abs(u_Y[0])+eps),
dx[1]/(std::abs(u_Y[1])+eps)),
- dx[2]/(std::abs(u_Y[2])+eps));
+ dx[2]/(std::abs(u_Y[2])+eps));
#else
+# if (defined WARPX_DIM_RZ)
+ Sx = dx[0];
+# else
Sx = std::min(dx[0]/(std::abs(u_X[0])+eps),
dx[2]/(std::abs(u_X[2])+eps));
+# endif
Sy = 1.0;
#endif
}
@@ -579,7 +610,7 @@ void
LaserParticleContainer::ComputeWeightMobility (Real Sx, Real Sy)
{
constexpr Real eps = 0.01;
- constexpr Real fac = 1.0/(2.0*3.1415926535897932*PhysConst::mu0*PhysConst::c*PhysConst::c*eps);
+ constexpr Real fac = 1.0/(2.0*MathConst::pi*PhysConst::mu0*PhysConst::c*PhysConst::c*eps);
weight = fac * wavelength * Sx * Sy / std::min(Sx,Sy) * e_max;
// The mobility is the constant of proportionality between the field to
@@ -612,14 +643,14 @@ LaserParticleContainer::calculate_laser_plane_coordinates (
Real * AMREX_RESTRICT const pplane_Xp,
Real * AMREX_RESTRICT const pplane_Yp)
{
- Real const * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr();
- Real const * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr();
- Real const * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr();
+ ParticleReal const * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr();
+ ParticleReal const * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr();
+ ParticleReal const * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr();
Real tmp_u_X_0 = u_X[0];
Real tmp_u_X_2 = u_X[2];
Real tmp_position_0 = position[0];
Real tmp_position_2 = position[2];
-#if (AMREX_SPACEDIM == 3)
+#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
Real tmp_u_X_1 = u_X[1];
Real tmp_u_Y_0 = u_Y[0];
Real tmp_u_Y_1 = u_Y[1];
@@ -630,7 +661,7 @@ LaserParticleContainer::calculate_laser_plane_coordinates (
amrex::ParallelFor(
np,
[=] AMREX_GPU_DEVICE (int i) {
-#if (AMREX_SPACEDIM == 3)
+#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
pplane_Xp[i] =
tmp_u_X_0 * (xp[i] - tmp_position_0) +
tmp_u_X_1 * (yp[i] - tmp_position_1) +
@@ -660,14 +691,14 @@ LaserParticleContainer::calculate_laser_plane_coordinates (
*/
void
LaserParticleContainer::update_laser_particle(
- const int np, Real * AMREX_RESTRICT const puxp, Real * AMREX_RESTRICT const puyp,
- Real * AMREX_RESTRICT const puzp, Real const * AMREX_RESTRICT const pwp,
+ const int np, ParticleReal * AMREX_RESTRICT const puxp, ParticleReal * AMREX_RESTRICT const puyp,
+ ParticleReal * AMREX_RESTRICT const puzp, ParticleReal const * AMREX_RESTRICT const pwp,
Real const * AMREX_RESTRICT const amplitude, const Real dt,
const int thread_num)
{
- Real * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr();
- Real * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr();
- Real * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr();
+ ParticleReal * const AMREX_RESTRICT xp = m_xp[thread_num].dataPtr();
+ ParticleReal * const AMREX_RESTRICT yp = m_yp[thread_num].dataPtr();
+ ParticleReal * const AMREX_RESTRICT zp = m_zp[thread_num].dataPtr();
Real tmp_p_X_0 = p_X[0];
Real tmp_p_X_1 = p_X[1];
Real tmp_p_X_2 = p_X[2];
@@ -702,7 +733,7 @@ LaserParticleContainer::update_laser_particle(
puzp[i] = gamma * vz;
// Push the the particle positions
xp[i] += vx * dt;
-#if (AMREX_SPACEDIM == 3)
+#if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
yp[i] += vy * dt;
#endif
zp[i] += vz * dt;
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-17 14:14:50 +0000