Skip to content

Commit 7bb4c8a

Browse files
RemiLeheAndrew Myers
RemiLehe
authored and
Andrew Myers
committed
Remove superfluous if condition in impact ionization (BLAST-WarpX#5895)
This `if` condition is already enclosed in an existing `if` condition: https://github.com/BLAST-WarpX/warpx/blob/development/Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H#L322
1 parent e9c7eb3 commit 7bb4c8a

File tree

1 file changed

+87
-92
lines changed

1 file changed

+87
-92
lines changed

‎Source/Particles/Collision/BinaryCollision/DSMC/SplitAndScatterFunc.H‎

Lines changed: 87 additions & 92 deletions
Original file line numberDiff line numberDiff line change
@@ -392,99 +392,94 @@ public:
392392
uy_p2 -= uCOM_y;
393393
uz_p2 -= uCOM_z;
394394

395-
if (mask[i] == int(ScatteringProcessType::IONIZATION)) {
396-
// calculate kinetic energy of the collision (in eV)
397-
const amrex::ParticleReal E1 = (
398-
0.5_prt * m1 * (ux1*ux1 + uy1*uy1 + uz1*uz1) / PhysConst::q_e
399-
);
400-
const amrex::ParticleReal E2 = (
401-
0.5_prt * m2 * (ux_p2*ux_p2 + uy_p2*uy_p2 + uz_p2*uz_p2) / PhysConst::q_e
402-
);
403-
const amrex::ParticleReal E_coll = E1 + E2;
404-
405-
// subtract the energy cost for ionization
406-
const amrex::ParticleReal E_out = (E_coll - ionization_energy) * PhysConst::q_e;
407-
408-
// Momentum and energy division after the ionization event
409-
// is done as follows:
410-
// Three numbers are generated that satisfy the triangle
411-
// inequality. These numbers will be scaled to give the
412-
// momentum for each particle (satisfying the triangle
413-
// inequality ensures that the momentum vectors can be
414-
// arranged such that the net linear momentum is zero).
415-
// Product 2 is definitely an ion, so we first choose its
416-
// proportion to be n_2 = min(E2 / E_coll, 0.5 * R), where
417-
// R is a random number between 0 and 1.
418-
// The other two numbers (n_0 & n_1) are then found
419-
// by choosing a random point on the ellipse characterized
420-
// by the semi-major and semi-minor axes
421-
// a = (1 - n_0) / 2.0
422-
// b = 0.5 * sqrt(1 - 2 * n_0)
423-
// The numbers are found by randomly sampling an x value
424-
// between -a and a, and finding the corresponding y value
425-
// that falls on the ellipse: y^2 = b^2 - b^2/a^2 * x^2.
426-
// Then n_0 = sqrt(y^2 + (x - n_2/2)^2) and
427-
// n_1 = 1 - n_2 - n_0.
428-
// Next, we need to find the number C, such that if
429-
// p_i = C * n_i, we would have E_0 + E_1 + E_2 = E_out
430-
// where E_i = p_i^2 / (2 M_i), i.e.,
431-
// C^2 * \sum_i [n_i^2 / (2 M_i)] = E_out
432-
// After C is determined the momentum vectors are arranged
433-
// such that the net linear momentum is zero.
434-
435-
const amrex::ParticleReal n_2 = std::min(E2 / E_coll, 0.5_prt * amrex::Random(engine));
436-
437-
// find ellipse semi-major and minor axis
438-
const amrex::ParticleReal a = 0.5_prt * (1.0_prt - n_2);
439-
const amrex::ParticleReal b = 0.5_prt * std::sqrt(1.0_prt - 2.0_prt * n_2);
440-
441-
// sample random x value and calculate y
442-
const amrex::ParticleReal x = (2._prt * amrex::Random(engine) - 1.0_prt) * a;
443-
const amrex::ParticleReal y2 = b*b - b*b/(a*a) * x*x;
444-
const amrex::ParticleReal n_0 = std::sqrt(y2 + x*x - x*n_2 + 0.25_prt*n_2*n_2);
445-
const amrex::ParticleReal n_1 = 1.0_prt - n_0 - n_2;
446-
447-
// calculate the value of C
448-
const amrex::ParticleReal C = std::sqrt(E_out / (
449-
n_0*n_0 / (2.0_prt * m1) + n_1*n_1 / (2.0_prt * m_ioniz_product1) + n_2*n_2 / (2.0_prt * m_ioniz_product2)
450-
));
451-
452-
// Now that appropriate momenta are set for each outgoing species
453-
// the directions for the velocity vectors must be chosen such
454-
// that the net linear momentum in the current frame is 0.
455-
// This is achieved by arranging the momentum vectors in
456-
// a triangle and finding the required angles between the vectors.
457-
const amrex::ParticleReal cos_alpha = (n_0*n_0 + n_1*n_1 - n_2*n_2) / (2.0_prt * n_0 * n_1);
458-
const amrex::ParticleReal sin_alpha = std::sqrt(1.0_prt - cos_alpha*cos_alpha);
459-
const amrex::ParticleReal cos_gamma = (n_0*n_0 + n_2*n_2 - n_1*n_1) / (2.0_prt * n_0 * n_2);
460-
const amrex::ParticleReal sin_gamma = std::sqrt(1.0_prt - cos_gamma*cos_gamma);
461-
462-
// choose random theta and phi values (orientation of the triangle)
463-
const amrex::ParticleReal Theta = amrex::Random(engine) * 2.0_prt * MathConst::pi;
464-
const amrex::ParticleReal phi = amrex::Random(engine) * MathConst::pi;
465-
466-
const amrex::ParticleReal cos_theta = std::cos(Theta);
467-
const amrex::ParticleReal sin_theta = std::sin(Theta);
468-
const amrex::ParticleReal cos_phi = std::cos(phi);
469-
const amrex::ParticleReal sin_phi = std::sin(phi);
470-
471-
// calculate the velocity components for each particle
472-
ux1 = C * n_0 / m1 * cos_theta * cos_phi;
473-
uy1 = C * n_0 / m1 * cos_theta * sin_phi;
474-
uz1 = -C * n_0 / m1 * sin_theta;
475-
476-
ux_p1 = C * n_1 / m_ioniz_product1 * (-cos_alpha * cos_theta * cos_phi - sin_alpha * sin_phi);
477-
uy_p1 = C * n_1 / m_ioniz_product1 * (-cos_alpha * cos_theta * sin_phi + sin_alpha * cos_phi);
478-
uz_p1 = C * n_1 / m_ioniz_product1 * (cos_alpha * sin_theta);
479-
480-
ux_p2 = C * n_2 / m_ioniz_product2 * (-cos_gamma * cos_theta * cos_phi + sin_gamma * sin_phi);
481-
uy_p2 = C * n_2 / m_ioniz_product2 * (-cos_gamma * cos_theta * sin_phi - sin_gamma * cos_phi);
482-
uz_p2 = C * n_2 / m_ioniz_product2 * (cos_gamma * sin_theta);
395+
// calculate kinetic energy of the collision (in eV)
396+
const amrex::ParticleReal E1 = (
397+
0.5_prt * m1 * (ux1*ux1 + uy1*uy1 + uz1*uz1) / PhysConst::q_e
398+
);
399+
const amrex::ParticleReal E2 = (
400+
0.5_prt * m2 * (ux_p2*ux_p2 + uy_p2*uy_p2 + uz_p2*uz_p2) / PhysConst::q_e
401+
);
402+
const amrex::ParticleReal E_coll = E1 + E2;
403+
404+
// subtract the energy cost for ionization
405+
const amrex::ParticleReal E_out = (E_coll - ionization_energy) * PhysConst::q_e;
406+
407+
// Momentum and energy division after the ionization event
408+
// is done as follows:
409+
// Three numbers are generated that satisfy the triangle
410+
// inequality. These numbers will be scaled to give the
411+
// momentum for each particle (satisfying the triangle
412+
// inequality ensures that the momentum vectors can be
413+
// arranged such that the net linear momentum is zero).
414+
// Product 2 is definitely an ion, so we first choose its
415+
// proportion to be n_2 = min(E2 / E_coll, 0.5 * R), where
416+
// R is a random number between 0 and 1.
417+
// The other two numbers (n_0 & n_1) are then found
418+
// by choosing a random point on the ellipse characterized
419+
// by the semi-major and semi-minor axes
420+
// a = (1 - n_0) / 2.0
421+
// b = 0.5 * sqrt(1 - 2 * n_0)
422+
// The numbers are found by randomly sampling an x value
423+
// between -a and a, and finding the corresponding y value
424+
// that falls on the ellipse: y^2 = b^2 - b^2/a^2 * x^2.
425+
// Then n_0 = sqrt(y^2 + (x - n_2/2)^2) and
426+
// n_1 = 1 - n_2 - n_0.
427+
// Next, we need to find the number C, such that if
428+
// p_i = C * n_i, we would have E_0 + E_1 + E_2 = E_out
429+
// where E_i = p_i^2 / (2 M_i), i.e.,
430+
// C^2 * \sum_i [n_i^2 / (2 M_i)] = E_out
431+
// After C is determined the momentum vectors are arranged
432+
// such that the net linear momentum is zero.
433+
434+
const amrex::ParticleReal n_2 = std::min(E2 / E_coll, 0.5_prt * amrex::Random(engine));
435+
436+
// find ellipse semi-major and minor axis
437+
const amrex::ParticleReal a = 0.5_prt * (1.0_prt - n_2);
438+
const amrex::ParticleReal b = 0.5_prt * std::sqrt(1.0_prt - 2.0_prt * n_2);
439+
440+
// sample random x value and calculate y
441+
const amrex::ParticleReal x = (2._prt * amrex::Random(engine) - 1.0_prt) * a;
442+
const amrex::ParticleReal y2 = b*b - b*b/(a*a) * x*x;
443+
const amrex::ParticleReal n_0 = std::sqrt(y2 + x*x - x*n_2 + 0.25_prt*n_2*n_2);
444+
const amrex::ParticleReal n_1 = 1.0_prt - n_0 - n_2;
445+
446+
// calculate the value of C
447+
const amrex::ParticleReal C = std::sqrt(E_out / (
448+
n_0*n_0 / (2.0_prt * m1) + n_1*n_1 / (2.0_prt * m_ioniz_product1) + n_2*n_2 / (2.0_prt * m_ioniz_product2)
449+
));
450+
451+
// Now that appropriate momenta are set for each outgoing species
452+
// the directions for the velocity vectors must be chosen such
453+
// that the net linear momentum in the current frame is 0.
454+
// This is achieved by arranging the momentum vectors in
455+
// a triangle and finding the required angles between the vectors.
456+
const amrex::ParticleReal cos_alpha = (n_0*n_0 + n_1*n_1 - n_2*n_2) / (2.0_prt * n_0 * n_1);
457+
const amrex::ParticleReal sin_alpha = std::sqrt(1.0_prt - cos_alpha*cos_alpha);
458+
const amrex::ParticleReal cos_gamma = (n_0*n_0 + n_2*n_2 - n_1*n_1) / (2.0_prt * n_0 * n_2);
459+
const amrex::ParticleReal sin_gamma = std::sqrt(1.0_prt - cos_gamma*cos_gamma);
460+
461+
// choose random theta and phi values (orientation of the triangle)
462+
const amrex::ParticleReal Theta = amrex::Random(engine) * 2.0_prt * MathConst::pi;
463+
const amrex::ParticleReal phi = amrex::Random(engine) * MathConst::pi;
464+
465+
const amrex::ParticleReal cos_theta = std::cos(Theta);
466+
const amrex::ParticleReal sin_theta = std::sin(Theta);
467+
const amrex::ParticleReal cos_phi = std::cos(phi);
468+
const amrex::ParticleReal sin_phi = std::sin(phi);
469+
470+
// calculate the velocity components for each particle
471+
ux1 = C * n_0 / m1 * cos_theta * cos_phi;
472+
uy1 = C * n_0 / m1 * cos_theta * sin_phi;
473+
uz1 = -C * n_0 / m1 * sin_theta;
474+
475+
ux_p1 = C * n_1 / m_ioniz_product1 * (-cos_alpha * cos_theta * cos_phi - sin_alpha * sin_phi);
476+
uy_p1 = C * n_1 / m_ioniz_product1 * (-cos_alpha * cos_theta * sin_phi + sin_alpha * cos_phi);
477+
uz_p1 = C * n_1 / m_ioniz_product1 * (cos_alpha * sin_theta);
478+
479+
ux_p2 = C * n_2 / m_ioniz_product2 * (-cos_gamma * cos_theta * cos_phi + sin_gamma * sin_phi);
480+
uy_p2 = C * n_2 / m_ioniz_product2 * (-cos_gamma * cos_theta * sin_phi - sin_gamma * cos_phi);
481+
uz_p2 = C * n_2 / m_ioniz_product2 * (cos_gamma * sin_theta);
483482

484-
}
485-
else {
486-
amrex::Abort("Unknown scattering process.");
487-
}
488483
// transform back to labframe
489484
ux1 += uCOM_x;
490485
uy1 += uCOM_y;

0 commit comments

Comments
 (0)