Skip to content

Commit

Permalink
Support combined single and multiple Coulomb scattering (#1230)
Browse files Browse the repository at this point in the history
* Start adding support for combined single and multiple Coulomb scattering

* Import polar angle limit and check for combined SS and MSC

* Move material properties to material data and separate common SS and MSC data

* Refactor Coulomb scattering data

- Add a new data/params class for Wentzel OK&VI data that will be used by both
  the single Coulomb scattering model and the Wentzel VI MSC model
- Import MscThetaLimit from EM parameters instead of the model PolarAngleLimit

* Update material test

* Update Wentzel distribution to support arbitrary cos theta min/max

* Check kin factor and cross sections in Wentzel helper test

* Defer WentzelVI cross section calculator

* Get mott factor from helper, finish importing angle limit factor, and clean up code

* Import nuclear form factor model from Geant4

* a_sq_factor can be zero with float

* Change (1 - generate_canonical) --> generate_canonical in Wentzel distribution and update tests
  • Loading branch information
amandalund authored May 23, 2024
1 parent 2da4678 commit 65ed16e
Show file tree
Hide file tree
Showing 50 changed files with 1,921 additions and 1,315 deletions.
4 changes: 4 additions & 0 deletions app/celer-sim/Runner.cc
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@
#include "celeritas/Types.hh"
#include "celeritas/Units.hh"
#include "celeritas/em/params/UrbanMscParams.hh"
#include "celeritas/em/params/WentzelOKVIParams.hh"
#include "celeritas/ext/GeantImporter.hh"
#include "celeritas/ext/GeantSetup.hh"
#include "celeritas/ext/RootFileManager.hh"
Expand Down Expand Up @@ -308,6 +309,9 @@ void Runner::build_core_params(RunnerInput const& inp,
params.cutoff = CutoffParams::from_import(
imported, params.particle, params.material);

// Construct shared data for Coulomb scattering
params.wentzel = WentzelOKVIParams::from_import(imported, params.material);

// Load physics: create individual processes with make_shared
params.physics = [&params, &inp, &imported] {
PhysicsParams::Input input;
Expand Down
5 changes: 5 additions & 0 deletions src/accel/SharedParams.cc
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include "geocel/GeantUtils.hh"
#include "geocel/g4/GeantGeoParams.hh"
#include "celeritas/Types.hh"
#include "celeritas/em/params/WentzelOKVIParams.hh"
#include "celeritas/ext/GeantImporter.hh"
#include "celeritas/ext/RootExporter.hh"
#include "celeritas/geo/GeoMaterialParams.hh"
Expand All @@ -58,6 +59,7 @@

#include "AlongStepFactory.hh"
#include "SetupOptions.hh"

#include "detail/HitManager.hh"
#include "detail/OffloadWriter.hh"

Expand Down Expand Up @@ -530,6 +532,9 @@ void SharedParams::initialize_core(SetupOptions const& options)
params.cutoff = CutoffParams::from_import(
*imported, params.particle, params.material);

// Construct shared data for Coulomb scattering
params.wentzel = WentzelOKVIParams::from_import(*imported, params.material);

// Load physics: create individual processes with make_shared
params.physics = [&params, &options, &imported] {
PhysicsParams::Input input;
Expand Down
1 change: 1 addition & 0 deletions src/celeritas/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ list(APPEND SOURCES
em/params/AtomicRelaxationParams.cc
em/params/FluctuationParams.cc
em/params/UrbanMscParams.cc
em/params/WentzelOKVIParams.cc
em/params/WentzelVIMscParams.cc
em/params/detail/MscParamsHelper.cc
em/process/BremsstrahlungProcess.cc
Expand Down
15 changes: 15 additions & 0 deletions src/celeritas/Types.cc
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,21 @@ char const* to_cstring(MscStepLimitAlgorithm value)
return to_cstring_impl(value);
}

//---------------------------------------------------------------------------//
/*!
* Get a string corresponding to the nuclear form factor model.
*/
char const* to_cstring(NuclearFormFactorType value)
{
static EnumStringMapper<NuclearFormFactorType> const to_cstring_impl{
"none",
"flat",
"exponential",
"gaussian",
};
return to_cstring_impl(value);
}

//---------------------------------------------------------------------------//
/*!
* Checks that the TrackOrder will sort tracks by actions applied at the given
Expand Down
14 changes: 14 additions & 0 deletions src/celeritas/Types.hh
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,17 @@ enum class MscStepLimitAlgorithm
size_,
};

//---------------------------------------------------------------------------//
//! Nuclear form factor model for Coulomb scattering
enum class NuclearFormFactorType
{
none,
flat,
exponential,
gaussian,
size_
};

//---------------------------------------------------------------------------//
// HELPER STRUCTS
//---------------------------------------------------------------------------//
Expand Down Expand Up @@ -196,6 +207,9 @@ char const* to_cstring(TrackOrder);
// Get a string corresponding to the MSC step limit algorithm
char const* to_cstring(MscStepLimitAlgorithm value);

// Get a string corresponding to the nuclear form factor model
char const* to_cstring(NuclearFormFactorType value);

// Checks that the TrackOrder will sort tracks by actions applied at the given
// ActionOrder
bool is_action_sorted(ActionOrder action, TrackOrder track);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file celeritas/em/data/MscData.hh
//! \file celeritas/em/data/CommonCoulombData.hh
//---------------------------------------------------------------------------//
#pragma once

Expand All @@ -19,7 +19,7 @@ namespace celeritas
*
* TODO these will probably be changed to a map over all particle IDs.
*/
struct MscIds
struct CoulombIds
{
ParticleId electron;
ParticleId positron;
Expand Down
107 changes: 9 additions & 98 deletions src/celeritas/em/data/CoulombScatteringData.hh
Original file line number Diff line number Diff line change
Expand Up @@ -9,120 +9,31 @@

#include "corecel/Macros.hh"
#include "corecel/Types.hh"
#include "corecel/data/Collection.hh"
#include "celeritas/Quantities.hh"
#include "celeritas/Types.hh"

namespace celeritas
{
//---------------------------------------------------------------------------//
/*!
* Particle and action ids used by CoulombScatteringModel.
*/
struct CoulombScatteringIds
{
ActionId action;
ParticleId electron;
ParticleId positron;
ParticleId proton;

explicit CELER_FUNCTION operator bool() const
{
// TODO: enable when protons are supported
return action && electron && positron /* && proton */;
}
};

//---------------------------------------------------------------------------//
/*!
* Per-element data used by the CoulombScatteringModel.
*
* The matrix of coefficients used to approximate the ratio of the Mott to
* Rutherford cross sections was developed in
* T. Lijian, H. Quing and L. Zhengming, Radiat. Phys. Chem. 45 (1995),
* 235-245
* and
* M. J. Boschini et al. arXiv:1111.4042
*/
struct CoulombScatteringElementData
{
//!@{
//! \name Dimensions for Mott coefficient matrices
static constexpr size_type num_mott_beta_bins = 6;
static constexpr size_type num_mott_theta_bins = 5;
static constexpr size_type num_mott_elements = 92;
//!@}

using BetaArray = Array<real_type, num_mott_beta_bins>;
using ThetaArray = Array<real_type, num_mott_theta_bins>;
using MottCoeffMatrix = Array<BetaArray, num_mott_theta_bins>;
#include "CommonCoulombData.hh"

//! Matrix of Mott coefficients [theta][beta]
MottCoeffMatrix mott_coeff;
};

//---------------------------------------------------------------------------//
/*!
* Supported models of nuclear form factors.
*/
enum class NuclearFormFactorType
namespace celeritas
{
none,
flat,
exponential,
gaussian
};

//---------------------------------------------------------------------------//
/*!
* Constant shared data used by the CoulombScatteringModel.
*/
template<Ownership W, MemSpace M>
struct CoulombScatteringData
{
template<class T>
using ElementItems = celeritas::Collection<T, W, M, ElementId>;
template<class T>
using IsotopeItems = celeritas::Collection<T, W, M, IsotopeId>;

// Ids
CoulombScatteringIds ids;

//! Constant prefactor for the squared momentum transfer [(MeV/c)^-2]
IsotopeItems<real_type> nuclear_form_prefactor;

// Per element form factors
ElementItems<CoulombScatteringElementData> elem_data;
// Particle IDs
CoulombIds ids;

// User-defined factor for the screening coefficient
real_type screening_factor;
// Action ID
ActionId action;

// Model for the form factor to use
NuclearFormFactorType form_factor_type;
//! Cosine of the maximum scattering polar angle
static CELER_CONSTEXPR_FUNCTION real_type cos_thetamax() { return -1; }

// Check if the data is initialized
explicit CELER_FUNCTION operator bool() const
{
return ids && !nuclear_form_prefactor.empty() && !elem_data.empty();
}

// Copy initialize from an existing CoulombScatteringData
template<Ownership W2, MemSpace M2>
CoulombScatteringData& operator=(CoulombScatteringData<W2, M2> const& other)
{
CELER_EXPECT(other);
ids = other.ids;
nuclear_form_prefactor = other.nuclear_form_prefactor;
elem_data = other.elem_data;
form_factor_type = other.form_factor_type;
screening_factor = other.screening_factor;
return *this;
}
explicit CELER_FUNCTION operator bool() const { return ids && action; }
};

using CoulombScatteringDeviceRef = DeviceCRef<CoulombScatteringData>;
using CoulombScatteringHostRef = HostCRef<CoulombScatteringData>;
using CoulombScatteringRef = NativeCRef<CoulombScatteringData>;

//---------------------------------------------------------------------------//
} // namespace celeritas
4 changes: 2 additions & 2 deletions src/celeritas/em/data/UrbanMscData.hh
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
#include "celeritas/Types.hh"
#include "celeritas/grid/XsGridData.hh"

#include "MscData.hh"
#include "CommonCoulombData.hh"

namespace celeritas
{
Expand Down Expand Up @@ -137,7 +137,7 @@ struct UrbanMscData
//// DATA ////

//! Particle IDs
MscIds ids;
CoulombIds ids;
//! Mass of of electron in MeV
units::MevMass electron_mass;
//! User-assignable options
Expand Down
123 changes: 123 additions & 0 deletions src/celeritas/em/data/WentzelOKVIData.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,123 @@
//----------------------------------*-C++-*----------------------------------//
// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
// See the top-level COPYRIGHT file for details.
// SPDX-License-Identifier: (Apache-2.0 OR MIT)
//---------------------------------------------------------------------------//
//! \file celeritas/em/data/WentzelOKVIData.hh
//---------------------------------------------------------------------------//
#pragma once

#include "corecel/Macros.hh"
#include "corecel/Types.hh"
#include "corecel/data/Collection.hh"
#include "celeritas/Constants.hh"
#include "celeritas/Quantities.hh"
#include "celeritas/Types.hh"

namespace celeritas
{
//---------------------------------------------------------------------------//
/*!
* Parameters used in both single Coulomb scattering and Wentzel VI MSC models.
*
* When the single Coulomb scattering and Wentzel VI MSC models are used
* together, the MSC model is used to sample scatterings with angles below the
* polar angle limit, and the single scattering model is used for angles above
* the limit.
*/
struct CoulombParameters
{
//! Whether to use combined single and multiple scattering
bool is_combined{true};
//! Polar angle limit between single and multiple scattering
real_type costheta_limit{-1};
//! Factor for the screening coefficient
real_type screening_factor{1};
//! Factor used to calculate the maximum scattering angle off of a nucleus
real_type a_sq_factor{0.5
* ipow<2>(constants::hbar_planck * constants::c_light
* units::femtometer)};
// Model for the form factor to use
NuclearFormFactorType form_factor_type{NuclearFormFactorType::exponential};

explicit CELER_FUNCTION operator bool() const
{
return costheta_limit >= -1 && costheta_limit <= 1
&& screening_factor > 0 && a_sq_factor >= 0
&& form_factor_type != NuclearFormFactorType::size_;
}
};

//---------------------------------------------------------------------------//
/*!
* Per-element data used by the Coulomb scattering and Wentzel VI models.
*
* The matrix of coefficients used to approximate the ratio of the Mott to
* Rutherford cross sections was developed in T. Lijian, H. Quing and L.
* Zhengming, Radiat. Phys. Chem. 45 (1995), 235-245 and M. J. Boschini et al.
* arXiv:1111.4042
*/
struct MottElementData
{
//!@{
//! \name Dimensions for Mott coefficient matrices
static constexpr size_type num_mott_beta_bins = 6;
static constexpr size_type num_mott_theta_bins = 5;
static constexpr size_type num_mott_elements = 92;
//!@}

using BetaArray = Array<real_type, num_mott_beta_bins>;
using ThetaArray = Array<real_type, num_mott_theta_bins>;
using MottCoeffMatrix = Array<BetaArray, num_mott_theta_bins>;

//! Matrix of Mott coefficients [theta][beta]
MottCoeffMatrix mott_coeff;
};

//---------------------------------------------------------------------------//
/*!
* Constant shared data used by the Coulomb scattering and Wentzel VI models.
*/
template<Ownership W, MemSpace M>
struct WentzelOKVIData
{
template<class T>
using ElementItems = celeritas::Collection<T, W, M, ElementId>;
template<class T>
using IsotopeItems = celeritas::Collection<T, W, M, IsotopeId>;
template<class T>
using MaterialItems = Collection<T, W, M, MaterialId>;

// User-assignable parameters
CoulombParameters params;

// Constant prefactor for the squared momentum transfer [(MeV/c)^-2]
IsotopeItems<real_type> nuclear_form_prefactor;

// Per element form factors
ElementItems<MottElementData> elem_data;

// Inverse effective A^2/3 [1/mass^2/3]
MaterialItems<real_type> inv_mass_cbrt_sq;

// Check if the data is initialized
explicit CELER_FUNCTION operator bool() const
{
return params && !elem_data.empty()
&& params.is_combined == !inv_mass_cbrt_sq.empty();
}

template<Ownership W2, MemSpace M2>
WentzelOKVIData& operator=(WentzelOKVIData<W2, M2> const& other)
{
CELER_EXPECT(other);
params = other.params;
nuclear_form_prefactor = other.nuclear_form_prefactor;
elem_data = other.elem_data;
inv_mass_cbrt_sq = other.inv_mass_cbrt_sq;
return *this;
}
};

//---------------------------------------------------------------------------//
} // namespace celeritas
Loading

0 comments on commit 65ed16e

Please sign in to comment.