aboutsummaryrefslogtreecommitdiff
path: root/Source/LaserParticleContainer.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'Source/LaserParticleContainer.cpp')
-rw-r--r--Source/LaserParticleContainer.cpp295
1 files changed, 269 insertions, 26 deletions
diff --git a/Source/LaserParticleContainer.cpp b/Source/LaserParticleContainer.cpp
index b8b448740..a76c2a33b 100644
--- a/Source/LaserParticleContainer.cpp
+++ b/Source/LaserParticleContainer.cpp
@@ -1,5 +1,8 @@
#include <limits>
+#include <cmath>
+#include <algorithm>
+#include <numeric>
#include <WarpX.H>
#include <WarpXConst.H>
@@ -7,11 +10,78 @@
#include <ParticleContainer.H>
#include <ParticleIterator.H>
+
+//
+// xxxxx need to make this work in 2D!
+//
+
+namespace
+{
+ Array<Real> CrossProduct (const Array<Real>& a, const Array<Real>& b)
+ {
+ return { a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0] };
+ }
+}
+
LaserParticleContainer::LaserParticleContainer (AmrCore* amr_core, int ispecies)
: WarpXParticleContainer(amr_core, ispecies)
{
- charge = PhysConst::q_e; // note that q_e is defined to be positive.
+ charge = 1.0;
mass = std::numeric_limits<Real>::max();
+
+ if (WarpX::use_laser)
+ {
+ ParmParse pp("laser");
+
+ // Parse the type of laser profile and set the corresponding flag `profile`
+ std::string laser_type_s;
+ pp.get("profile", laser_type_s);
+ std::transform(laser_type_s.begin(), laser_type_s.end(), laser_type_s.begin(), ::tolower);
+ if (laser_type_s == "gaussian") {
+ profile = laser_t::Gaussian;
+ } else {
+ BoxLib::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("e_max", e_max);
+ pp.get("wavelength", wavelength);
+
+ 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);
+ }
+
+ // 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 };
+
+ // 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);
+ if (std::abs(dp) > 1.e-14) {
+ BoxLib::Abort("Laser plane vector is not perpendicular to the main polarization vector");
+ }
+
+ p_Y = CrossProduct(nvec, p_X); // The second polarization vector
+
+#if BL_SPACEDIM == 3
+ u_X = p_X;
+ u_Y = p_Y;
+#else
+ u_X = CrossProduct({0., 1., 0.}, nvec);
+ u_Y = {0., 1., 0.};
+#endif
+ }
}
void
@@ -25,14 +95,135 @@ LaserParticleContainer::AllocData ()
void
LaserParticleContainer::InitData ()
{
-
+ m_particles.resize(GDB().finestLevel()+1);
+
+ const int lev = 0;
+
+ const Geometry& geom = GDB().Geom(lev);
+ const Real* dx = geom.CellSize();
+ const RealBox& prob_domain = geom.ProbDomain();
+
+ // spacing of laser particles in the laser plane
+ const Real eps = dx[0]*1.e-50;
+ const Real S_X = std::min(std::min(dx[0]/(std::abs(p_X[0])+eps),
+ dx[1]/(std::abs(p_X[1])+eps)),
+ dx[2]/(std::abs(p_X[2])+eps));
+ const Real S_Y = std::min(std::min(dx[0]/(std::abs(p_Y[0])+eps),
+ dx[1]/(std::abs(p_Y[1])+eps)),
+ dx[2]/(std::abs(p_Y[2])+eps));
+
+ const Real particle_weight = ComputeWeight(S_X, S_Y);
+ // The mobility is the constant of proportionality between the field to
+ // be emitted, and the corresponding velocity that the particles need to have.
+ mobility = (S_X * S_Y)/( particle_weight * PhysConst::mu0 * PhysConst::c * PhysConst::c);
+
+ // Give integer coordinates in the laser plane, return the real coordinates in the "lab" frame
+ auto Transform = [&](int i, int j) -> Array<Real>
+ {
+ return { position[0] + (S_X*i)*p_X[0] + (S_Y*j)*p_Y[0],
+ position[1] + (S_X*i)*p_X[1] + (S_Y*j)*p_Y[1],
+ position[2] + (S_X*i)*p_X[2] + (S_Y*j)*p_Y[2] };
+ };
+
+ // Given the "lab" frame coordinates, return the real coordinates in the laser plane coordinates
+ auto InverseTransform = [&](const Array<Real>& pos) -> Array<Real>
+ {
+ return { p_X[0]*pos[0] + p_X[1]*pos[1] + p_X[2]*pos[2],
+ p_Y[0]*pos[0] + p_Y[1]*pos[1] + p_Y[2]*pos[2],
+ nvec[0]*pos[0] + nvec[1]*pos[1] + nvec[2]*pos[2] };
+ };
+
+ Array<int> plane_lo(2, std::numeric_limits<int>::max());
+ Array<int> plane_hi(2, std::numeric_limits<int>::min());
+ {
+ auto compute_min_max = [&](Real x, Real y, Real z)
+ {
+ const Array<Real> pos_plane = InverseTransform({x, y, z});
+ int i = pos_plane[0]/S_X;
+ int j = pos_plane[1]/S_Y;
+ plane_lo[0] = std::min(plane_lo[0], i);
+ plane_lo[1] = std::min(plane_lo[1], j);
+ plane_hi[0] = std::max(plane_hi[0], i+1);
+ plane_hi[1] = std::max(plane_hi[1], j+1);
+ };
+
+ const Real* prob_lo = prob_domain.lo();
+ const Real* prob_hi = prob_domain.hi();
+ compute_min_max(prob_lo[0], prob_lo[1], prob_lo[2]);
+ compute_min_max(prob_hi[0], prob_lo[1], prob_lo[2]);
+ compute_min_max(prob_lo[0], prob_hi[1], prob_lo[2]);
+ compute_min_max(prob_hi[0], prob_hi[1], prob_lo[2]);
+ compute_min_max(prob_lo[0], prob_lo[1], prob_hi[2]);
+ compute_min_max(prob_hi[0], prob_lo[1], prob_hi[2]);
+ compute_min_max(prob_lo[0], prob_hi[1], prob_hi[2]);
+ compute_min_max(prob_hi[0], prob_hi[1], prob_hi[2]);
+ }
+
+ const int nprocs = ParallelDescriptor::NProcs();
+ const int myproc = ParallelDescriptor::MyProc();
+
+ const Box plane_box {IntVect(D_DECL(plane_lo[0],plane_lo[1],0)),
+ IntVect(D_DECL(plane_hi[0],plane_hi[1],0))};
+ BoxArray plane_ba {plane_box};
+ {
+ IntVect chunk(plane_box.size());
+ const int min_size = 8;
+ while (plane_ba.size() < nprocs && chunk[0] > min_size && chunk[1] > min_size)
+ {
+ for (int j = 1; j >= 0 ; j--)
+ {
+ chunk[j] /= 2;
+
+ if (plane_ba.size() < nprocs) {
+ plane_ba.maxSize(chunk);
+ }
+ }
+ }
+ }
+
+ Array<Real> particle_x, particle_y, particle_z, particle_w;
+
+ const DistributionMapping plane_dm {plane_ba, nprocs};
+ const Array<int>& procmap = plane_dm.ProcessorMap();
+ for (int i = 0, n = plane_ba.size(); i < n; ++i)
+ {
+ if (procmap[i] == myproc)
+ {
+ const Box& bx = plane_ba[i];
+ for (IntVect cell = bx.smallEnd(); cell <= bx.bigEnd(); bx.next(cell))
+ {
+ const Array<Real>& pos = Transform(cell[0], cell[1]);
+ if (prob_domain.contains(pos.data()))
+ {
+ for (int k = 0; k<2; ++k) {
+ particle_x.push_back(pos[0]);
+ particle_y.push_back(pos[1]);
+ particle_z.push_back(pos[2]);
+ }
+ particle_w.push_back( particle_weight);
+ particle_w.push_back(-particle_weight);
+ }
+ }
+ }
+ }
+ const int np = particle_z.size();
+ Array<Real> particle_ux(np, 0.0);
+ Array<Real> particle_uy(np, 0.0);
+ Array<Real> particle_uz(np, 0.0);
+
+ if (ParallelDescriptor::IOProcessor()) {
+ std::cout << "Adding laser particles\n";
+ }
+ AddNParticles(np, particle_x.data(), particle_y.data(), particle_z.data(),
+ particle_ux.data(), particle_uy.data(), particle_uz.data(),
+ 1, particle_w.data(), 1);
}
void
LaserParticleContainer::Evolve (int lev,
const MultiFab&, const MultiFab&, const MultiFab&,
const MultiFab&, const MultiFab&, const MultiFab&,
- MultiFab& jx, MultiFab& jy, MultiFab& jz, Real dt)
+ MultiFab& jx, MultiFab& jy, MultiFab& jz, Real t, Real dt)
{
BL_PROFILE("Laser::Evolve()");
BL_PROFILE_VAR_NS("Laser::Evolve::Copy", blp_copy);
@@ -57,17 +248,17 @@ LaserParticleContainer::Evolve (int lev,
long ngy_j = 0;
long ngz_j = ngx_j;
#endif
-
+
BL_ASSERT(OnSameGrids(lev,jx));
{
- Array<Real> xp, yp, zp, wp, uxp, uyp, uzp, giv;
+ Array<Real> xp, yp, zp, wp, uxp, uyp, uzp, giv,
+ plane_Xp, plane_Yp, amplitude_E;
PartIterInfo info {lev, do_tiling, tile_size};
for (PartIter pti(*this, info); pti.isValid(); ++pti)
{
const int gid = pti.index();
- const Box& tbx = pti.tilebox();
const Box& vbx = pti.validbox();
const long np = pti.numParticles();
@@ -84,6 +275,9 @@ LaserParticleContainer::Evolve (int lev,
uyp.resize(np);
uzp.resize(np);
giv.resize(np);
+ plane_Xp.resize(np);
+ plane_Yp.resize(np);
+ amplitude_E.resize(np);
//
// copy data from particle container to temp arrays
@@ -99,24 +293,36 @@ LaserParticleContainer::Evolve (int lev,
yp[i] = std::numeric_limits<Real>::quiet_NaN();
zp[i] = p.m_pos[1];
#endif
- wp[i] = p.m_data[PIdx::w];
- uxp[i] = p.m_data[PIdx::ux];
- uyp[i] = p.m_data[PIdx::uy];
- uzp[i] = p.m_data[PIdx::uz];
+ wp[i] = p.m_data[PIdx::w];
+ // Note: the particle momentum is not copied, since it is
+ // overwritten at each timestep, for the laser particles
+
+ // Find the coordinates of the particles in the emission plane
+#if (BL_SPACEDIM == 3)
+ plane_Xp[i] = u_X[0]*(xp[i] - position[0])
+ + u_X[1]*(yp[i] - position[1])
+ + u_X[2]*(zp[i] - position[2]);
+ plane_Yp[i] = u_Y[0]*(xp[i] - position[0])
+ + u_Y[1]*(yp[i] - position[1])
+ + u_Y[2]*(zp[i] - position[2]);
+#elif (BL_SPACEDIM == 2)
+ plane_Xp[i] = u_X[0]*(xp[i] - position[0])
+ + u_X[2]*(zp[i] - position[2]);
+ plane_Yp[i] = 0;
+#endif
});
BL_PROFILE_VAR_STOP(blp_copy);
-
const Box& box = BoxLib::enclosedCells(ba[gid]);
BL_ASSERT(box == vbx);
#if (BL_SPACEDIM == 3)
long nx = box.length(0);
long ny = box.length(1);
- long nz = box.length(2);
+ long nz = box.length(2);
#elif (BL_SPACEDIM == 2)
long nx = box.length(0);
long ny = 0;
- long nz = box.length(1);
+ long nz = box.length(1);
#endif
RealBox grid_box = RealBox( box, gm.CellSize(), gm.ProbLo() );
#if (BL_SPACEDIM == 3)
@@ -128,12 +334,42 @@ LaserParticleContainer::Evolve (int lev,
//
// Particle Push
//
- BL_PROFILE_VAR_START(blp_pxr_pp);
- warpx_laser_pusher(&np, xp.data(), yp.data(), zp.data(),
- uxp.data(), uyp.data(), uzp.data(), giv.data(),
- &this->charge, &dt,
- &WarpX::laser_pusher_algo);
- BL_PROFILE_VAR_STOP(blp_pxr_pp);
+ BL_PROFILE_VAR_START(blp_pxr_pp);
+ // Calculate the laser amplitude to be emitted,
+ // at the position of the emission plane
+ if (profile == laser_t::Gaussian) {
+ warpx_gaussian_laser( &np, plane_Xp.data(), plane_Yp.data(),
+ &t, &wavelength, &e_max, &profile_waist, &profile_duration,
+ &profile_t_peak, &profile_focal_distance, amplitude_E.data() );
+ }
+ // Calculate the corresponding momentum and position for the particles
+ {
+ Real v_over_c, vx, vy, vz, gamma, sign_charge;
+
+ pti.foreach([&](int i, ParticleType& p) {
+ // Calculate the velocity according to the amplitude of E
+ sign_charge = std::copysign( 1.0, wp[i] );
+ v_over_c = sign_charge * mobility * amplitude_E[i];
+ giv[i] = std::sqrt( 1 - v_over_c * v_over_c );
+ gamma = 1./giv[i];
+ // The velocity is along the laser polarization p_X
+ vx = PhysConst::c * v_over_c * p_X[0];
+ vy = PhysConst::c * v_over_c * p_X[1];
+ vz = PhysConst::c * v_over_c * p_X[2];
+ // Get the corresponding momenta
+ uxp[i] = gamma * vx;
+ uyp[i] = gamma * vy;
+ uzp[i] = gamma * vz;
+ // Push the the particle positions
+ xp[i] += vx * dt;
+#if (BL_SPACEDIM == 3)
+ yp[i] += vy * dt;
+#endif
+ zp[i] += vz * dt;
+ });
+
+ }
+ BL_PROFILE_VAR_STOP(blp_pxr_pp);
//
// Current Deposition
@@ -142,13 +378,12 @@ LaserParticleContainer::Evolve (int lev,
long lvect = 8;
BL_PROFILE_VAR_START(blp_pxr_cd);
warpx_current_deposition(jxfab.dataPtr(), jyfab.dataPtr(), jzfab.dataPtr(),
- &np, xp.data(), yp.data(), zp.data(),
- uxp.data(), uyp.data(), uzp.data(),
- giv.data(), wp.data(), &this->charge,
- &xyzmin[0], &xyzmin[1], &xyzmin[2],
+ &np, xp.data(), yp.data(), zp.data(),
+ uxp.data(), uyp.data(), uzp.data(),
+ giv.data(), wp.data(), &this->charge,
+ &xyzmin[0], &xyzmin[1], &xyzmin[2],
&dt, &dx[0], &dx[1], &dx[2], &nx, &ny, &nz,
- &ngx_j, &ngy_j, &ngz_j,
- &WarpX::nox,&WarpX::noy,&WarpX::noz,
+ &ngx_j, &ngy_j, &ngz_j, &WarpX::nox,&WarpX::noy,&WarpX::noz,
&lvect,&WarpX::current_deposition_algo);
BL_PROFILE_VAR_STOP(blp_pxr_cd);
@@ -171,6 +406,14 @@ LaserParticleContainer::Evolve (int lev,
p.m_data[PIdx::uz] = uzp[i];
});
BL_PROFILE_VAR_STOP(blp_copy);
- }
+ }
}
}
+
+Real
+LaserParticleContainer::ComputeWeight (Real Sx, Real Sy) const
+{
+ constexpr Real eps = 0.1;
+ constexpr Real fac = 1.0/(2.0*3.1415926535897932*PhysConst::mu0*PhysConst::c*PhysConst::c*eps);
+ return fac * wavelength * Sx * Sy / std::min(Sx,Sy) * e_max;
+}
generated by cgit v1.2.3 (git 2.46.0) at 2025-08-26 18:24:41 +0000