Skip to content

Remove clight factor from u #1311

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
MaxThevenet merged 10 commits into from
Nov 25, 2025
Merged

Remove clight factor from u #1311

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
bca7c31
Remove clight factor from u
AlexanderSinn Nov 12, 2025
bd7d7ae
fix openPMD output
AlexanderSinn Nov 12, 2025
6131df6
fix collisions
AlexanderSinn Nov 12, 2025
77793e5
fix adaptive timestep
AlexanderSinn Nov 12, 2025
c38a00e
try wrong factor with RR
AlexanderSinn Nov 12, 2025
d164ef5
try 0.9
AlexanderSinn Nov 12, 2025
6f304b0
reset rr checksum
AlexanderSinn Nov 12, 2025
d00a71f
Merge branch 'development' into Remove_clight_factor_from_u
AlexanderSinn Nov 21, 2025
4621712
fix radiation reactions in normalized units
AlexanderSinn Nov 21, 2025
df8b72f
remove proper velocity from comments
AlexanderSinn Nov 25, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/diagnostics/OpenPMDWriter.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -437,8 +437,8 @@ OpenPMDWriter::SetupPos (openPMD::ParticleSpecies& currSpecies, BeamParticleCont
// calculate the multiplier to convert from Hipace to SI units
double hipace_to_SI_pos = 1.;
double hipace_to_SI_weight = 1.;
double hipace_to_SI_momentum = beam.m_mass;
double hipace_to_unitSI_momentum = beam.m_mass;
double hipace_to_SI_momentum = beam.m_mass * phys_const_SI.c;
double hipace_to_unitSI_momentum = beam.m_mass * phys_const_SI.c;
double hipace_to_SI_charge = 1.;
double hipace_to_SI_mass = 1.;

Expand Down
7 changes: 3 additions & 4 deletions src/particles/beam/BeamParticleContainer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -505,7 +505,6 @@ BeamParticleContainer::InSituComputeDiags (int islice)

const amrex::Real insitu_radius_sq = m_insitu_radius * m_insitu_radius;
const PhysConst phys_const = get_phys_const();
const amrex::Real clight_inv = 1.0_rt/phys_const.c;
const auto ptd = getBeamSlice(WhichBeamSlice::This).getParticleTileData();

amrex::TypeMultiplier<amrex::ReduceOps, amrex::ReduceOpSum[m_insitu_nrp + m_insitu_nip]> reduce_op;
Expand All @@ -518,9 +517,9 @@ BeamParticleContainer::InSituComputeDiags (int islice)
const amrex::Real x = ptd.pos(0, ip);
const amrex::Real y = ptd.pos(1, ip);
const amrex::Real z = ptd.pos(2, ip);
const amrex::Real ux = ptd.rdata(BeamIdx::ux)[ip] * clight_inv; // proper velocity to u
const amrex::Real uy = ptd.rdata(BeamIdx::uy)[ip] * clight_inv;
const amrex::Real uz = ptd.rdata(BeamIdx::uz)[ip] * clight_inv;
const amrex::Real ux = ptd.rdata(BeamIdx::ux)[ip];
const amrex::Real uy = ptd.rdata(BeamIdx::uy)[ip];
const amrex::Real uz = ptd.rdata(BeamIdx::uz)[ip];
const amrex::Real w = ptd.rdata(BeamIdx::w)[ip];

const amrex::Real uz_inv = uz == 0._rt ? 0._rt : 1._rt / uz;
Expand Down
72 changes: 32 additions & 40 deletions src/particles/beam/BeamParticleContainerInit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,6 @@ namespace
* \param[in] weight weight of the single particle
* \param[in] pid particle ID to be assigned to the particle
* \param[in] ip index of the particle
* \param[in] speed_of_light speed of light in the current units
* \param[in] enforceBC functor to enforce the boundary condition
*/
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
Expand All @@ -46,20 +45,20 @@ namespace
const amrex::Real& ux, const amrex::Real& uy, const amrex::Real& uz,
const amrex::Real& sx, const amrex::Real& sy, const amrex::Real& sz,
const amrex::Real& weight, const amrex::Long pid, const amrex::Long ip,
const amrex::Real& speed_of_light, const EnforceBC& enforceBC, bool do_spin) noexcept
const EnforceBC& enforceBC, bool do_spin) noexcept
{
amrex::Real xp = x;
amrex::Real yp = y;
amrex::Real uxp = ux * speed_of_light;
amrex::Real uyp = uy * speed_of_light;
amrex::Real uxp = ux;
amrex::Real uyp = uy;
if (enforceBC(ptd, ip, xp, yp, uxp, uyp, BeamIdx::w)) return;

ptd.rdata(BeamIdx::x )[ip] = xp;
ptd.rdata(BeamIdx::y )[ip] = yp;
ptd.rdata(BeamIdx::z )[ip] = z;
ptd.rdata(BeamIdx::ux )[ip] = uxp;
ptd.rdata(BeamIdx::uy )[ip] = uyp;
ptd.rdata(BeamIdx::uz )[ip] = uz * speed_of_light;
ptd.rdata(BeamIdx::x )[ip] = xp;
ptd.rdata(BeamIdx::y )[ip] = yp;
ptd.rdata(BeamIdx::z )[ip] = z;
ptd.rdata(BeamIdx::ux)[ip] = uxp;
ptd.rdata(BeamIdx::uy)[ip] = uyp;
ptd.rdata(BeamIdx::uz)[ip] = uz;
if (do_spin) {
ptd.m_runtime_rdata[0][ip] = sx;
ptd.m_runtime_rdata[1][ip] = sy;
Expand All @@ -83,7 +82,6 @@ namespace
* \param[in] weight weight of the single particle
* \param[in] pid particle ID to be assigned to the particle at index 0
* \param[in] ip index of the particle
* \param[in] speed_of_light speed of light in the current units
* \param[in] enforceBC functor to enforce the boundary condition
* \param[in] is_valid if the particle is valid
*/
Expand All @@ -92,22 +90,22 @@ namespace
const BeamTile::ParticleTileDataType& ptd, const amrex::Real x,
const amrex::Real y, const amrex::Real z, const amrex::Real ux, const amrex::Real uy,
const amrex::Real uz, const amrex::Real weight, const amrex::Long pid,
const amrex::Long ip, const amrex::Real speed_of_light, const EnforceBC& enforceBC,
const amrex::Long ip, const EnforceBC& enforceBC,
const bool is_valid=true) noexcept
{
amrex::Real xp = x;
amrex::Real yp = y;
amrex::Real uxp = ux * speed_of_light;
amrex::Real uyp = uy * speed_of_light;
amrex::Real uxp = ux;
amrex::Real uyp = uy;
if (enforceBC(ptd, ip, xp, yp, uxp, uyp, BeamIdx::w)) return;

ptd.rdata(BeamIdx::x )[ip] = xp;
ptd.rdata(BeamIdx::y )[ip] = yp;
ptd.rdata(BeamIdx::z )[ip] = z;
ptd.rdata(BeamIdx::ux )[ip] = uxp;
ptd.rdata(BeamIdx::uy )[ip] = uyp;
ptd.rdata(BeamIdx::uz )[ip] = uz * speed_of_light;
ptd.rdata(BeamIdx::w )[ip] = is_valid ? std::abs(weight) : amrex::Real{0};
ptd.rdata(BeamIdx::x )[ip] = xp;
ptd.rdata(BeamIdx::y )[ip] = yp;
ptd.rdata(BeamIdx::z )[ip] = z;
ptd.rdata(BeamIdx::ux)[ip] = uxp;
ptd.rdata(BeamIdx::uy)[ip] = uyp;
ptd.rdata(BeamIdx::uz)[ip] = uz;
ptd.rdata(BeamIdx::w )[ip] = is_valid ? std::abs(weight) : amrex::Real{0};

ptd.idata(BeamIdx::nsubcycles)[ip] = 0;
ptd.idata(BeamIdx::mr_level)[ip] = 0;
Expand Down Expand Up @@ -298,8 +296,6 @@ InitBeamFixedPPCSlice (const int islice, const int which_beam_slice)
const uint64_t pid = m_id64;
m_id64 += num_to_add;

const amrex::Real speed_of_light = get_phys_const().c;

const auto enforceBC = EnforceBC();

amrex::ParallelForRNG(to2D(slice_box),
Expand Down Expand Up @@ -345,7 +341,7 @@ InitBeamFixedPPCSlice (const int islice, const int which_beam_slice)
const amrex::Real weight = density * scale_fac;

AddOneBeamParticleSlice(ptd, x, y, z, u[0], u[1], u[2], weight,
pid, pidx, speed_of_light, enforceBC, true);
pid, pidx, enforceBC, true);

++pidx;
}
Expand Down Expand Up @@ -404,8 +400,6 @@ InitBeamFixedWeightSlice (int slice, int which_slice)

if (num_to_add == 0) return;

const amrex::Real clight = get_phys_const().c;

auto& particle_tile = getBeamSlice(which_slice);

// Access particles' SoA
Expand Down Expand Up @@ -461,21 +455,21 @@ InitBeamFixedWeightSlice (int slice, int which_slice)
{
AddOneBeamParticleSlice(ptd, cental_x_pos+x, cental_y_pos+y,
z_central, u[0], u[1], u[2], weight,
pid, i, clight, enforceBC, is_valid);
pid, i, enforceBC, is_valid);

} else {
AddOneBeamParticleSlice(ptd, cental_x_pos+x, cental_y_pos+y,
z_central, u[0], u[1], u[2], weight,
pid, 4*i, clight, enforceBC, is_valid);
pid, 4*i, enforceBC, is_valid);
AddOneBeamParticleSlice(ptd, cental_x_pos-x, cental_y_pos+y,
z_central, -u[0], u[1], u[2], weight,
pid, 4*i+1, clight, enforceBC, is_valid);
pid, 4*i+1, enforceBC, is_valid);
AddOneBeamParticleSlice(ptd, cental_x_pos+x, cental_y_pos-y,
z_central, u[0], -u[1], u[2], weight,
pid, 4*i+2, clight, enforceBC, is_valid);
pid, 4*i+2, enforceBC, is_valid);
AddOneBeamParticleSlice(ptd, cental_x_pos-x, cental_y_pos-y,
z_central, -u[0], -u[1], u[2], weight,
pid, 4*i+3, clight, enforceBC, is_valid);
pid, 4*i+3, enforceBC, is_valid);
}
});

Expand Down Expand Up @@ -618,7 +612,6 @@ InitBeamFixedWeightPDFSlice (int slice, int which_slice)
// Access particles' SoA
const auto ptd = particle_tile.getParticleTileData();

const amrex::Real clight = get_phys_const().c;
const bool do_symmetrize = m_do_symmetrize;
const bool peak_density_is_specified = m_peak_density_is_specified;
const amrex::Real z_foc = m_z_foc;
Expand Down Expand Up @@ -685,17 +678,17 @@ InitBeamFixedWeightPDFSlice (int slice, int which_slice)
if (!do_symmetrize)
{
AddOneBeamParticleSlice(ptd, x_mean+x, y_mean+y,
z, ux, uy, uz, weight, pid, i, clight, enforceBC, is_valid);
z, ux, uy, uz, weight, pid, i, enforceBC, is_valid);

} else {
AddOneBeamParticleSlice(ptd, x_mean+x, y_mean+y,
z, ux, uy, uz, weight, pid, 4*i, clight, enforceBC, is_valid);
z, ux, uy, uz, weight, pid, 4*i, enforceBC, is_valid);
AddOneBeamParticleSlice(ptd, x_mean-x, y_mean+y,
z, -ux, uy, uz, weight, pid, 4*i+1, clight, enforceBC, is_valid);
z, -ux, uy, uz, weight, pid, 4*i+1, enforceBC, is_valid);
AddOneBeamParticleSlice(ptd, x_mean+x, y_mean-y,
z, ux, -uy, uz, weight, pid, 4*i+2, clight, enforceBC, is_valid);
z, ux, -uy, uz, weight, pid, 4*i+2, enforceBC, is_valid);
AddOneBeamParticleSlice(ptd, x_mean-x, y_mean-y,
z, -ux, -uy, uz, weight, pid, 4*i+3, clight, enforceBC, is_valid);
z, -ux, -uy, uz, weight, pid, 4*i+3, enforceBC, is_valid);
}
});

Expand Down Expand Up @@ -756,7 +749,6 @@ InitBeamFromList3D ()

const auto ptd = particle_tile.getParticleTileData();
const auto enforceBC = EnforceBC();
const amrex::Real clight = get_phys_const().c;

const uint64_t pid = m_id64;
m_id64 += m_num_particles_list;
Expand All @@ -770,7 +762,7 @@ InitBeamFromList3D ()
do_spin_tracking ? p_sy[i] : 0.,
do_spin_tracking ? p_sz[i] : 0.,
p_w[i],
pid, i, clight, enforceBC, do_spin_tracking);
pid, i, enforceBC, do_spin_tracking);
});

amrex::Gpu::streamSynchronize();
Expand Down Expand Up @@ -1213,7 +1205,7 @@ InitBeamFromFile (const std::string input_file,
do_spin_tracking ? static_cast<amrex::Real>(s_y_ptr[i]) : 0.,
do_spin_tracking ? static_cast<amrex::Real>(s_z_ptr[i]) : 0.,
static_cast<amrex::Real>(w_w_ptr[i] * unit_ww),
pid, i, phys_const.c, enforceBC, do_spin_tracking);
pid, i, enforceBC, do_spin_tracking);
});

amrex::Gpu::streamSynchronize();
Expand Down
19 changes: 10 additions & 9 deletions src/particles/collisions/ComputeTemperature.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ AMREX_GPU_HOST_DEVICE
T_R ComputeTemperature (
T_index const Is, T_index const Ie, T_index const * AMREX_RESTRICT I,
T_R const * AMREX_RESTRICT ux, T_R const * AMREX_RESTRICT uy, T_R const * AMREX_RESTRICT psi,
T_R const m, T_R const clight, T_R const inv_c2, bool is_beam_coll )
T_R const m, T_R const clight, bool is_beam_coll )
{
using namespace amrex::literals;

Expand All @@ -27,18 +27,19 @@ T_R ComputeTemperature (
for (int i = Is; i < (int) Ie; ++i)
{
// particle's Lorentz factor
gm = is_beam_coll ? std::sqrt( 1._rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]]
+ psi[I[i]]*psi[I[i]])*inv_c2 )
: (1.0_rt + (ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]])*inv_c2
gm = is_beam_coll ? std::sqrt( 1._rt + ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]]
+ psi[I[i]]*psi[I[i]] )
: (1.0_rt + ux[I[i]]*ux[I[i]] + uy[I[i]]*uy[I[i]]
+ psi[I[i]]*psi[I[i]]) / (2.0_rt * psi[I[i]] );
const amrex::Real uz = is_beam_coll ? psi[I[i]] : clight * (gm - psi[I[i]]);
const amrex::Real uz = is_beam_coll ? psi[I[i]] : gm - psi[I[i]];
const amrex::Real gmi = clight / gm;
us = ( ux[ I[i] ] * ux[ I[i] ] +
uy[ I[i] ] * uy[ I[i] ] +
uz * uz);
vx += ux[ I[i] ] / gm;
vy += uy[ I[i] ] / gm;
vz += uz / gm;
vs += us / gm / gm;
vx += ux[ I[i] ] * gmi;
vy += uy[ I[i] ] * gmi;
vz += uz * gmi;
vs += us * gmi * gmi;
}

vx = vx / N; vy = vy / N;
Expand Down
10 changes: 3 additions & 7 deletions src/particles/collisions/CoulombCollision.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,6 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
bool normalized_units = Hipace::m_normalized_units;

const amrex::Real clight = cst.c;
const amrex::Real inv_c = 1.0_rt / cst.c;
const amrex::Real inv_c2 = 1.0_rt / ( cst.c * cst.c );
constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c;
constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c );

Expand Down Expand Up @@ -135,7 +133,7 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
indices1, indices1,
ux1, uy1, psi1, ux1, uy1, psi1, w1, w1, ion_lev1, ion_lev1,
q1, q1, m1, m1, -1.0_rt, -1.0_rt, can_ionize1, can_ionize1,
dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
dt, CoulombLog, inv_dV, clight, inv_c_SI, inv_c2_SI,
normalized_units, background_density_SI, is_same_species, false, engine );
}
);
Expand Down Expand Up @@ -226,7 +224,7 @@ CoulombCollision::doPlasmaPlasmaCoulombCollision (
indices1, indices2,
ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev1, ion_lev2,
q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2,
dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
dt, CoulombLog, inv_dV, clight, inv_c_SI, inv_c2_SI,
normalized_units, background_density_SI, is_same_species, false, engine );
}
);
Expand All @@ -253,8 +251,6 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
bool normalized_units = Hipace::m_normalized_units;

const amrex::Real clight = cst.c;
const amrex::Real inv_c = 1.0_rt / cst.c;
const amrex::Real inv_c2 = 1.0_rt / ( cst.c * cst.c );
constexpr amrex::Real inv_c_SI = 1.0_rt / PhysConstSI::c;
constexpr amrex::Real inv_c2_SI = 1.0_rt / ( PhysConstSI::c * PhysConstSI::c );

Expand Down Expand Up @@ -339,7 +335,7 @@ CoulombCollision::doBeamPlasmaCoulombCollision (
indices1, indices2,
ux1, uy1, psi1, ux2, uy2, psi2, w1, w2, ion_lev2, ion_lev2, // passing ion_lev2 for beam particles, will never be used
q1, q2, m1, m2, -1.0_rt, -1.0_rt, can_ionize1, can_ionize2,
dt, CoulombLog, inv_dV, clight, inv_c, inv_c_SI, inv_c2, inv_c2_SI,
dt, CoulombLog, inv_dV, clight, inv_c_SI, inv_c2_SI,
normalized_units, background_density_SI, false, true, engine );
}
);
Expand Down
Loading
Loading