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.

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

Store a single process-integrated energy loss and range table per particle #1536

Merged
merged 3 commits into from
Dec 6, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
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
5 changes: 2 additions & 3 deletions src/celeritas/alongstep/detail/FluctELoss.hh
Original file line number Diff line number Diff line change
Expand Up @@ -81,9 +81,8 @@ CELER_FUNCTION bool FluctELoss::is_applicable(CoreTrackView const& track) const
if (track.make_sim_view().status() == TrackStatus::errored)
return false;

// Energy loss grid ID will be 'false' if inapplicable
auto ppid = track.make_physics_view().eloss_ppid();
return static_cast<bool>(ppid);
// Energy loss grid ID is 'false'
return static_cast<bool>(track.make_physics_view().energy_loss_grid());
}

//---------------------------------------------------------------------------//
Expand Down
5 changes: 2 additions & 3 deletions src/celeritas/alongstep/detail/MeanELoss.hh
Original file line number Diff line number Diff line change
Expand Up @@ -53,9 +53,8 @@ CELER_FUNCTION bool MeanELoss::is_applicable(CoreTrackView const& track) const
if (track.make_sim_view().status() == TrackStatus::errored)
return false;

// Energy loss grid ID will be 'false' if inapplicable
auto ppid = track.make_physics_view().eloss_ppid();
return static_cast<bool>(ppid);
// Energy loss grid ID is 'false'
return static_cast<bool>(track.make_physics_view().energy_loss_grid());
}

//---------------------------------------------------------------------------//
Expand Down
11 changes: 3 additions & 8 deletions src/celeritas/em/msc/detail/UrbanMscHelper.hh
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,6 @@
#include "celeritas/grid/EnergyLossCalculator.hh"
#include "celeritas/grid/InverseRangeCalculator.hh"
#include "celeritas/grid/RangeCalculator.hh"
#include "celeritas/grid/ValueGridType.hh"
#include "celeritas/grid/XsCalculator.hh"
#include "celeritas/phys/ParticleTrackView.hh"
#include "celeritas/phys/PhysicsTrackView.hh"
Expand Down Expand Up @@ -145,10 +144,8 @@ CELER_FUNCTION real_type UrbanMscHelper::calc_msc_mfp(Energy energy) const
CELER_FUNCTION auto
UrbanMscHelper::calc_inverse_range(real_type step) const -> Energy
{
auto range_gid
= physics_.value_grid(ValueGridType::range, physics_.eloss_ppid());
auto range_to_energy
= physics_.make_calculator<InverseRangeCalculator>(range_gid);
auto range_to_energy = physics_.make_calculator<InverseRangeCalculator>(
physics_.range_grid());
return range_to_energy(step);
}

Expand All @@ -172,11 +169,9 @@ UrbanMscHelper::calc_end_energy(real_type step) const -> Energy
real_type range = physics_.dedx_range();
if (step <= range * shared_.params.dtrl())
{
auto eloss_gid = physics_.value_grid(ValueGridType::energy_loss,
physics_.eloss_ppid());
// Assume constant energy loss rate over the step
real_type dedx = physics_.make_calculator<EnergyLossCalculator>(
eloss_gid)(particle_.energy());
physics_.energy_loss_grid())(particle_.energy());

return particle_.energy() - Energy{step * dedx};
}
Expand Down
1 change: 0 additions & 1 deletion src/celeritas/phys/Model.hh
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
#include "celeritas/Types.hh"
#include "celeritas/global/ActionInterface.hh" // IWYU pragma: export
#include "celeritas/grid/ValueGridBuilder.hh"
#include "celeritas/grid/ValueGridType.hh"

#include "Applicability.hh" // IWYU pragma: export

Expand Down
29 changes: 13 additions & 16 deletions src/celeritas/phys/PhysicsData.hh
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,6 @@
#include "celeritas/em/data/AtomicRelaxationData.hh"
#include "celeritas/em/data/EPlusGGData.hh"
#include "celeritas/em/data/LivermorePEData.hh"
#include "celeritas/grid/ValueGridType.hh"
#include "celeritas/grid/XsGridData.hh"
#include "celeritas/neutron/data/NeutronElasticData.hh"

Expand Down Expand Up @@ -118,22 +117,20 @@ struct IntegralXsProcess
/*!
* Processes for a single particle type.
*
* Each index should be accessed with type ParticleProcessId. The "tables" are
* a fixed-size number of ItemRange references to ValueTables. The first index
* of the table (hard-coded) corresponds to ValueGridType; the second index is
* a ParticleProcessId. So the cross sections for ParticleProcessId{2} would
* be \code tables[ValueGridType::macro_xs][2] \endcode. This
* awkward access is encapsulated by the PhysicsTrackView. \c integral_xs will
* only be assigned if the integral approach is used and the particle has
* continuous-discrete processes.
* Each index should be accessed with type ParticleProcessId. \c macro_xs
* stores the cross section tables for each process, while \c energy_loss and
* \c range are the process-integrated dE/dx and range for the particle. c
amandalund marked this conversation as resolved.
Show resolved Hide resolved
* integral_xs will only be assigned if the integral approach is used and the
* particle has continuous-discrete processes.
*/
struct ProcessGroup
{
ItemRange<ProcessId> processes; //!< Processes that apply [ppid]
ValueGridArray<ItemRange<ValueTable>> tables; //!< [vgt][ppid]
ItemRange<IntegralXsProcess> integral_xs; //!< [ppid]
ItemRange<ModelGroup> models; //!< Model applicability [ppid]
ParticleProcessId eloss_ppid{}; //!< Process with de/dx and range tables
ItemRange<IntegralXsProcess> integral_xs; //!< [ppid]
ItemRange<ValueTable> macro_xs; //!< [ppid]
ValueTableId energy_loss; //!< Process-integrated energy loss
ValueTableId range; //!< Process-integrated range
bool has_at_rest{}; //!< Whether the particle type has an at-rest process

//! True if assigned and valid
Expand Down Expand Up @@ -306,16 +303,16 @@ struct PhysicsParamsScalars
/*!
* Persistent shared physics data.
*
* This includes macroscopic cross section, energy loss, and range tables
* ordered by [particle][process][material][energy].
* This includes macroscopic cross section tables ordered by
* [particle][process][material][energy] and process-integrated energy loss and
* range tables ordered by [particle][material][energy].
*
* So the first applicable process (ProcessId{0}) for an arbitrary particle
* (ParticleId{1}) in material 2 (MaterialId{2}) will have the following
* ID and cross section grid: \code
ProcessId proc_id = params.particle[1].processes[0];
const UniformGridData& grid
=
params.particle[1].table[int(ValueGridType::macro_xs)][0].material[2].log_energy;
= params.particle[1].macro_xs[0].material[2].log_energy;
* \endcode
*/
template<Ownership W, MemSpace M>
Expand Down
101 changes: 49 additions & 52 deletions src/celeritas/phys/PhysicsParams.cc
Original file line number Diff line number Diff line change
Expand Up @@ -454,6 +454,7 @@ void PhysicsParams::build_xs(Options const& opts,

using UPGridBuilder = Process::UPConstGridBuilder;
using Energy = Applicability::Energy;
using VGT = ValueGridType;

ValueGridInserter insert_grid(&data->reals, &data->value_grids);
auto value_tables = make_builder(&data->value_tables);
Expand All @@ -477,12 +478,11 @@ void PhysicsParams::build_xs(Options const& opts,
= data->model_groups[process_groups.models];
CELER_ASSERT(processes.size() == model_groups.size());

// Material-dependent physics tables, one per particle-process
ValueGridArray<std::vector<ValueTable>> temp_tables;
for (auto& vec : temp_tables)
{
vec.resize(processes.size());
}
// Material-dependent physics tables, one cross section table per
// particle-process and one dedx/range table per particle
std::vector<ValueTable> xs_table(processes.size());
ValueTable eloss_table;
ValueTable range_table;

// Processes with dE/dx and macro xs tables
std::vector<IntegralXsProcess> temp_integral_xs(processes.size());
Expand All @@ -501,11 +501,9 @@ void PhysicsParams::build_xs(Options const& opts,
Process const& proc = *this->process(processes[pp_idx]);

// Grid IDs for each grid type, each material
ValueGridArray<std::vector<ValueGridId>> temp_grid_ids;
for (auto& vec : temp_grid_ids)
{
vec.resize(mats.size());
}
std::vector<ValueGridId> xs_grid_ids(mats.size());
std::vector<ValueGridId> eloss_grid_ids(mats.size());
std::vector<ValueGridId> range_grid_ids(mats.size());

// Energy of maximum cross section for each material
std::vector<real_type> energy_max_xs;
Expand All @@ -517,9 +515,9 @@ void PhysicsParams::build_xs(Options const& opts,
}

// Loop over materials
for (auto mat_id : range(MaterialId{mats.size()}))
for (auto mat_idx : range(MaterialId::size_type{mats.size()}))
{
applic.material = mat_id;
applic.material = MaterialId(mat_idx);

// Construct step limit builders
auto builders = proc.step_limits(applic);
Expand All @@ -532,11 +530,10 @@ void PhysicsParams::build_xs(Options const& opts,
"have at least one)");

// Construct grids
for (auto vgt : range(ValueGridType::size_))
{
temp_grid_ids[vgt][mat_id.get()]
= build_grid(builders[vgt]);
}
xs_grid_ids[mat_idx] = build_grid(builders[VGT::macro_xs]);
eloss_grid_ids[mat_idx]
= build_grid(builders[VGT::energy_loss]);
range_grid_ids[mat_idx] = build_grid(builders[VGT::range]);

if (processes[pp_idx] == data->hardwired.positron_annihilation)
{
Expand All @@ -547,11 +544,10 @@ void PhysicsParams::build_xs(Options const& opts,
{
// Annihilation cross section is maximum at zero and
// decreases with increasing energy
energy_max_xs[mat_id.get()] = 0;
energy_max_xs[mat_idx] = 0;
}
}
else if (auto grid_id
= temp_grid_ids[ValueGridType::macro_xs][mat_id.get()])
else if (auto grid_id = xs_grid_ids[mat_idx])
{
auto const& grid_data = data->value_grids[grid_id];
auto data_ref = make_const_ref(*data);
Expand Down Expand Up @@ -579,37 +575,39 @@ void PhysicsParams::build_xs(Options const& opts,
}
}
CELER_ASSERT(e_max > 0);
energy_max_xs[mat_id.get()] = e_max;
energy_max_xs[mat_idx] = e_max;
}
}

// Index of the energy loss process that stores the de/dx and
// range tables
if (temp_grid_ids[ValueGridType::energy_loss][mat_id.get()]
&& temp_grid_ids[ValueGridType::range][mat_id.get()])
{
// Only one particle-process should have energy loss tables
CELER_ASSERT(!process_groups.eloss_ppid
|| pp_idx == process_groups.eloss_ppid.get());
process_groups.eloss_ppid = ParticleProcessId{pp_idx};
}
}

// Outer loop over grid types
for (auto vgt : range(ValueGridType::size_))
{
if (!std::any_of(temp_grid_ids[vgt].begin(),
temp_grid_ids[vgt].end(),
[](ValueGridId id) { return bool(id); }))
{
continue;
}
// Check if any material has value grids
auto has_grids = [](std::vector<ValueGridId> const& vec_id) {
return std::any_of(vec_id.begin(),
vec_id.end(),
[](ValueGridId id) { return bool(id); });
};

// Construct value grid table
ValueTable& temp_table = temp_tables[vgt][pp_idx];
temp_table.grids = value_grid_ids.insert_back(
temp_grid_ids[vgt].begin(), temp_grid_ids[vgt].end());
CELER_ASSERT(temp_table.grids.size() == mats.size());
// Construct value grid tables
if (has_grids(xs_grid_ids))
{
xs_table[pp_idx].grids = value_grid_ids.insert_back(
xs_grid_ids.begin(), xs_grid_ids.end());
CELER_ASSERT(xs_table[pp_idx].grids.size() == mats.size());
}
if (has_grids(eloss_grid_ids))
{
CELER_VALIDATE(!eloss_table && !range_table,
<< "more than one process for particle ID "
<< particle_id.get()
<< " has energy loss tables");

CELER_ASSERT(has_grids(range_grid_ids));
eloss_table.grids = value_grid_ids.insert_back(
eloss_grid_ids.begin(), eloss_grid_ids.end());
range_table.grids = value_grid_ids.insert_back(
range_grid_ids.begin(), range_grid_ids.end());
CELER_ASSERT(eloss_table.grids.size() == mats.size()
&& range_table.grids.size() == mats.size());
}

// Store the energies of the maximum cross sections
Expand All @@ -627,11 +625,10 @@ void PhysicsParams::build_xs(Options const& opts,
temp_integral_xs.begin(), temp_integral_xs.end());

// Construct value tables
for (auto vgt : range(ValueGridType::size_))
{
process_groups.tables[vgt] = value_tables.insert_back(
temp_tables[vgt].begin(), temp_tables[vgt].end());
}
process_groups.macro_xs
= value_tables.insert_back(xs_table.begin(), xs_table.end());
process_groups.energy_loss = value_tables.push_back(eloss_table);
process_groups.range = value_tables.push_back(range_table);
}
}

Expand Down
13 changes: 3 additions & 10 deletions src/celeritas/phys/PhysicsStepUtils.hh
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@
#include "celeritas/grid/InverseRangeCalculator.hh"
#include "celeritas/grid/RangeCalculator.hh"
#include "celeritas/grid/SplineXsCalculator.hh"
#include "celeritas/grid/ValueGridType.hh"
#include "celeritas/grid/XsCalculator.hh"
#include "celeritas/mat/MaterialTrackView.hh"
#include "celeritas/random/Selector.hh"
Expand All @@ -42,8 +41,6 @@ calc_physics_step_limit(MaterialTrackView const& material,
{
CELER_EXPECT(physics.has_interaction_mfp());

using VGT = ValueGridType;

/*! \todo For particles with decay, macro XS calculation will incorporate
* decay probability, dividing decay constant by speed to become 1/len to
* compete with interactions.
Expand Down Expand Up @@ -86,9 +83,8 @@ calc_physics_step_limit(MaterialTrackView const& material,
else
{
limit.step = physics.interaction_mfp() / total_macro_xs;
if (auto ppid = physics.eloss_ppid())
if (auto grid_id = physics.range_grid())
{
auto grid_id = physics.value_grid(VGT::range, ppid);
auto calc_range = physics.make_calculator<RangeCalculator>(grid_id);
real_type range = calc_range(particle.energy());
// Save range for the current step and reuse it elsewhere
Expand Down Expand Up @@ -180,20 +176,17 @@ calc_mean_energy_loss(ParticleTrackView const& particle,
real_type step)
{
CELER_EXPECT(step > 0);
CELER_EXPECT(physics.eloss_ppid());
using Energy = ParticleTrackView::Energy;
using VGT = ValueGridType;
static_assert(Energy::unit_type::value()
== EnergyLossCalculator::Energy::unit_type::value(),
"Incompatible energy types");

auto ppid = physics.eloss_ppid();
Energy const pre_step_energy = particle.energy();

// Calculate the sum of energy loss rate over all processes.
Energy eloss;
{
auto grid_id = physics.value_grid(VGT::energy_loss, ppid);
auto grid_id = physics.energy_loss_grid();
CELER_ASSERT(grid_id);

size_type order = physics.scalars().spline_eloss_order;
Expand All @@ -218,7 +211,7 @@ calc_mean_energy_loss(ParticleTrackView const& particle,
// approximation is probably wrong. Use the definition of the range as
// the integral of 1/loss to back-calculate the actual energy loss
// along the curve given the actual step.
auto grid_id = physics.value_grid(VGT::range, ppid);
auto grid_id = physics.range_grid();
CELER_ASSERT(grid_id);

// Use the range limit stored from calc_physics_step_limit
Expand Down
Loading
Loading