Skip to content
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

Per-slice fixed_weight beam initialization #1008

Merged
merged 5 commits into from
Aug 29, 2023
Merged

Per-slice fixed_weight beam initialization #1008

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
88afb2f
Per slice fixed_weight beam initialization
AlexanderSinn Aug 9, 2023
5974ca8
some fixes
AlexanderSinn Aug 10, 2023
c82e20a
add missing resize to fix next_deposition CI test
AlexanderSinn Aug 27, 2023
25d5be2
reset CI because RNG changed
AlexanderSinn Aug 27, 2023
7a1aafb
fix weight
AlexanderSinn Aug 27, 2023
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
2 changes: 1 addition & 1 deletion examples/gaussian_weight/analysis.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@
assert(np.abs((np.average(xp)-x_avg)) < 5e-7)
assert(np.abs((np.average(yp)-y_avg)/y_avg) < .03)

assert( np.abs((np.average(zp)-z_avg)/z_avg) < .035)
assert( np.abs((np.average(zp)-z_avg)/z_avg) < .05)
assert(np.abs((np.std(xp)-x_std)/x_std) < .03)
assert(np.abs((np.std(yp)-y_std)/y_std) < .03)
assert(np.abs((np.std(zp)-z_std)/z_std) < .03)
36 changes: 19 additions & 17 deletions src/particles/beam/BeamParticleContainer.H
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,17 +87,8 @@ public:
void InitBeamFixedPPCSlice (const int islice, const int which_beam_slice);

/** Initialize a beam with a fix number of particles, and fixed weight */
void InitBeamFixedWeight (int num_to_add,
const GetInitialMomentum& get_momentum,
const amrex::ParserExecutor<1>& pos_mean_x,
const amrex::ParserExecutor<1>& pos_mean_y,
const amrex::Real pos_mean_z,
const amrex::RealVect pos_std,
const amrex::Real total_charge,
const amrex::Real z_foc,
const bool do_symmetrize,
const bool can, const amrex::Real zmin, const amrex::Real zmax,
const amrex::Real radius);
void InitBeamFixedWeight3D ();
void InitBeamFixedWeightSlice (const int islice, const int which_beam_slice);

#ifdef HIPACE_USE_OPENPMD
/** Checks the input file first to determine its Datatype
Expand Down Expand Up @@ -213,37 +204,48 @@ private:
std::string m_name; /**< name of the species */
/** injection type, fixed_width or fixed_ppc */
std::string m_injection_type;
uint64_t m_id64 = 1; /**< 64 bit ID to initialize many particles without overflowing */
/** Min longitudinal particle position of the beam */
amrex::Real m_zmin = -std::numeric_limits<amrex::Real>::infinity();
/** Max longitudinal particle position of the beam */
amrex::Real m_zmax = std::numeric_limits<amrex::Real>::infinity();
amrex::Real m_radius {std::numeric_limits<amrex::Real>::infinity()}; /**< Radius of the beam */
GetInitialMomentum m_get_momentum {}; /**< momentum profile of the beam */

// fixed_ppc:

amrex::IntVect m_ppc {1, 1, 1}; /**< Number of particles per cell in each direction */
amrex::Real m_zmin; /**< Min longitudinal particle position of the beam */
amrex::Real m_zmax; /**< Max longitudinal particle position of the beam */
amrex::Real m_radius {std::numeric_limits<amrex::Real>::infinity()}; /**< Radius of the beam */
amrex::RealVect m_position_mean {0., 0., 0.}; /**< mean position of the beam */
amrex::Real m_min_density {0.}; /**< minimum density at which beam particles are generated */
amrex::IntVect m_random_ppc {0, 0, 0}; /**< if the cell position is random in each direction */
GetInitialDensity m_get_density {}; /**< density profile of the beam */
GetInitialMomentum m_get_momentum {}; /**< momentum profile of the beam */
/** Density parser for fixed-ppc beam. Owns data for m_density_func */
amrex::Parser m_density_parser;
uint64_t m_id64 = 1; /**< 64 bit ID to initialize many particles without overflowing */

// fixed_weight:

bool m_can_profile = false;
/** Average x position of the fixed-weight beam depending on z */
amrex::Parser m_pos_mean_x_parser;
/** Average x position of the fixed-weight beam depending on z */
amrex::ParserExecutor<1> m_pos_mean_x_func;
/** Average y position of the fixed-weight beam depending on z */
amrex::Parser m_pos_mean_y_parser;
/* Average y position of the fixed-weight beam depending on z */
amrex::ParserExecutor<1> m_pos_mean_y_func;
/* Average z position of the fixed-weight beam depending on z */
amrex::Real m_pos_mean_z = 0;
/** Width of the Gaussian beam. Only used for a fixed-weight beam */
amrex::RealVect m_position_std {0., 0., 0.};
/** Distance at which the beam is focused, starting from its initial position */
amrex::Real m_z_foc {0.};
amrex::Real m_duz_per_uz0_dzeta {0.}; /**< relative energy spread per dzeta */
int m_num_particles; /**< Number of particles for fixed-weigth Gaussian beam */
amrex::Long m_num_particles; /**< Number of particles for fixed-weight Gaussian beam */
amrex::Real m_total_charge; /**< Total beam charge for fixed-weight Gaussian beam */
amrex::Real m_density; /**< Peak density for fixed-weight Gaussian beam */
bool m_do_symmetrize {0}; /**< Option to symmetrize the beam */
/** Array for the z position of all beam particles */
amrex::PODVector<amrex::Real, amrex::PolymorphicArenaAllocator<amrex::Real>> m_z_array {};
Comment on lines +247 to +248
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In principle, could this be always float if memory is a problem?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In that we would also loose precision. I think we would be better off using inverr to make the z distribution per slice.


// from_file:

Expand Down
46 changes: 24 additions & 22 deletions src/particles/beam/BeamParticleContainer.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,22 +104,19 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom)

} else if (m_injection_type == "fixed_weight") {

bool can = false;
amrex::Real zmin = -std::numeric_limits<amrex::Real>::infinity();
amrex::Real zmax = std::numeric_limits<amrex::Real>::infinity();
std::string profile = "gaussian";
queryWithParser(pp, "profile", profile);
queryWithParser(pp, "radius", m_radius);
if (profile == "can") {
can = true;
getWithParser(pp, "zmin", zmin);
getWithParser(pp, "zmax", zmax);
m_can_profile = true;
getWithParser(pp, "zmin", m_zmin);
getWithParser(pp, "zmax", m_zmax);
} else if (profile == "gaussian") {
queryWithParser(pp, "zmin", zmin);
queryWithParser(pp, "zmax", zmax);
queryWithParser(pp, "zmin", m_zmin);
queryWithParser(pp, "zmax", m_zmax);
AMREX_ALWAYS_ASSERT_WITH_MESSAGE( !m_do_salame ||
(zmin != -std::numeric_limits<amrex::Real>::infinity() &&
zmax != std::numeric_limits<amrex::Real>::infinity()),
(m_zmin != -std::numeric_limits<amrex::Real>::infinity() &&
m_zmax != std::numeric_limits<amrex::Real>::infinity()),
"For the SALAME algorithm it is mandatory to either use a 'can' profile or "
"'zmin' and 'zmax' with a gaussian profile");
} else {
Expand All @@ -129,10 +126,9 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom)
std::array<std::string, 3> pos_mean_arr{"","",""};
getWithParser(pp, "position_mean", pos_mean_arr);
queryWithParser(pp, "z_foc", m_z_foc);
auto pos_mean_x = makeFunctionWithParser<1>(pos_mean_arr[0], m_pos_mean_x_parser, {"z"});
auto pos_mean_y = makeFunctionWithParser<1>(pos_mean_arr[1], m_pos_mean_y_parser, {"z"});
amrex::Real pos_mean_z = 0;
Parser::fillWithParser(pos_mean_arr[2], pos_mean_z);
m_pos_mean_x_func = makeFunctionWithParser<1>(pos_mean_arr[0], m_pos_mean_x_parser, {"z"});
m_pos_mean_y_func = makeFunctionWithParser<1>(pos_mean_arr[1], m_pos_mean_y_parser, {"z"});
Parser::fillWithParser(pos_mean_arr[2], m_pos_mean_z);

getWithParser(pp, "position_std", m_position_std);
getWithParser(pp, "num_particles", m_num_particles);
Expand Down Expand Up @@ -161,10 +157,13 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom)
m_total_charge *= geom.InvCellSize(0)*geom.InvCellSize(1)*geom.InvCellSize(2);
}

const GetInitialMomentum get_momentum(m_name);
InitBeamFixedWeight(m_num_particles, get_momentum, pos_mean_x, pos_mean_y, pos_mean_z,
m_position_std, m_total_charge, m_z_foc, m_do_symmetrize, can, zmin, zmax, m_radius);
m_total_num_particles = getBeamInitSlice().size();
m_get_momentum = GetInitialMomentum{m_name};
InitBeamFixedWeight3D();
m_total_num_particles = m_num_particles;
if (Hipace::HeadRank()) {
m_init_sorter.sortParticlesByBox(m_z_array.dataPtr(), m_z_array.size(),
m_initialize_on_cpu, geom);
}

} else if (m_injection_type == "from_file") {
#ifdef HIPACE_USE_OPENPMD
Expand All @@ -182,6 +181,11 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom)
geom, m_plasma_density, m_num_iteration, m_species_name,
species_specified);
m_total_num_particles = getBeamInitSlice().size();
if (Hipace::HeadRank()) {
m_init_sorter.sortParticlesByBox(
getBeamInitSlice().GetStructOfArrays().GetRealData(BeamIdx::z).dataPtr(),
getBeamInitSlice().size(), m_initialize_on_cpu, geom);
}
#else
amrex::Abort("beam particle injection via external_file requires openPMD support: "
"Add HiPACE_OPENPMD=ON when compiling HiPACE++.\n");
Expand All @@ -200,10 +204,6 @@ BeamParticleContainer::InitData (const amrex::Geometry& geom)
amrex::ParallelDescriptor::Communicator());
#endif

if (Hipace::HeadRank() && m_injection_type != "fixed_ppc") {
m_init_sorter.sortParticlesByBox(*this, geom);
}

if (m_insitu_period > 0) {
#ifdef HIPACE_USE_OPENPMD
AMREX_ALWAYS_ASSERT_WITH_MESSAGE(m_insitu_file_prefix !=
Expand Down Expand Up @@ -275,6 +275,8 @@ BeamParticleContainer::intializeSlice (int slice, int which_slice) {

if (m_injection_type == "fixed_ppc") {
InitBeamFixedPPCSlice(slice, which_slice);
} else if (m_injection_type == "fixed_weight") {
InitBeamFixedWeightSlice(slice, which_slice);
} else {
const int num_particles = m_init_sorter.m_box_counts_cpu[slice];

Expand Down
Loading