Skip to content

Commit

Permalink
Refactor MSC params and add Wentzel VI params and data (celeritas-pro…
Browse files Browse the repository at this point in the history 
…ject#1129)

* Add Wentzel VI data

- Move common MSC data into separate structs
- Import geom_factor
- Check imported Wentzel VI data in test

* Add Wentzel VI params and MSC params interface

* Move common MSC params to physics

* Add helper class for constructing MSC params and address other feedback

* Remove geometry factor (currently unused)

* Fix test for float

* Change fact --> factor

* Add min/max energy methods to XsCalculator
  • Loading branch information
@amandalund
amandalund committed Mar 17, 2024
1 parent 291393e commit d600808
Show file tree
Hide file tree
Showing 23 changed files with 761 additions and 322 deletions.
2 changes: 2 additions & 0 deletions src/celeritas/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,8 @@ list(APPEND SOURCES
em/AtomicRelaxationParams.cc
em/FluctuationParams.cc
em/UrbanMscParams.cc
em/WentzelVIMscParams.cc
em/detail/MscParamsHelper.cc
em/detail/Utils.cc
em/process/BremsstrahlungProcess.cc
em/process/ComptonProcess.cc
Expand Down
245 changes: 59 additions & 186 deletions src/celeritas/em/UrbanMscParams.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,20 +22,19 @@
#include "corecel/sys/ScopedMem.hh"
#include "celeritas/Quantities.hh"
#include "celeritas/grid/PolyEvaluator.hh"
#include "celeritas/grid/ValueGridBuilder.hh"
#include "celeritas/grid/ValueGridInserter.hh"
#include "celeritas/grid/XsGridData.hh"
#include "celeritas/grid/XsCalculator.hh"
#include "celeritas/io/ImportData.hh"
#include "celeritas/io/ImportProcess.hh"
#include "celeritas/mat/MaterialParams.hh"
#include "celeritas/mat/MaterialView.hh"
#include "celeritas/phys/ImportedProcessAdapter.hh"
#include "celeritas/phys/PDGNumber.hh"
#include "celeritas/phys/ParticleParams.hh"
#include "celeritas/phys/ParticleView.hh"

#include "data/UrbanMscData.hh"

#include "detail/MscParamsHelper.hh"

namespace celeritas
{
//---------------------------------------------------------------------------//
Expand All @@ -45,26 +44,18 @@ namespace celeritas
std::shared_ptr<UrbanMscParams>
UrbanMscParams::from_import(ParticleParams const& particles,
MaterialParams const& materials,
ImportData const& import)
ImportData const& data)
{
auto is_urban = [](ImportMscModel const& imm) {
return imm.model_class == ImportModelClass::urban_msc;
};
if (!std::any_of(
import.msc_models.begin(), import.msc_models.end(), is_urban))
if (!std::any_of(data.msc_models.begin(), data.msc_models.end(), is_urban))
{
// No Urban MSC present
return nullptr;
}

Options opts;
opts.lambda_limit = import.em_params.msc_lambda_limit;
opts.safety_fact = import.em_params.msc_safety_factor;
opts.range_fact = import.em_params.msc_range_factor;
opts.step_limit_algorithm = import.em_params.msc_step_algorithm;

return std::make_shared<UrbanMscParams>(
particles, materials, import.msc_models, opts);
particles, materials, data.msc_models);
}

//---------------------------------------------------------------------------//
Expand All @@ -73,202 +64,84 @@ UrbanMscParams::from_import(ParticleParams const& particles,
*/
UrbanMscParams::UrbanMscParams(ParticleParams const& particles,
MaterialParams const& materials,
VecImportMscModel const& mdata_vec,
Options options)
VecImportMscModel const& mdata_vec)
{
using units::MevEnergy;

ScopedMem record_mem("UrbanMscParams.construct");

HostVal<UrbanMscData> host_data;

host_data.ids.electron = particles.find(pdg::electron());
host_data.ids.positron = particles.find(pdg::positron());
CELER_VALIDATE(host_data.ids.electron && host_data.ids.positron,
<< "missing e-/e+ (required for Urban MSC)");
MscParamsHelper helper(
particles, materials, mdata_vec, ImportModelClass::urban_msc);
helper.build_ids(&host_data.ids);
helper.build_xs(&host_data.xs, &host_data.reals);

// Save electron mass
host_data.electron_mass = particles.get(host_data.ids.electron).mass();

// TODO: change IDs to a vector for all particles. This model should apply
// to muons and charged hadrons as well
if (particles.find(pdg::mu_minus()) || particles.find(pdg::mu_plus())
|| particles.find(pdg::proton()))
{
CELER_LOG(warning) << "Multiple scattering is not implemented for for "
"particles other than electron and positron";
}
// Coefficients for scaled Z
static Array<double, 2> const a_coeff{{0.87, 0.70}};
static Array<double, 2> const b_coeff{{2.0 / 3, 1.0 / 2}};

// Save parameters
CELER_VALIDATE(options.lambda_limit > 0,
<< "invalid lambda_limit=" << options.lambda_limit
<< " (should be positive)");
CELER_VALIDATE(options.safety_fact >= 0.1,
<< "invalid safety_fact=" << options.safety_fact
<< " (should be >= 0.1)");
CELER_VALIDATE(options.range_fact > 0 && options.range_fact < 1,
<< "invalid range_fact=" << options.range_fact
<< " (should be within 0 < limit < 1)");
host_data.params.lambda_limit = options.lambda_limit;
host_data.params.range_fact = options.range_fact;
host_data.params.safety_fact = options.safety_fact;
host_data.params.step_limit_algorithm = options.step_limit_algorithm;
if (host_data.params.step_limit_algorithm
== MscStepLimitAlgorithm::distance_to_boundary)
{
CELER_LOG(warning) << "Unsupported step limit algorithm '"
<< to_cstring(host_data.params.step_limit_algorithm)
<< "': defaulting to '"
<< to_cstring(MscStepLimitAlgorithm::safety) << "'";
host_data.params.step_limit_algorithm = MscStepLimitAlgorithm::safety;
}
// Builders
auto mdata = make_builder(&host_data.material_data);
auto pmdata = make_builder(&host_data.par_mat_data);
mdata.reserve(materials.num_materials());
pmdata.reserve(2 * materials.num_materials());

// Filter MSC data by model and particle type
std::vector<ImportMscModel const*> urban_data(particles.size(), nullptr);
for (ImportMscModel const& imm : mdata_vec)
for (auto mat_id : range(MaterialId{materials.num_materials()}))
{
// Filter out other MSC models
if (imm.model_class != ImportModelClass::urban_msc)
continue;
auto&& mat = materials.get(mat_id);

// Filter out unused particles
PDGNumber pdg{imm.particle_pdg};
ParticleId pid = pdg ? particles.find(pdg) : ParticleId{};
if (!pid)
continue;
// Build material-dependent data
mdata.push_back(UrbanMscParams::calc_material_data(mat));

if (!urban_data[pid.get()])
{
// Save data
urban_data[pid.get()] = &imm;
}
else
// Build particle-dependent data
Array<ParticleId, 2> const par_ids{{particles.find(pdg::electron()),
particles.find(pdg::positron())}};
double const zeff = mat.zeff();
for (size_type p : range(par_ids.size()))
{
// Warn: possibly multiple physics lists or different models in
// different regions
CELER_LOG(warning)
<< "duplicate " << to_cstring(imm.model_class)
<< " physics data for particle " << particles.id_to_label(pid)
<< ": ignoring all but the first encountered model";
}
}

auto get_scaled_xs = [&urban_data, &particles](ParticleId pid) {
CELER_ASSERT(pid < urban_data.size());
ImportMscModel const* imm = urban_data[pid.unchecked_get()];
CELER_VALIDATE(imm,
<< "missing Urban MSC physics data for particle "
<< particles.id_to_label(pid));
return &imm->xs_table;
};
UrbanMscParMatData this_pm;

{
// Set initial high/low energy limits
auto const& phys_vec
= get_scaled_xs(host_data.ids.electron)->physics_vectors;
CELER_ASSERT(!phys_vec.empty());
CELER_ASSERT(!phys_vec[0].x.empty());
host_data.params.low_energy_limit = MevEnergy(phys_vec[0].x.front());
host_data.params.high_energy_limit = MevEnergy(phys_vec[0].x.back());
}

{
// Particle-dependent data
Array<ParticleId, 2> const par_ids{
{host_data.ids.electron, host_data.ids.positron}};
Array<ImportPhysicsTable const*, 2> const xs_tables{{
get_scaled_xs(par_ids[0]),
get_scaled_xs(par_ids[1]),
}};
CELER_ASSERT(xs_tables[0]->x_units == ImportUnits::mev);
CELER_ASSERT(xs_tables[0]->y_units == ImportUnits::mev_2_per_cm);
CELER_ASSERT(xs_tables[1]->x_units == ImportUnits::mev);
CELER_ASSERT(xs_tables[1]->y_units == ImportUnits::mev_2_per_cm);
// Calculate scaled zeff
this_pm.scaled_zeff = a_coeff[p] * fastpow(zeff, b_coeff[p]);

// Coefficients for scaled Z
static Array<double, 2> const a_coeff{{0.87, 0.70}};
static Array<double, 2> const b_coeff{{2.0 / 3, 1.0 / 2}};

// Builders
auto mdata = make_builder(&host_data.material_data);
auto pmdata = make_builder(&host_data.par_mat_data);
mdata.reserve(materials.num_materials());
pmdata.reserve(2 * materials.num_materials());

// TODO: simplify when refactoring ValueGridInserter, etc
ValueGridInserter::XsGridCollection xgc;
ValueGridInserter vgi{&host_data.reals, &xgc};

for (auto mat_id : range(MaterialId{materials.num_materials()}))
{
auto&& mat = materials.get(mat_id);

// Build material-dependent data
mdata.push_back(UrbanMscParams::calc_material_data(mat));

// Build particle-dependent data
double const zeff = mat.zeff();
for (size_type p : range(par_ids.size()))
// Compute the maximum distance that particles can travel
// (different for electrons, hadrons)
if (par_ids[p] == host_data.ids.electron
|| par_ids[p] == host_data.ids.positron)
{
UrbanMscParMatData this_pm;

// Calculate scaled zeff
this_pm.scaled_zeff = a_coeff[p] * fastpow(zeff, b_coeff[p]);

// Compute the maximum distance that particles can travel
// (different for electrons, hadrons)
if (par_ids[p] == host_data.ids.electron
|| par_ids[p] == host_data.ids.positron)
{
// Electrons and positrons
this_pm.d_over_r = 9.6280e-1 - 8.4848e-2 * std::sqrt(zeff)
+ 4.3769e-3 * zeff;
CELER_ASSERT(0 < this_pm.d_over_r && this_pm.d_over_r <= 1);
}
else
{
// Muons and charged hadrons
this_pm.d_over_r = 1.15 - 9.76e-4 * zeff;
CELER_ASSERT(0 < this_pm.d_over_r);
}

// Get the cross section data for this particle and material
ImportPhysicsVector const& pvec
= xs_tables[p]->physics_vectors[mat_id.unchecked_get()];
CELER_ASSERT(pvec.vector_type == ImportPhysicsVectorType::log);

// Check that the limits are the same for all materials and
// particles; otherwise we need to change
// `UrbanMsc::is_applicable` to look up the particle and
// material
CELER_VALIDATE(host_data.params.low_energy_limit
== MevEnergy(pvec.x.front())
&& host_data.params.high_energy_limit
== MevEnergy(pvec.x.back()),
<< "multiple scattering cross section energy "
"limits are inconsistent across particles "
"and/or materials");

// To reuse existing code (TODO: simplify when refactoring)
// use the value grid builder to construct the grid entry in a
// temporary container and then copy it into the pm data.
auto vgb = ValueGridLogBuilder::from_geant(make_span(pvec.x),
make_span(pvec.y));
auto grid_id = vgb->build(vgi);
CELER_ASSERT(grid_id.get() == pmdata.size());
this_pm.xs = xgc[grid_id];
pmdata.push_back(this_pm);
CELER_ASSERT(host_data.at(mat_id, par_ids[p]).get() + 1
== host_data.par_mat_data.size());
// Electrons and positrons
this_pm.d_over_r = 9.6280e-1 - 8.4848e-2 * std::sqrt(zeff)
+ 4.3769e-3 * zeff;
CELER_ASSERT(0 < this_pm.d_over_r && this_pm.d_over_r <= 1);
}
else
{
// Muons and charged hadrons
this_pm.d_over_r = 1.15 - 9.76e-4 * zeff;
CELER_ASSERT(0 < this_pm.d_over_r);
}
pmdata.push_back(this_pm);
CELER_ASSERT(
host_data.at<UrbanMscParMatData>(mat_id, par_ids[p]).get() + 1
== host_data.par_mat_data.size());
}
CELER_ASSERT(host_data);
}

// Move to mirrored data, copying to device
mirror_ = CollectionMirror<UrbanMscData>{std::move(host_data)};
// Save high/low energy limits
XsCalculator calc_xs(host_data.xs[ItemId<XsGridData>(0)],
make_const_ref(host_data).reals);
host_data.params.low_energy_limit = calc_xs.energy_min();
host_data.params.high_energy_limit = calc_xs.energy_max();

CELER_ENSURE(this->mirror_);
CELER_ASSERT(host_data);

// Move to mirrored data, copying to device
data_ = CollectionMirror<UrbanMscData>{std::move(host_data)};
CELER_ENSURE(data_);
}

//---------------------------------------------------------------------------//
Expand Down
30 changes: 7 additions & 23 deletions src/celeritas/em/UrbanMscParams.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@
#pragma once

#include <memory>
#include <vector>

#include "corecel/data/CollectionMirror.hh"
#include "corecel/data/ParamsDataInterface.hh"
Expand All @@ -29,10 +28,6 @@ struct ImportMscModel;
* Construct and store data for Urban multiple scattering.
*
* Multiple scattering is used by the along-step kernel(s).
*
* TODO: UrbanMsc is the only MSC model presently in Celeritas, but we should
* extend this to be an interface since it's an interchangeable component of
* the along-step kernel(s).
*/
class UrbanMscParams final : public ParamsDataInterface<UrbanMscData>
{
Expand All @@ -42,41 +37,30 @@ class UrbanMscParams final : public ParamsDataInterface<UrbanMscData>
using VecImportMscModel = std::vector<ImportMscModel>;
//!@}

//! MSC configuration options
struct Options
{
real_type lambda_limit{1 * units::millimeter}; //!< Lambda limit
real_type range_fact{0.04}; //!< Range factor for e-/e+
real_type safety_fact{0.6}; //!< Safety factor
MscStepLimitAlgorithm step_limit_algorithm{
MscStepLimitAlgorithm::safety};
};

public:
// Construct if MSC process data is present, else return nullptr
static std::shared_ptr<UrbanMscParams>
from_import(ParticleParams const& particles,
MaterialParams const& materials,
ImportData const& import);
ImportData const& data);

// Construct from process data
inline UrbanMscParams(ParticleParams const& particles,
MaterialParams const& materials,
VecImportMscModel const& mdata,
Options options);
UrbanMscParams(ParticleParams const& particles,
MaterialParams const& materials,
VecImportMscModel const& mdata);

// TODO: possible "applicability" interface used for constructing
// along-step kernels?

//! Access UrbanMsc data on the host
HostRef const& host_ref() const final { return mirror_.host_ref(); }
HostRef const& host_ref() const final { return data_.host_ref(); }

//! Access UrbanMsc data on the device
DeviceRef const& device_ref() const final { return mirror_.device_ref(); }
DeviceRef const& device_ref() const final { return data_.device_ref(); }

private:
// Host/device storage and reference
CollectionMirror<UrbanMscData> mirror_;
CollectionMirror<UrbanMscData> data_;

static UrbanMscMaterialData
calc_material_data(MaterialView const& material_view);
Expand Down

0 comments on commit d600808

Please sign in to comment.