From 706496358a6bbbd2a5035932c9f0abdc8fd72f77 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 9 Feb 2021 16:28:23 -0500 Subject: [PATCH 1/6] Update spack scripts --- scripts/dev/env/celeritas-darwin.yaml | 5 +++-- scripts/dev/env/celeritas-linux.yaml | 1 + 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/scripts/dev/env/celeritas-darwin.yaml b/scripts/dev/env/celeritas-darwin.yaml index 5eae1fbd6e..5290bf48e3 100644 --- a/scripts/dev/env/celeritas-darwin.yaml +++ b/scripts/dev/env/celeritas-darwin.yaml @@ -2,6 +2,7 @@ spack: specs: - cmake + - doxygen - geant4 - git - git-lfs @@ -13,7 +14,7 @@ spack: - openmpi - python - - root +aqua + - root ~aqua - swig - veccore @@ -22,7 +23,7 @@ spack: concretization: together packages: root: - variants: ~davix ~examples ~x ~opengl ~tbb ~rootfit ~math ~gsl cxxstd=14 + variants: ~davix ~examples ~opengl ~x ~tbb ~rootfit ~math ~gsl cxxstd=14 all: providers: blas: [openblas] diff --git a/scripts/dev/env/celeritas-linux.yaml b/scripts/dev/env/celeritas-linux.yaml index 3cadb97cbe..59d33aba5c 100644 --- a/scripts/dev/env/celeritas-linux.yaml +++ b/scripts/dev/env/celeritas-linux.yaml @@ -2,6 +2,7 @@ spack: specs: - cmake - cuda + - doxygen - geant4 - git - git-lfs From 0d55a4f7d7081264f589839f40adf6a49a8d275d Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Tue, 9 Feb 2021 16:24:47 -0500 Subject: [PATCH 2/6] Add size accessors to material params --- src/physics/material/MaterialParams.hh | 6 ++++++ test/physics/material/Material.test.cc | 3 +++ 2 files changed, 9 insertions(+) diff --git a/src/physics/material/MaterialParams.hh b/src/physics/material/MaterialParams.hh index a6cc4536f0..236dd1df52 100644 --- a/src/physics/material/MaterialParams.hh +++ b/src/physics/material/MaterialParams.hh @@ -66,6 +66,12 @@ class MaterialParams // Construct with a vector of material definitions explicit MaterialParams(const Input& inp); + //! Number of material definitions + MaterialId::value_type size() const { return matnames_.size(); } + + //! Number of distinct elements definitions + ElementId::value_type num_elements() const { return elnames_.size(); } + // Get element name inline const std::string& id_to_label(ElementId id) const; diff --git a/test/physics/material/Material.test.cc b/test/physics/material/Material.test.cc index 6f638b5d04..eb07172c76 100644 --- a/test/physics/material/Material.test.cc +++ b/test/physics/material/Material.test.cc @@ -122,6 +122,9 @@ TEST_F(MaterialTest, params) { ASSERT_TRUE(params); + EXPECT_EQ(3, params->size()); + EXPECT_EQ(4, params->num_elements()); + EXPECT_EQ(MaterialId{0}, params->find("NaI")); EXPECT_EQ(MaterialId{1}, params->find("hard vacuum")); EXPECT_EQ(MaterialId{2}, params->find("H2")); From 10da811d72bc63db4fd95323579f0c7d15c3115e Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Feb 2021 17:33:40 -0500 Subject: [PATCH 3/6] Replace step limit builders with an array --- src/physics/base/Process.hh | 16 ++++++++-------- src/physics/em/ComptonProcess.cc | 2 +- test/physics/base/MockModel.cc | 2 +- test/physics/base/MockProcess.cc | 33 +++++++++++++++++++------------- test/physics/base/MockProcess.hh | 5 ++++- 5 files changed, 34 insertions(+), 24 deletions(-) diff --git a/src/physics/base/Process.hh b/src/physics/base/Process.hh index f8ef048487..d19be819aa 100644 --- a/src/physics/base/Process.hh +++ b/src/physics/base/Process.hh @@ -12,6 +12,7 @@ #include "Applicability.hh" #include "ModelIdGenerator.hh" #include "Types.hh" +#include "PhysicsInterface.hh" #include "physics/grid/ValueGridBuilder.hh" namespace celeritas @@ -30,6 +31,12 @@ class Model; * * Each process has an interaction ("post step doit") and may have both energy * loss and range limiters. + * + * The StepLimitBuilders is a fixed-size array corresponding to the physics + * interface enum \c ValueGridType : + * - macro_xs: Cross section [1/cm] + * - energy_loss: dE/dx [MeV/cm] + * - range: Range limit [cm] */ class Process { @@ -39,16 +46,9 @@ class Process using SPConstModel = std::shared_ptr; using UPConstGridBuilder = std::unique_ptr; using VecModel = std::vector; + using StepLimitBuilders = ValueGridArray; //!@} - //! Array builders for along-step energy loss and other quantities - struct StepLimitBuilders - { - UPConstGridBuilder macro_xs; //!< Cross section [1/cm] - UPConstGridBuilder energy_loss; //!< dE/dx [MeV/cm] - UPConstGridBuilder range; //!< Range limit [cm] - }; - public: // Virtual destructor for polymorphic deletion virtual ~Process() = 0; diff --git a/src/physics/em/ComptonProcess.cc b/src/physics/em/ComptonProcess.cc index fcc058672e..6f98d34a13 100644 --- a/src/physics/em/ComptonProcess.cc +++ b/src/physics/em/ComptonProcess.cc @@ -58,7 +58,7 @@ auto ComptonProcess::step_limits(Applicability range) const -> StepLimitBuilders // refined energies, switch based on the given applicability or hope that // the input has already preprocessed the energy dependence. StepLimitBuilders builders; - builders.macro_xs + builders[size_type(ValueGridType::macro_xs)] = std::make_unique(ValueGridXsBuilder::from_geant( make_span(lo.x), make_span(lo.y), make_span(hi.x), make_span(hi.y))); return builders; diff --git a/test/physics/base/MockModel.cc b/test/physics/base/MockModel.cc index 101c56a32f..c2d40466e5 100644 --- a/test/physics/base/MockModel.cc +++ b/test/physics/base/MockModel.cc @@ -34,7 +34,7 @@ void MockModel::interact(const ModelInteractPointers&) const std::string MockModel::label() const { std::ostringstream os; - os << "MockModel(p=" << applic_.particle.get() + os << "MockModel(" << id_.get() << ", p=" << applic_.particle.get() << ", emin=" << applic_.lower.value() << ", emax=" << applic_.upper.value() << ")"; return os.str(); diff --git a/test/physics/base/MockProcess.cc b/test/physics/base/MockProcess.cc index d05713e013..35bc4a6243 100644 --- a/test/physics/base/MockProcess.cc +++ b/test/physics/base/MockProcess.cc @@ -18,9 +18,10 @@ namespace celeritas_test MockProcess::MockProcess(Input data) : data_(std::move(data)) { CELER_EXPECT(data_.materials); + CELER_EXPECT(!data_.label.empty()); CELER_EXPECT(!data_.applic.empty()); CELER_EXPECT(data_.interact); - CELER_EXPECT(data_.xs >= 0); + CELER_EXPECT(data_.xs >= celeritas::zero_quantity()); CELER_EXPECT(data_.energy_loss >= 0); CELER_EXPECT(data_.range >= 0); } @@ -49,24 +50,30 @@ auto MockProcess::step_limits(Applicability range) const -> StepLimitBuilders real_type numdens = mat.number_density(); StepLimitBuilders builders; - if (data_.xs > 0) + if (data_.xs > zero_quantity()) { - real_type value = data_.xs * numdens; - builders.macro_xs = std::make_unique( - range.lower.value(), range.upper.value(), VecReal{value, value}); + real_type value = unit_cast(data_.xs) * numdens; + builders[size_type(ValueGridType::macro_xs)] + = std::make_unique(range.lower.value(), + range.upper.value(), + VecReal{value, value}); } if (data_.energy_loss > 0) { - real_type value = data_.energy_loss * numdens; - builders.energy_loss = std::make_unique( - range.lower.value(), range.upper.value(), VecReal{value, value}); + real_type value = data_.energy_loss * numdens; + builders[size_type(ValueGridType::energy_loss)] + = std::make_unique(range.lower.value(), + range.upper.value(), + VecReal{value, value}); } if (data_.range > 0) { - builders.range = std::make_unique( - range.lower.value(), - range.upper.value(), - VecReal{range.lower.value(), range.upper.value()}); + builders[size_type(ValueGridType::range)] + = std::make_unique( + range.lower.value(), + range.upper.value(), + VecReal{data_.range * range.lower.value(), + data_.range * range.upper.value()}); } return builders; @@ -75,7 +82,7 @@ auto MockProcess::step_limits(Applicability range) const -> StepLimitBuilders //---------------------------------------------------------------------------// std::string MockProcess::label() const { - return "mock"; + return data_.label; } //---------------------------------------------------------------------------// diff --git a/test/physics/base/MockProcess.hh b/test/physics/base/MockProcess.hh index 7c6deece6e..11b9c73b5b 100644 --- a/test/physics/base/MockProcess.hh +++ b/test/physics/base/MockProcess.hh @@ -12,6 +12,7 @@ #include #include #include "physics/base/Model.hh" +#include "physics/base/Units.hh" #include "physics/material/MaterialParams.hh" namespace celeritas_test @@ -37,6 +38,7 @@ class MockProcess : public celeritas::Process //!@{ //! Type aliases using real_type = celeritas::real_type; + using BarnMicroXs = celeritas::Quantity; using Applicability = celeritas::Applicability; using ModelIdGenerator = celeritas::ModelIdGenerator; using VecApplicability = std::vector; @@ -47,9 +49,10 @@ class MockProcess : public celeritas::Process struct Input { SPConstMaterials materials; + std::string label; VecApplicability applic; //!< Applicablity per model ModelCallback interact; //!< MockModel::interact callback - real_type xs{}; //!< Constant per atom [cm^2] + BarnMicroXs xs{}; //!< Constant per atom [bn] real_type energy_loss{}; //!< Constant per atom [MeV/cm / (1/cm^3)] real_type range{}; //!< Linearly increasing [cm/MeV] }; From d770ad9833111914dddce4bd03fcc35d4c110593 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Feb 2021 17:34:37 -0500 Subject: [PATCH 4/6] Implement backend for general physics operations --- src/physics/base/Applicability.hh | 13 +- src/physics/base/ModelInterface.hh | 4 +- src/physics/base/PhysicsInterface.hh | 136 ++++++-- src/physics/base/PhysicsParams.cc | 319 +++++++++++++++++- src/physics/base/PhysicsParams.hh | 99 +++--- src/physics/base/PhysicsStepUtils.i.hh | 17 +- src/physics/base/PhysicsTrackView.hh | 24 +- src/physics/base/PhysicsTrackView.i.hh | 63 ++-- src/physics/base/Types.hh | 2 +- src/physics/grid/ValueGridBuilder.cc | 21 +- src/physics/grid/ValueGridBuilder.hh | 17 +- test/CMakeLists.txt | 4 +- test/physics/base/Physics.test.cc | 375 ++++++++++++++++----- test/physics/base/PhysicsStepUtils.test.cc | 75 +++++ test/physics/base/PhysicsTestBase.cc | 159 +++++++++ test/physics/base/PhysicsTestBase.hh | 92 +++++ 16 files changed, 1200 insertions(+), 220 deletions(-) create mode 100644 test/physics/base/PhysicsStepUtils.test.cc create mode 100644 test/physics/base/PhysicsTestBase.cc create mode 100644 test/physics/base/PhysicsTestBase.hh diff --git a/src/physics/base/Applicability.hh b/src/physics/base/Applicability.hh index d7cf7c696f..5b258c2565 100644 --- a/src/physics/base/Applicability.hh +++ b/src/physics/base/Applicability.hh @@ -34,10 +34,13 @@ namespace celeritas */ struct Applicability { - MaterialId material{}; - ParticleId particle{}; - units::MevEnergy lower = zero_quantity(); - units::MevEnergy upper = max_quantity(); + using EnergyUnits = units::Mev; + using Energy = Quantity; + + MaterialId material{}; + ParticleId particle{}; + Energy lower = zero_quantity(); + Energy upper = max_quantity(); //! Range for a particle at rest static inline Applicability at_rest(ParticleId id) @@ -53,7 +56,7 @@ struct Applicability //! Whether applicability is in a valid state inline explicit operator bool() const { - return static_cast(particle); + return static_cast(particle) && lower < upper; } }; diff --git a/src/physics/base/ModelInterface.hh b/src/physics/base/ModelInterface.hh index c822d83b0d..2c333df64f 100644 --- a/src/physics/base/ModelInterface.hh +++ b/src/physics/base/ModelInterface.hh @@ -26,7 +26,7 @@ struct ModelInteractParams { ParticleParamsData particle; MaterialParamsData material; - PhysicsParamsPointers physics; + PhysicsParamsData physics; //! True if valid CELER_FUNCTION operator bool() const @@ -46,7 +46,7 @@ struct ModelInteractState { ParticleStateData particle; MaterialStateData material; - PhysicsStatePointers physics; + PhysicsStateData physics; Span direction; RngStatePointers rng; diff --git a/src/physics/base/PhysicsInterface.hh b/src/physics/base/PhysicsInterface.hh index 5cda14e9e6..016dcfa6ab 100644 --- a/src/physics/base/PhysicsInterface.hh +++ b/src/physics/base/PhysicsInterface.hh @@ -7,24 +7,31 @@ //---------------------------------------------------------------------------// #pragma once -#include "base/Span.hh" +#include "base/Array.hh" +#include "base/Pie.hh" #include "Types.hh" #include "physics/grid/XsGridInterface.hh" #include "physics/em/detail/LivermorePE.hh" #include "physics/em/detail/EPlusGG.hh" #include "physics/material/Types.hh" +#ifndef __CUDA_ARCH__ +# include "base/PieBuilder.hh" +#endif + namespace celeritas { //---------------------------------------------------------------------------// +// TYPES +//---------------------------------------------------------------------------// //! Currently all value grids are cross section grids -using ValueGrid = XsGridData; +using ValueGrid = XsGridData; +using ValueGridId = OpaqueId; +using ValueTableId = OpaqueId; -//---------------------------------------------------------------------------// -// PARAMS //---------------------------------------------------------------------------// //! Hardcoded types of grid data -enum class PhysicsTableType +enum class ValueGridType { macro_xs, //!< Interaction cross sections energy_loss, //!< Energy loss per unit length @@ -32,6 +39,11 @@ enum class PhysicsTableType size_ //!< Sentinel value }; +template +using ValueGridArray = Array; + +//---------------------------------------------------------------------------// +// PARAMS //---------------------------------------------------------------------------// /*! * Energy-dependent model IDs for a single process and particle type. @@ -43,8 +55,8 @@ enum class PhysicsTableType */ struct ModelGroup { - Span energy; //!< Energy grid bounds [MeV] - Span model; //!< Corresponding models + PieSlice energy; //!< Energy grid bounds [MeV] + PieSlice model; //!< Corresponding models //! True if assigned explicit CELER_FUNCTION operator bool() const @@ -64,7 +76,7 @@ struct ModelGroup */ struct ValueTable { - Span material; //!< Value grid by material index + PieSlice material; //!< Value grid by material index //! True if assigned explicit CELER_FUNCTION operator bool() const { return !material.empty(); } @@ -75,18 +87,17 @@ struct ValueTable * Processes for a single particle type. * * Each index should be accessed with type ParticleProcessId. The "tables" are - * a fixed-size number of Span references to ValueTables. The first index of - * the table (hard-coded) corresponds to PhysicsTableType; the second index is + * a fixed-size number of PieSlice 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[size_type(PhysicsTableType::macro_xs)][2] \endcode. This + * be \code tables[size_type(ValueGridType::macro_xs)][2] \endcode. This * awkward access is encapsulated by the PhysicsTrackView. */ struct ProcessGroup { - Span processes; //!< Processes that apply - - Array, size_type(PhysicsTableType::size_)> tables; //!< Data - Span models; //!< Model applicability + PieSlice processes; //!< Processes that apply [ppid] + ValueGridArray> tables; //!< [vgt][ppid] + PieSlice models; //!< Model applicability [ppid] //! True if assigned and valid explicit CELER_FUNCTION operator bool() const @@ -126,24 +137,66 @@ struct HardwiredModels ProcessId proc_id = params.particle[1].processes[0]; const UniformGridData& grid = - params.particle[1].table[int(PhysicsTableType::macro_xs)][0].material[2].log_energy; + params.particle[1].table[int(ValueGridType::macro_xs)][0].material[2].log_energy; * \endcode */ -struct PhysicsParamsPointers +template +struct PhysicsParamsData { - Span particle; - HardwiredModels hardwired; - size_type max_particle_processes{}; + template + using Data = Pie; + template + using ParticleData = Pie; + + // Backend storage + Data reals; + Data model_ids; + Data value_grids; + Data value_grid_ids; + Data process_ids; + Data value_tables; + Data model_groups; + ParticleData process_groups; + + HardwiredModels hardwired; + ProcessId::value_type max_particle_processes{}; //// USER-CONFIGURABLE CONSTANTS //// real_type scaling_min_range{}; //!< rho [cm] real_type scaling_fraction{}; //!< alpha [unitless] // real_type max_eloss_fraction{}; //!< For scaled range calculation + //// MEMBER FUNCTIONS //// + //! True if assigned explicit CELER_FUNCTION operator bool() const { - return !particle.empty() && max_particle_processes; + return !process_groups.empty() && max_particle_processes + && scaling_min_range > 0 && scaling_fraction > 0; + } + + //! Assign from another set of data + template + PhysicsParamsData& operator=(const PhysicsParamsData& other) + { + CELER_EXPECT(other); + + reals = other.reals; + model_ids = other.model_ids; + value_grids = other.value_grids; + value_grid_ids = other.value_grid_ids; + process_ids = other.process_ids; + value_tables = other.value_tables; + model_groups = other.model_groups; + process_groups = other.process_groups; + + hardwired = other.hardwired; + max_particle_processes = other.max_particle_processes; + + scaling_min_range = other.scaling_min_range; + scaling_fraction = other.scaling_fraction; + + return *this; } }; @@ -186,17 +239,52 @@ struct PhysicsTrackInitializer * greatest number of element components of any material in the problem. This * can be used for the physics to calculate microscopic cross sections. */ -struct PhysicsStatePointers +template +struct PhysicsStateData { - Span state; //!< Track state [track] - Span per_process_xs; //!< XS [track][particle process] + template + using StateData = celeritas::StatePie; + template + using Data = celeritas::Pie; + + StateData state; //!< Track state [track] + Data per_process_xs; //!< XS [track][particle process] //! True if assigned explicit CELER_FUNCTION operator bool() const { return !state.empty(); } //! State size CELER_FUNCTION size_type size() const { return state.size(); } + + //! Assign from another set of states + template + PhysicsStateData& operator=(PhysicsStateData& other) + { + CELER_EXPECT(other); + state = other.state; + per_process_xs = other.per_process_xs; + return *this; + } }; +#ifndef __CUDA_ARCH__ +//---------------------------------------------------------------------------// +/*! + * Resize a material state in host code. + */ +template +inline void +resize(PhysicsStateData* data, + const PhysicsParamsData& params, + size_type size) +{ + CELER_EXPECT(size > 0); + CELER_EXPECT(params.max_particle_processes > 0); + make_pie_builder(&data->state).resize(size); + make_pie_builder(&data->per_process_xs) + .resize(size * params.max_particle_processes); +} +#endif + //---------------------------------------------------------------------------// } // namespace celeritas diff --git a/src/physics/base/PhysicsParams.cc b/src/physics/base/PhysicsParams.cc index ea81c63b54..33c81aa8f2 100644 --- a/src/physics/base/PhysicsParams.cc +++ b/src/physics/base/PhysicsParams.cc @@ -8,13 +8,36 @@ #include "PhysicsParams.hh" #include +#include +#include +#include "base/Assert.hh" #include "base/Range.hh" #include "base/VectorUtils.hh" +#include "comm/Logger.hh" #include "ParticleParams.hh" +#include "physics/grid/ValueGridInserter.hh" #include "physics/material/MaterialParams.hh" namespace celeritas { +namespace +{ +const char* to_cstring(ValueGridType grid) +{ + switch (grid) + { + case ValueGridType::macro_xs: + return "macro_xs"; + case ValueGridType::energy_loss: + return "energy_loss"; + case ValueGridType::range: + return "range"; + default: + return "[INVALID]"; + } +} +} // namespace + //---------------------------------------------------------------------------// /*! * Construct with processes and helper classes. @@ -28,8 +51,46 @@ PhysicsParams::PhysicsParams(Input inp) : processes_(std::move(inp.processes)) CELER_EXPECT(inp.particles); CELER_EXPECT(inp.materials); - // Resize: processes for each particle - processes_by_particle_.resize(inp.particles->size()); + // Emit models for associated proceses + models_ = this->build_models(); + + // Construct data + HostValue host_data; + this->build_options(inp.options, &host_data); + this->build_ids(*inp.particles, &host_data); + this->build_xs(*inp.materials, &host_data); + + CELER_LOG(debug) + << "Constructed physics sizes:" + << "\n reals: " << host_data.reals.size() + << "\n model_ids: " << host_data.model_ids.size() + << "\n value_grids: " << host_data.value_grids.size() + << "\n value_grid_ids: " << host_data.value_grid_ids.size() + << "\n process_ids: " << host_data.process_ids.size() + << "\n value_tables: " << host_data.value_tables.size() + << "\n model_groups: " << host_data.model_groups.size() + << "\n process_groups: " << host_data.process_groups.size(); + + data_ = PieMirror{std::move(host_data)}; +} + +//---------------------------------------------------------------------------// +/*! + * Get the list of process IDs that apply to a particle type. + */ +auto PhysicsParams::processes(ParticleId id) const -> SpanConstProcessId +{ + CELER_EXPECT(id < this->num_processes()); + const auto& data = this->host_pointers(); + return data.process_ids[data.process_groups[id].processes]; +} + +//---------------------------------------------------------------------------// +// HELPER FUNCTIONS +//---------------------------------------------------------------------------// +auto PhysicsParams::build_models() const -> VecModel +{ + VecModel models; // Construct models, assigning each model ID ModelIdGenerator next_model_id; @@ -37,23 +98,261 @@ PhysicsParams::PhysicsParams(Input inp) : processes_(std::move(inp.processes)) { auto new_models = processes_[process_idx]->build_models(next_model_id); CELER_ASSERT(!new_models.empty()); - for (const SPConstModel& model : new_models) + for (SPConstModel& model : new_models) { + CELER_ASSERT(model); ModelId model_id = next_model_id(); CELER_ASSERT(model->model_id() == model_id); - // TODO: Add processes to particle-per-process + // Save model and the process that it belongs to + models.push_back({std::move(model), ProcessId{process_idx}}); + } + } + + CELER_ENSURE(!models.empty()); + return models; +} + +//---------------------------------------------------------------------------// +/*! + * Construct on-device physics options. + */ +void PhysicsParams::build_options(const Options& opts, HostValue* data) const +{ + CELER_VALIDATE( + opts.max_step_over_range > 0, + "Non-positive max_step_over_range=" << opts.max_step_over_range); + CELER_VALIDATE(opts.min_range > 0, + "Non-positive min_range=" << opts.min_range); + data->scaling_min_range = opts.min_range; + data->scaling_fraction = opts.max_step_over_range; +} + +//---------------------------------------------------------------------------// +/*! + * Construct particle -> process -> model mappings. + */ +void PhysicsParams::build_ids(const ParticleParams& particles, + HostValue* data) const +{ + CELER_EXPECT(data); + CELER_EXPECT(!models_.empty()); + using ModelRange = std::tuple; + + // Note: use map to keep ProcessId sorted + std::vector>> particle_models( + particles.size()); + + // Construct particle -> process -> model map + for (auto model_idx : range(this->num_models())) + { + const Model& m = *models_[model_idx].first; + const ProcessId process_id = models_[model_idx].second; + for (const Applicability& applic : m.applicability()) + { + if (applic.material) + { + CELER_NOT_IMPLEMENTED("Material-dependent models"); + } + CELER_VALIDATE(applic.particle < particles.size(), + "Invalid particle ID"); + CELER_ASSERT(applic.lower < applic.upper); + particle_models[applic.particle.get()][process_id].push_back( + {applic.lower.value(), + applic.upper.value(), + ModelId{model_idx}}); + } + } + + auto process_groups = make_pie_builder(&data->process_groups); + auto process_ids = make_pie_builder(&data->process_ids); + auto model_groups = make_pie_builder(&data->model_groups); + auto model_ids = make_pie_builder(&data->model_ids); + auto reals = make_pie_builder(&data->reals); + + process_groups.reserve(particle_models.size()); + + // Loop over particle IDs, set ProcessGroup + for (auto particle_idx : range(particles.size())) + { + auto& process_to_models = particle_models[particle_idx]; + if (process_to_models.empty()) + { + CELER_LOG(warning) + << "No processes are defined for particle '" + << particles.id_to_label(ParticleId{particle_idx}); } + data->max_particle_processes = std::max( + data->max_particle_processes, process_to_models.size()); - // TODO: Loop over materials and applicable particles, and call - // step_limits with each particle type and material ID to get physics - // tables + std::vector temp_processes; + std::vector temp_model_groups; + temp_processes.reserve(process_to_models.size()); + temp_model_groups.reserve(process_to_models.size()); + for (auto& pid_models : process_to_models) + { + // Add process ID + temp_processes.push_back(pid_models.first); - // Move models to the end of the vector - celeritas::move_extend(std::move(new_models), &models_); + std::vector& models = pid_models.second; + CELER_ASSERT(!models.empty()); + + // Construct model group + std::vector temp_energy_grid; + std::vector temp_models; + temp_energy_grid.reserve(models.size() + 1); + temp_models.reserve(models.size()); + + // Sort, and add the first grid point + std::sort(models.begin(), models.end()); + temp_energy_grid.push_back(std::get<0>(models[0])); + + for (const ModelRange& r : models) + { + CELER_VALIDATE( + temp_energy_grid.back() == std::get<0>(r), + "Models for process '" + << this->process(pid_models.first).label() + << "' of particle type '" + << particles.id_to_label(ParticleId{particle_idx}) + << "' are discontinuous in energy"); + temp_energy_grid.push_back(std::get<1>(r)); + temp_models.push_back(std::get<2>(r)); + } + + ModelGroup mgroup; + mgroup.energy = reals.insert_back(temp_energy_grid.begin(), + temp_energy_grid.end()); + mgroup.model = model_ids.insert_back(temp_models.begin(), + temp_models.end()); + CELER_ASSERT(mgroup); + temp_model_groups.push_back(mgroup); + } + + ProcessGroup pgroup; + pgroup.processes = process_ids.insert_back(temp_processes.begin(), + temp_processes.end()); + pgroup.models = model_groups.insert_back(temp_model_groups.begin(), + temp_model_groups.end()); + // NOTE: data tables will be assigned later + CELER_ASSERT(pgroup); + process_groups.push_back(pgroup); } - CELER_NOT_IMPLEMENTED("constructing physics params"); + // TODO: hardwired models + + CELER_ENSURE(*data); +} + +//---------------------------------------------------------------------------// +/*! + * Construct cross section data. + */ +void PhysicsParams::build_xs(const MaterialParams& mats, HostValue* data) const +{ + CELER_EXPECT(*data); + + using UPGridBuilder = Process::UPConstGridBuilder; + + ValueGridInserter insert_grid(&data->reals, &data->value_grids); + auto value_tables = make_pie_builder(&data->value_tables); + auto value_grid_ids = make_pie_builder(&data->value_grid_ids); + auto build_grid + = [insert_grid](const UPGridBuilder& builder) -> ValueGridId { + return builder ? builder->build(insert_grid) : ValueGridId{}; + }; + + Applicability applic; + for (auto particle_idx : range(data->process_groups.size())) + { + applic.particle = ParticleId{particle_idx}; + + // Processes for this particle + ProcessGroup& process_group + = data->process_groups[ParticleId{particle_idx}]; + Span processes + = data->process_ids[process_group.processes]; + Span model_groups + = data->model_groups[process_group.models]; + CELER_ASSERT(processes.size() == model_groups.size()); + + // Material-dependent physics tables, one per particle-process + ValueGridArray> temp_tables; + for (auto& vec : temp_tables) + { + vec.resize(processes.size()); + } + + // Loop over per-particle processes + for (auto pp_idx : range(processes.size())) + { + // Get energy bounds for this process + Span energy_grid + = data->reals[model_groups[pp_idx].energy]; + applic.lower = Applicability::Energy{energy_grid.front()}; + applic.upper = Applicability::Energy{energy_grid.back()}; + CELER_ASSERT(applic.lower < applic.upper); + + const Process& proc = this->process(processes[pp_idx]); + + // Grid IDs for each grid type, each material + ValueGridArray> temp_grid_ids; + for (auto& vec : temp_grid_ids) + { + vec.resize(mats.size()); + } + + // Loop over materials + for (auto mat_idx : range(mats.size())) + { + applic.material = MaterialId{mat_idx}; + + // Construct step limit builders + auto builders = proc.step_limits(applic); + CELER_VALIDATE( + std::any_of(builders.begin(), + builders.end(), + [](const UPGridBuilder& p) { return bool(p); }), + "Process '" << proc.label() + << "' has neither interaction nor energy " + "loss"); + + // Construct grids + for (auto vgt : range(size_type(ValueGridType::size_))) + { + temp_grid_ids[vgt][mat_idx] = build_grid(builders[vgt]); + } + } + + // Outer loop over grid types + for (auto vgt : range(size_type(ValueGridType::size_))) + { + if (!std::any_of(temp_grid_ids[vgt].begin(), + temp_grid_ids[vgt].end(), + [](ValueGridId id) { return bool(id); })) + { + // Skip this table type since it's not present for any + // material for this particle process + CELER_LOG(debug) << "No " << to_cstring(ValueGridType(vgt)) + << " for process " << proc.label(); + continue; + } + + // Construct value grid table + ValueTable& temp_table = temp_tables[vgt][pp_idx]; + temp_table.material = value_grid_ids.insert_back( + temp_grid_ids[vgt].begin(), temp_grid_ids[vgt].end()); + CELER_ASSERT(temp_table.material.size() == mats.size()); + } + } + + // Construct value tables + for (auto vgt : range(size_type(ValueGridType::size_))) + { + process_group.tables[vgt] = value_tables.insert_back( + temp_tables[vgt].begin(), temp_tables[vgt].end()); + } + } } //---------------------------------------------------------------------------// diff --git a/src/physics/base/PhysicsParams.hh b/src/physics/base/PhysicsParams.hh index fac41cecce..f52b25c676 100644 --- a/src/physics/base/PhysicsParams.hh +++ b/src/physics/base/PhysicsParams.hh @@ -9,8 +9,9 @@ #include #include -#include "base/DeviceVector.hh" +#include "base/PieMirror.hh" #include "base/Types.hh" +#include "base/Units.hh" #include "Model.hh" #include "Process.hh" #include "PhysicsInterface.hh" @@ -32,6 +33,13 @@ class ParticleParams; * * During construction it constructs models and their corresponding list of * \c ModelId values, as well as the tables of cross section data. + * + * Input options are: + * - \c min_range: below this value, there is no extra transformation from + * particle range to step length. + * - \c max_step_over_range: at higher energies (longer ranges), gradually + * decrease the minimum step length until it's this fraction of the tabulated + * range. */ class PhysicsParams { @@ -43,13 +51,17 @@ class PhysicsParams using SPConstProcess = std::shared_ptr; using VecProcess = std::vector; using SpanConstProcessId = Span; + using HostRef + = PhysicsParamsData; + using DeviceRef + = PhysicsParamsData; //!@} //! Global physics configuration options struct Options { - real_type max_step_over_range; //!< alpha_r, limit step range at high E - real_type min_step; //!< rho_R, minimum range at low E + real_type min_range = 1 * units::millimeter; //!< rho_R + real_type max_step_over_range = 0.2; //!< alpha_r }; //! Physics parameter construction arguments @@ -69,16 +81,16 @@ class PhysicsParams //// HOST ACCESSORS //// //! Number of models - size_type num_models() const { return models_.size(); } + ModelId::value_type num_models() const { return models_.size(); } //! Number of processes - size_type num_processes() const { return processes_.size(); } + ProcessId::value_type num_processes() const { return processes_.size(); } - //! Number of particle types - size_type num_particles() const { return processes_by_particle_.size(); } + // Number of particle types + inline ParticleId::value_type num_particles() const; - //! Maximum number of processes that apply to any one particle - size_type max_particle_processes() const { return max_pbp_; } + // Maximum number of processes that apply to any one particle + inline ProcessId::value_type max_particle_processes() const; // Get a model inline const Model& model(ModelId) const; @@ -86,67 +98,72 @@ class PhysicsParams // Get a process inline const Process& process(ProcessId) const; - // Get the processes that apply to a particlular particle - inline SpanConstProcessId processes(ParticleId) const; + // Get the processes that apply to a particular particle + SpanConstProcessId processes(ParticleId) const; - //// DEVICE ACCESSORS //// + //! Access material properties on the host + const HostRef& host_pointers() const { return data_.host(); } - // Get managed data - PhysicsParamsPointers device_pointers() const; + //! Access material properties on the device + const DeviceRef& device_pointers() const { return data_.device(); } private: using SPConstModel = std::shared_ptr; + using VecModel = std::vector>; + using HostValue = PhysicsParamsData; - //// HOST DATA //// - - VecProcess processes_; - std::vector> processes_by_particle_; - std::vector models_; - size_type max_pbp_ = 0; + // Host metadata/access + VecProcess processes_; + VecModel models_; - //// DEVICE DATA //// + // Host/device storage and reference + PieMirror data_; - DeviceVector grid_inputs_; - DeviceVector xsgrid_data_; - DeviceVector value_grids_; - DeviceVector value_tables_; - DeviceVector process_ids_; - DeviceVector model_ids_; - DeviceVector model_energy_; - DeviceVector model_groups_; - DeviceVector process_groups_; + private: + VecModel build_models() const; + void build_options(const Options& opts, HostValue* data) const; + void build_ids(const ParticleParams& particles, HostValue* data) const; + void build_xs(const MaterialParams& mats, HostValue* data) const; }; //---------------------------------------------------------------------------// // INLINE FUNCTIONS //---------------------------------------------------------------------------// /*! - * Get a model. + * Number of particle types. */ -const Model& PhysicsParams::model(ModelId id) const +auto PhysicsParams::num_particles() const -> ParticleId::value_type { - CELER_EXPECT(id < this->num_models()); - return *models_[id.get()]; + return this->host_pointers().process_ids.size(); } //---------------------------------------------------------------------------// /*! - * Get a process. + * Number of particle types. */ -const Process& PhysicsParams::process(ProcessId id) const +auto PhysicsParams::max_particle_processes() const -> ProcessId::value_type { - CELER_EXPECT(id < this->num_processes()); - return *processes_[id.get()]; + return this->host_pointers().max_particle_processes; +} + +//---------------------------------------------------------------------------// +/*! + * Get a model. + */ +const Model& PhysicsParams::model(ModelId id) const +{ + CELER_EXPECT(id < this->num_models()); + return *models_[id.get()].first; } //---------------------------------------------------------------------------// /*! - * Get the list of process IDs that apply to a particle type. + * Get a process. */ -auto PhysicsParams::processes(ParticleId id) const -> SpanConstProcessId +const Process& PhysicsParams::process(ProcessId id) const { CELER_EXPECT(id < this->num_processes()); - return make_span(processes_by_particle_[id.get()]); + return *processes_[id.get()]; } //---------------------------------------------------------------------------// diff --git a/src/physics/base/PhysicsStepUtils.i.hh b/src/physics/base/PhysicsStepUtils.i.hh index a58a34c63b..868148340b 100644 --- a/src/physics/base/PhysicsStepUtils.i.hh +++ b/src/physics/base/PhysicsStepUtils.i.hh @@ -10,7 +10,6 @@ #include "base/Algorithms.hh" #include "base/NumericLimits.hh" #include "base/Range.hh" -#include "physics/grid/PhysicsGridCalculator.hh" #include "Types.hh" namespace celeritas @@ -25,7 +24,7 @@ CELER_FUNCTION real_type calc_tabulated_physics_step( CELER_EXPECT(physics.has_interaction_mfp()); constexpr real_type inf = numeric_limits::infinity(); - using PTT = PhysicsTableType; + using VGT = ValueGridType; // Loop over all processes that apply to this track (based on particle // type) and calculate cross section and particle range. @@ -36,12 +35,12 @@ CELER_FUNCTION real_type calc_tabulated_physics_step( { const ParticleProcessId ppid{pp_idx}; real_type process_xs = 0; - if (const auto* table = physics.table(PTT::macro_xs, ppid)) + if (auto grid_id = physics.value_grid(VGT::macro_xs, ppid)) { // Calculate macroscopic cross section for this process, then // accumulate it into the total cross section and save the cross // section for later. - PhysicsGridCalculator calc_xs(*table, {}); + auto calc_xs = physics.make_calculator(grid_id); process_xs = calc_xs(particle.energy()); total_macro_xs += process_xs; } @@ -61,9 +60,9 @@ CELER_FUNCTION real_type calc_tabulated_physics_step( } physics.per_process_xs(ppid) = process_xs; - if (const auto* table = physics.table(PTT::range, ppid)) + if (auto grid_id = physics.value_grid(VGT::range, ppid)) { - PhysicsGridCalculator calc_range(*table, {}); + auto calc_range = physics.make_calculator(grid_id); real_type process_range = calc_range(particle.energy()); // TODO: scale range by sqrt(particle.energy() / minKE) // if < minKE?? @@ -92,7 +91,7 @@ CELER_FUNCTION real_type calc_energy_loss(const ParticleTrackView& particle, { CELER_EXPECT(step >= 0); - using PTT = PhysicsTableType; + using VGT = ValueGridType; // Loop over all processes that apply to this track (based on particle // type) and calculate cross section and particle range. @@ -101,9 +100,9 @@ CELER_FUNCTION real_type calc_energy_loss(const ParticleTrackView& particle, range(physics.num_particle_processes())) { const ParticleProcessId ppid{pp_idx}; - if (const auto* table = physics.table(PTT::energy_loss, ppid)) + if (auto grid_id = physics.value_grid(VGT::energy_loss, ppid)) { - PhysicsGridCalculator calc_eloss_rate(*table, {}); + auto calc_eloss_rate = physics.make_calculator(grid_id); total_eloss_rate += calc_eloss_rate(particle.energy()); } } diff --git a/src/physics/base/PhysicsTrackView.hh b/src/physics/base/PhysicsTrackView.hh index d3dcd91c65..b2c365a3c9 100644 --- a/src/physics/base/PhysicsTrackView.hh +++ b/src/physics/base/PhysicsTrackView.hh @@ -10,6 +10,9 @@ #include "PhysicsInterface.hh" #include "base/Macros.hh" #include "base/Types.hh" +#include "physics/base/Units.hh" +#include "physics/grid/GridIdFinder.hh" +#include "physics/grid/PhysicsGridCalculator.hh" #include "physics/material/Types.hh" #include "Types.hh" @@ -27,9 +30,13 @@ class PhysicsTrackView public: //!@{ //! Type aliases - using SpanConstProcessId = Span; - using PhysicsGridPointers = ValueGrid; - using Initializer_t = PhysicsTrackInitializer; + using Initializer_t = PhysicsTrackInitializer; + using PhysicsParamsPointers + = PhysicsParamsData; + using PhysicsStatePointers + = PhysicsStateData; + + using ModelFinder = GridIdFinder; //!@} public: @@ -82,17 +89,22 @@ class PhysicsTrackView inline CELER_FUNCTION ProcessId process(ParticleProcessId) const; // Get table, null if not present for this particle/material/type - inline CELER_FUNCTION const PhysicsGridPointers* - table(PhysicsTableType table, ParticleProcessId) const; + inline CELER_FUNCTION ValueGridId value_grid(ValueGridType table, + ParticleProcessId) const; // Models that apply to the given process ID - inline CELER_FUNCTION const ModelGroup& models(ParticleProcessId) const; + inline CELER_FUNCTION + ModelFinder make_model_finder(ParticleProcessId) const; //// STATIC FUNCTIONS (depend only on params data) //// // Calculate scaled step range inline CELER_FUNCTION real_type range_to_step(real_type range) const; + // Construct a grid calculator from a physics table + inline CELER_FUNCTION + PhysicsGridCalculator make_calculator(ValueGridId) const; + //// SCRATCH SPACE //// // Access scratch space for particle-process cross section calculations diff --git a/src/physics/base/PhysicsTrackView.i.hh b/src/physics/base/PhysicsTrackView.i.hh index fe2dda8b10..00036446cc 100644 --- a/src/physics/base/PhysicsTrackView.i.hh +++ b/src/physics/base/PhysicsTrackView.i.hh @@ -162,7 +162,7 @@ CELER_FUNCTION ParticleProcessId::value_type CELER_FUNCTION ProcessId PhysicsTrackView::process(ParticleProcessId ppid) const { CELER_EXPECT(ppid < this->num_particle_processes()); - return this->process_group().processes[ppid.get()]; + return params_.process_ids[this->process_group().processes[ppid.get()]]; } //---------------------------------------------------------------------------// @@ -176,42 +176,54 @@ CELER_FUNCTION ProcessId PhysicsTrackView::process(ParticleProcessId ppid) const * associated value (e.g. if the table type is "energy_loss" and the process is * not a slowing-down process). */ -CELER_FUNCTION auto PhysicsTrackView::table(PhysicsTableType table_type, - ParticleProcessId ppid) const - -> const PhysicsGridPointers* +CELER_FUNCTION auto PhysicsTrackView::value_grid(ValueGridType table_type, + ParticleProcessId ppid) const + -> ValueGridId { - CELER_EXPECT(int(table_type) < int(PhysicsTableType::size_)); + CELER_EXPECT(int(table_type) < int(ValueGridType::size_)); CELER_EXPECT(ppid < this->num_particle_processes()); - const ValueTable& table + ValueTableId table_id = this->process_group().tables[int(table_type)][ppid.get()]; + CELER_ASSERT(table_id); + const ValueTable& table = params_.value_tables[table_id]; if (!table) - return nullptr; + return {}; // No table for this process CELER_EXPECT(material_ < table.material.size()); - CELER_ENSURE(table.material[material_.get()]); - return &table.material[material_.get()]; + auto grid_id_ref = table.material[material_.get()]; + if (!grid_id_ref) + return {}; // No table for this particular material + + return params_.value_grid_ids[grid_id_ref]; } //---------------------------------------------------------------------------// /*! * Models that apply to the given process ID. */ -CELER_FUNCTION const ModelGroup& - PhysicsTrackView::models(ParticleProcessId ppid) const +CELER_FUNCTION auto +PhysicsTrackView::make_model_finder(ParticleProcessId ppid) const + -> ModelFinder { CELER_EXPECT(ppid < this->num_particle_processes()); - return this->process_group().models[ppid.get()]; + const ModelGroup& mg + = params_.model_groups[this->process_group().models[ppid.get()]]; + return ModelFinder(params_.reals[mg.energy], params_.model_ids[mg.model]); } //---------------------------------------------------------------------------// /*! * Calculate scaled step range. * - * This is the updated step function given by Eq. 7.4: \f[ + * This is the updated step function given by Eq. 7.4 of Geant4 Physics + * Reference Manual, Release 10.6: \f[ s = \alpha r + \rho (1 - \alpha) (2 - \frac{\rho}{r}) \f] - * where alpha is \c max_step_over_range and rho is \c min_step + * where alpha is \c scaling_fraction and rho is \c scaling_min_range . + * + * Below scaling_min_range, no step scaling is applied, but the step can still + * be arbitrarily small. */ CELER_FUNCTION real_type PhysicsTrackView::range_to_step(real_type range) const { @@ -227,6 +239,17 @@ CELER_FUNCTION real_type PhysicsTrackView::range_to_step(real_type range) const return range; } +//---------------------------------------------------------------------------// +/*! + * Access data for constructing PhysicsGridCalculator. + */ +CELER_FUNCTION PhysicsGridCalculator +PhysicsTrackView::make_calculator(ValueGridId id) const +{ + CELER_EXPECT(id < params_.value_grids.size()); + return PhysicsGridCalculator(params_.value_grids[id], params_.reals); +} + //---------------------------------------------------------------------------// /*! * Access scratch space for particle-process cross section calculations. @@ -239,7 +262,7 @@ CELER_FUNCTION real_type& CELER_EXPECT(ppid < this->num_particle_processes()); auto idx = thread_.get() * params_.max_particle_processes + ppid.get(); CELER_ENSURE(idx < states_.per_process_xs.size()); - return states_.per_process_xs[idx]; + return states_.per_process_xs[PieId(idx)]; } //---------------------------------------------------------------------------// @@ -252,7 +275,7 @@ real_type PhysicsTrackView::per_process_xs(ParticleProcessId ppid) const CELER_EXPECT(ppid < this->num_particle_processes()); auto idx = thread_.get() * params_.max_particle_processes + ppid.get(); CELER_ENSURE(idx < states_.per_process_xs.size()); - return states_.per_process_xs[idx]; + return states_.per_process_xs[PieId(idx)]; } //---------------------------------------------------------------------------// @@ -279,20 +302,20 @@ CELER_FUNCTION ProcessId PhysicsTrackView::eplusgg_process_id() const //! Get the thread-local state (mutable) CELER_FUNCTION PhysicsTrackState& PhysicsTrackView::state() { - return states_.state[thread_.get()]; + return states_.state[thread_]; } //! Get the thread-local state (const) CELER_FUNCTION const PhysicsTrackState& PhysicsTrackView::state() const { - return states_.state[thread_.get()]; + return states_.state[thread_]; } //! Get the group of processes that apply to the particle CELER_FUNCTION const ProcessGroup& PhysicsTrackView::process_group() const { - CELER_EXPECT(particle_ < params_.particle.size()); - return params_.particle[particle_.get()]; + CELER_EXPECT(particle_ < params_.process_groups.size()); + return params_.process_groups[particle_]; } //---------------------------------------------------------------------------// diff --git a/src/physics/base/Types.hh b/src/physics/base/Types.hh index 6961905c8f..18e8f08204 100644 --- a/src/physics/base/Types.hh +++ b/src/physics/base/Types.hh @@ -23,7 +23,7 @@ using ModelId = OpaqueId; using ProcessId = OpaqueId; //! Opaque index of a process applicable to a single particle type -using ParticleProcessId = OpaqueId; +using ParticleProcessId = OpaqueId; //! Opaque index of electron subshell using SubshellId = OpaqueId; diff --git a/src/physics/grid/ValueGridBuilder.cc b/src/physics/grid/ValueGridBuilder.cc index b741012ec1..4ae7c3e9e0 100644 --- a/src/physics/grid/ValueGridBuilder.cc +++ b/src/physics/grid/ValueGridBuilder.cc @@ -135,7 +135,7 @@ ValueGridXsBuilder::ValueGridXsBuilder(real_type emin, /*! * Construct on device. */ -void ValueGridXsBuilder::build(ValueGridInserter insert) const +auto ValueGridXsBuilder::build(ValueGridInserter insert) const -> ValueGridId { auto log_energy = UniformGridData::from_bounds(log_emin_, log_emax_, xs_.size()); @@ -145,9 +145,10 @@ void ValueGridXsBuilder::build(ValueGridInserter insert) const auto prime_index = grid.find(log_eprime_); CELER_ASSERT(soft_equal(grid[prime_index], log_eprime_)); - insert(UniformGridData::from_bounds(log_emin_, log_emax_, xs_.size()), - prime_index, - make_span(xs_)); + return insert( + UniformGridData::from_bounds(log_emin_, log_emax_, xs_.size()), + prime_index, + make_span(xs_)); } //---------------------------------------------------------------------------// @@ -170,10 +171,11 @@ ValueGridLogBuilder::ValueGridLogBuilder(real_type emin, /*! * Construct on device. */ -void ValueGridLogBuilder::build(ValueGridInserter insert) const +auto ValueGridLogBuilder::build(ValueGridInserter insert) const -> ValueGridId { - insert(UniformGridData::from_bounds(log_emin_, log_emax_, xs_.size()), - make_span(xs_)); + return insert( + UniformGridData::from_bounds(log_emin_, log_emax_, xs_.size()), + make_span(xs_)); } //---------------------------------------------------------------------------// @@ -210,10 +212,13 @@ ValueGridGenericBuilder::ValueGridGenericBuilder(VecReal grid, VecReal value) /*! * Construct grid data in the given mutable insert. */ -void ValueGridGenericBuilder::build(ValueGridInserter insert) const +auto ValueGridGenericBuilder::build(ValueGridInserter insert) const + -> ValueGridId { insert({make_span(grid_), grid_interp_}, {make_span(value_), value_interp_}); + CELER_NOT_IMPLEMENTED("generic grids"); + return {}; } //---------------------------------------------------------------------------// diff --git a/src/physics/grid/ValueGridBuilder.hh b/src/physics/grid/ValueGridBuilder.hh index e25aa834fc..6fa89bb200 100644 --- a/src/physics/grid/ValueGridBuilder.hh +++ b/src/physics/grid/ValueGridBuilder.hh @@ -9,6 +9,7 @@ #include #include +#include "base/Pie.hh" #include "base/Span.hh" #include "base/Types.hh" @@ -24,12 +25,18 @@ class ValueGridInserter; */ class ValueGridBuilder { + public: + //!@{ + //! Type aliases + using ValueGridId = PieId; + //!@} + public: //! Virtual destructor for polymorphic deletion virtual ~ValueGridBuilder() = 0; //! Construct the grid given a mutable reference to a store - virtual void build(ValueGridInserter) const = 0; + virtual ValueGridId build(ValueGridInserter) const = 0; }; //---------------------------------------------------------------------------// @@ -62,7 +69,7 @@ class ValueGridXsBuilder final : public ValueGridBuilder VecReal xs); // Construct in the given store - void build(ValueGridInserter) const final; + ValueGridId build(ValueGridInserter) const final; private: real_type log_emin_; @@ -83,6 +90,7 @@ class ValueGridLogBuilder final : public ValueGridBuilder //!@{ //! Type aliases using VecReal = std::vector; + using Id = PieId; //!@} public: @@ -90,7 +98,7 @@ class ValueGridLogBuilder final : public ValueGridBuilder ValueGridLogBuilder(real_type emin, real_type emax, VecReal value); // Construct in the given store - void build(ValueGridInserter) const final; + ValueGridId build(ValueGridInserter) const final; private: real_type log_emin_; @@ -108,6 +116,7 @@ class ValueGridGenericBuilder final : public ValueGridBuilder //!@{ //! Type aliases using VecReal = std::vector; + using Id = PieId; //!@} public: @@ -121,7 +130,7 @@ class ValueGridGenericBuilder final : public ValueGridBuilder ValueGridGenericBuilder(VecReal grid, VecReal value); // Construct in the given store - void build(ValueGridInserter) const final; + ValueGridId build(ValueGridInserter) const final; private: VecReal grid_; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c4bb7d079f..c46a47dde5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -165,13 +165,15 @@ add_library(CeleritasPhysicsTest physics/SecondaryIO.cc physics/base/MockModel.cc physics/base/MockProcess.cc + physics/base/PhysicsTestBase.cc ) target_link_libraries(CeleritasPhysicsTest PRIVATE celeritas CeleritasTest) set(CELERITASTEST_LINK_LIBRARIES CeleritasPhysicsTest) celeritas_setup_tests(SERIAL PREFIX physics/base) celeritas_cudaoptional_test(physics/base/Particle) -celeritas_add_test(physics/base/Physics.test.cc ${_not_impl}) +celeritas_add_test(physics/base/Physics.test.cc) +celeritas_add_test(physics/base/PhysicsStepUtils.test.cc ${_not_impl}) celeritas_setup_tests(SERIAL PREFIX physics/grid) celeritas_add_test(physics/grid/GridIdFinder.test.cc) diff --git a/test/physics/base/Physics.test.cc b/test/physics/base/Physics.test.cc index d35910bc95..4ac0d9f91c 100644 --- a/test/physics/base/Physics.test.cc +++ b/test/physics/base/Physics.test.cc @@ -5,130 +5,327 @@ //---------------------------------------------------------------------------// //! \file Physics.test.cc //---------------------------------------------------------------------------// -#include "physics/base/PhysicsTrackView.hh" -#include "physics/base/PhysicsStepUtils.hh" #include "physics/base/PhysicsParams.hh" +#include "physics/base/PhysicsTrackView.hh" #include "celeritas_test.hh" -#include "MockProcess.hh" -#include "physics/material/MaterialParams.hh" +#include "base/Range.hh" +#include "base/PieStateStore.hh" #include "physics/base/ParticleParams.hh" -// #include "Physics.test.hh" +#include "PhysicsTestBase.hh" using namespace celeritas; using namespace celeritas_test; -using celeritas::units::MevEnergy; //---------------------------------------------------------------------------// -// TEST HARNESS BASE +// PHYSICS BASE CLASS +//---------------------------------------------------------------------------// + +class PhysicsParamsTest : public PhysicsTestBase +{ +}; + +TEST_F(PhysicsParamsTest, accessors) +{ + const PhysicsParams& p = *this->physics(); + + EXPECT_EQ(5, p.num_processes()); + EXPECT_EQ(2 + 1 + 3 + 2 + 2, p.num_models()); + EXPECT_EQ(3, p.max_particle_processes()); + + // Test process names after construction + std::vector process_names; + for (auto process_idx : range(p.num_processes())) + { + process_names.push_back(p.process(ProcessId{process_idx}).label()); + } + const std::string expected_process_names[] + = {"scattering", "absorption", "purrs", "hisses", "meows"}; + EXPECT_VEC_EQ(expected_process_names, process_names); + + // Test model names after construction + std::vector model_names; + for (auto model_idx : range(p.num_models())) + { + model_names.push_back(p.model(ModelId{model_idx}).label()); + } + const std::string expected_model_names[] + = {"MockModel(0, p=0, emin=1e-06, emax=100)", + "MockModel(1, p=1, emin=1, emax=100)", + "MockModel(2, p=0, emin=1e-06, emax=100)", + "MockModel(3, p=1, emin=0.001, emax=1)", + "MockModel(4, p=1, emin=1, emax=10)", + "MockModel(5, p=1, emin=10, emax=100)", + "MockModel(6, p=2, emin=0.001, emax=1)", + "MockModel(7, p=2, emin=1, emax=100)", + "MockModel(8, p=1, emin=0.001, emax=10)", + "MockModel(9, p=2, emin=0.001, emax=10)"}; + EXPECT_VEC_EQ(expected_model_names, model_names); + + // Test host-accessible process map + std::vector process_map; + for (auto particle_idx : range(this->particles()->size())) + { + std::string prefix + = this->particles()->id_to_label(ParticleId{particle_idx}); + prefix.push_back(':'); + for (ProcessId process_id : p.processes(ParticleId{particle_idx})) + { + process_map.push_back(prefix + process_names[process_id.get()]); + } + } + const std::string expected_process_map[] = {"gamma:scattering", + "gamma:absorption", + "celeriton:scattering", + "celeriton:purrs", + "celeriton:meows", + "anti-celeriton:hisses", + "anti-celeriton:meows"}; + EXPECT_VEC_EQ(expected_process_map, process_map); +} + +//---------------------------------------------------------------------------// +// PHYSICS TRACK VIEW (HOST) //---------------------------------------------------------------------------// -class PhysicsTest : public celeritas::Test +class PhysicsTrackViewHostTest : public PhysicsTestBase { - protected: + using Base = PhysicsTestBase; + + public: + //!@{ + //! Type aliases + using StateStore = PieStateStore; + using ParamsHostRef + = PhysicsParamsData; + using MevEnergy = celeritas::units::MevEnergy; + //!@} + void SetUp() override { - using namespace celeritas::units; + Base::SetUp(); - // Create materials + // Make one state per particle + auto state_size = this->particles()->size(); + + CELER_ASSERT(this->physics()); + params_ref = this->physics()->host_pointers(); + state = StateStore(*this->physics(), state_size); + + // Save mapping of process label -> ID + for (auto process_idx : range(this->physics()->num_processes())) { - MaterialParams::Input inp; - inp.elements = {{1, AmuMass{1.0}, "celerogen"}, - {4, AmuMass{10.0}, "celerinium"}}; - inp.materials.push_back({1e20, - 300, - MatterState::gas, - {{ElementId{0}, 1.0}}, - "lo density celerogen"}); - inp.materials.push_back({1e21, - 300, - MatterState::liquid, - {{ElementId{0}, 1.0}}, - "hi density celerogen"}); - inp.materials.push_back({1e23, - 300, - MatterState::solid, - {{ElementId{1}, 1.0}}, - "solid celerinium"}); - materials = std::make_shared(std::move(inp)); + ProcessId id{process_idx}; + process_names[this->physics()->process(id).label()] = id; } - // Create particles + } + + PhysicsTrackView make_track_view(const char* particle, MaterialId mid) + { + CELER_EXPECT(particle && mid); + + auto pid = this->particles()->find(particle); + CELER_ASSERT(pid); + CELER_ASSERT(pid.get() < state.size()); + + ThreadId tid((pid.get() + 1) % state.size()); + + // Construct (thread depends on particle here to shake things up) and + // initialize + PhysicsTrackView phys(params_ref, state.ref(), pid, mid, tid); + phys = PhysicsTrackInitializer{}; + + return phys; + } + + ParticleProcessId + find_ppid(const PhysicsTrackView& track, const char* label) const + { + auto iter = process_names.find(label); + CELER_VALIDATE(iter != process_names.end(), + "No process named " << label); + ProcessId pid = iter->second; + for (auto pp_idx : range(track.num_particle_processes())) { - namespace pdg = celeritas::pdg; - - constexpr auto zero = celeritas::zero_quantity(); - constexpr auto stable - = celeritas::ParticleDef::stable_decay_constant(); - - ParticleParams::Input inp; - inp.push_back({"gamma", pdg::gamma(), zero, zero, stable}); - inp.push_back({"celeriton", - PDGNumber{1337}, - MevMass{1}, - ElementaryCharge{1}, - stable}); - inp.push_back({"anti-celeriton", - PDGNumber{-1337}, - MevMass{1}, - ElementaryCharge{-1}, - stable}); - particles = std::make_shared(std::move(inp)); + if (track.process(ParticleProcessId{pp_idx}) == pid) + return ParticleProcessId{pp_idx}; } + return {}; } - std::shared_ptr materials; - std::shared_ptr particles; - std::shared_ptr physics; + ParamsHostRef params_ref; + StateStore state; + std::map process_names; }; -TEST_F(PhysicsTest, accessors) +TEST_F(PhysicsTrackViewHostTest, accessors) { - // PTestInput input; - // input.num_threads = 0; - // auto result = p_test(input); - // PRINT_EXPECTED(result.foo); - // EXPECT_VEC_SOFT_EQ(expected_foo, result.foo); + PhysicsTrackView gamma = this->make_track_view("gamma", MaterialId{0}); + PhysicsTrackView celer = this->make_track_view("celeriton", MaterialId{1}); + const PhysicsTrackView& gamma_cref = gamma; + + // Interaction MFP + { + EXPECT_FALSE(gamma_cref.has_interaction_mfp()); + + gamma.interaction_mfp(1.234); + celer.interaction_mfp(2.345); + EXPECT_DOUBLE_EQ(1.234, gamma_cref.interaction_mfp()); + EXPECT_DOUBLE_EQ(2.345, celer.interaction_mfp()); + } + + // Cross sections + { + const PhysicsTrackView& gamma_cref = gamma; + + gamma.per_process_xs(ParticleProcessId{0}) = 1.2; + gamma.per_process_xs(ParticleProcessId{1}) = 10.0; + celer.per_process_xs(ParticleProcessId{0}) = 100.0; + EXPECT_DOUBLE_EQ(1.2, gamma_cref.per_process_xs(ParticleProcessId{0})); + EXPECT_DOUBLE_EQ(10.0, gamma_cref.per_process_xs(ParticleProcessId{1})); + EXPECT_DOUBLE_EQ(100.0, celer.per_process_xs(ParticleProcessId{0})); + } } -//---------------------------------------------------------------------------// -// HOST TESTS -//---------------------------------------------------------------------------// +TEST_F(PhysicsTrackViewHostTest, processes) +{ + // Gamma + { + const PhysicsTrackView phys + = this->make_track_view("gamma", MaterialId{0}); + + EXPECT_EQ(2, phys.num_particle_processes()); + const ParticleProcessId scat_ppid{0}; + const ParticleProcessId abs_ppid{1}; + EXPECT_EQ(ProcessId{0}, phys.process(scat_ppid)); + EXPECT_EQ(ProcessId{1}, phys.process(abs_ppid)); + } + + // Celeriton + { + const PhysicsTrackView phys + = this->make_track_view("celeriton", MaterialId{0}); + + EXPECT_EQ(3, phys.num_particle_processes()); + const ParticleProcessId scat_ppid{0}; + const ParticleProcessId purr_ppid{1}; + const ParticleProcessId meow_ppid{2}; + EXPECT_EQ(ProcessId{0}, phys.process(scat_ppid)); + EXPECT_EQ(ProcessId{2}, phys.process(purr_ppid)); + EXPECT_EQ(ProcessId{4}, phys.process(meow_ppid)); + } + + // Anti-celeriton + { + const PhysicsTrackView phys + = this->make_track_view("anti-celeriton", MaterialId{1}); + + EXPECT_EQ(2, phys.num_particle_processes()); + const ParticleProcessId hiss_ppid{0}; + const ParticleProcessId meow_ppid{1}; + EXPECT_EQ(ProcessId{3}, phys.process(hiss_ppid)); + EXPECT_EQ(ProcessId{4}, phys.process(meow_ppid)); + } +} -class PhysicsTestHost : public PhysicsTest +TEST_F(PhysicsTrackViewHostTest, value_grids) { - protected: - void SetUp() override { CELER_NOT_IMPLEMENTED("physics step tests"); } + std::vector grid_ids; - struct + for (const char* particle : {"gamma", "celeriton", "anti-celeriton"}) { - ParticleParamsData particle; - PhysicsParamsPointers physics; - } params; - struct + for (auto mat_idx : range(this->materials()->size())) + { + const PhysicsTrackView phys + = this->make_track_view(particle, MaterialId{mat_idx}); + + for (auto pp_idx : range(phys.num_particle_processes())) + { + for (ValueGridType vgt : range(ValueGridType::size_)) + { + auto id = phys.value_grid(vgt, ParticleProcessId{pp_idx}); + grid_ids.push_back(id ? id.get() : -1); + } + } + } + } + + // Grid IDs should be unique if they exist. Gammas should have fewer + // because there aren't any slowing down/range limiters. + const int expected_grid_ids[] + = {0, -1, -1, 3, -1, -1, 1, -1, -1, 4, -1, -1, 2, -1, -1, 5, + -1, -1, 6, -1, -1, 9, 10, 11, 18, 19, 20, 7, -1, -1, 12, 13, + 14, 21, 22, 23, 8, -1, -1, 15, 16, 17, 24, 25, 26, 27, 28, 29, + 36, 37, 38, 30, 31, 32, 39, 40, 41, 33, 34, 35, 42, 43, 44}; + EXPECT_VEC_EQ(expected_grid_ids, grid_ids); +} + +TEST_F(PhysicsTrackViewHostTest, calc_xs) +{ + // Cross sections: same across particle types, constant in energy, scale + // according to material number density + std::vector xs; + for (const char* particle : {"gamma", "celeriton"}) { - ParticleStateData particle; - PhysicsStatePointers physics; - } states; -}; + for (auto mat_idx : range(this->materials()->size())) + { + const PhysicsTrackView phys + = this->make_track_view(particle, MaterialId{mat_idx}); + auto scat_ppid = this->find_ppid(phys, "scattering"); + auto id = phys.value_grid(ValueGridType::macro_xs, scat_ppid); + ASSERT_TRUE(id); + auto calc_xs = phys.make_calculator(id); + xs.push_back(calc_xs(MevEnergy{1.0})); + } + } -TEST_F(PhysicsTestHost, calc_tabulated_physics_step) + const double expected_xs[] = {0.0001, 0.001, 0.1, 0.0001, 0.001, 0.1}; + EXPECT_VEC_SOFT_EQ(expected_xs, xs); +} + +TEST_F(PhysicsTrackViewHostTest, calc_range) { - ParticleTrackView particle(params.particle, states.particle, ThreadId{0}); - PhysicsTrackView physics(params.physics, - states.physics, - ParticleId{0}, - MaterialId{0}, - ThreadId{0}); - - // XXX add tests for a variety of energy ranges and multiple material IDs + // Default range and scaling + EXPECT_SOFT_EQ(0.1 * units::centimeter, params_ref.scaling_min_range); + EXPECT_SOFT_EQ(0.2, params_ref.scaling_fraction); + std::vector range; + std::vector step; + + // Range: increases with energy, constant with material + for (const char* particle : {"celeriton", "anti-celeriton"}) { - states.particle.state[ThreadId{0}].energy = MevEnergy{1e3}; - real_type step - = celeritas::calc_tabulated_physics_step(particle, physics); - EXPECT_EQ(0, step); + const PhysicsTrackView phys + = this->make_track_view(particle, MaterialId{0}); + auto meow_ppid = this->find_ppid(phys, "meows"); + ASSERT_TRUE(meow_ppid); + auto id = phys.value_grid(ValueGridType::range, meow_ppid); + ASSERT_TRUE(id); + auto calc_range = phys.make_calculator(id); + for (real_type energy : {0.01, 1.0, 1e2}) + { + range.push_back(calc_range(MevEnergy{energy})); + step.push_back(phys.range_to_step(range.back())); + } } + + const double expected_range[] = {0.04, 4, 40, 0.04, 4, 40}; + const double expected_step[] = {0.04, 0.958, 8.1598, 0.04, 0.958, 8.1598}; + EXPECT_VEC_SOFT_EQ(expected_range, range); + EXPECT_VEC_SOFT_EQ(expected_step, step); } -TEST_F(PhysicsTestHost, calc_energy_loss) {} +TEST_F(PhysicsTrackViewHostTest, model_finder) +{ + const PhysicsTrackView phys + = this->make_track_view("celeriton", MaterialId{0}); + auto purr_ppid = this->find_ppid(phys, "purrs"); + ASSERT_TRUE(purr_ppid); + auto find_model = phys.make_model_finder(purr_ppid); -TEST_F(PhysicsTestHost, select_model) {} + // See expected_model_names above + EXPECT_FALSE(find_model(MevEnergy{0.999e-3})); + EXPECT_EQ(3, find_model(MevEnergy{0.5}).unchecked_get()); + EXPECT_EQ(4, find_model(MevEnergy{5}).unchecked_get()); + EXPECT_EQ(5, find_model(MevEnergy{50}).unchecked_get()); + EXPECT_FALSE(find_model(MevEnergy{100.1})); +} diff --git a/test/physics/base/PhysicsStepUtils.test.cc b/test/physics/base/PhysicsStepUtils.test.cc new file mode 100644 index 0000000000..7096e48ea2 --- /dev/null +++ b/test/physics/base/PhysicsStepUtils.test.cc @@ -0,0 +1,75 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file PhysicsStepUtils.test.cc +//---------------------------------------------------------------------------// +#include "physics/base/PhysicsStepUtils.hh" + +#include "base/PieStateStore.hh" +#include "physics/base/ParticleParams.hh" +#include "physics/base/PhysicsParams.hh" +#include "celeritas_test.hh" +#include "PhysicsTestBase.hh" + +using namespace celeritas; +using namespace celeritas_test; +using celeritas::units::MevEnergy; + +//---------------------------------------------------------------------------// +// TEST HARNESS +//---------------------------------------------------------------------------// + +class PhysicsStepUtilsTest : public PhysicsTestBase +{ + using Base = PhysicsTestBase; + + protected: + using ParticleStateStore = PieStateStore; + using PhysicsStateStore = PieStateStore; + + void SetUp() override + { + Base::SetUp(); + + // Construct state for a single host thread + par_state = ParticleStateStore(*this->particles(), 1); + phys_state = PhysicsStateStore(*this->physics(), 1); + } + + ParticleStateStore par_state; + PhysicsStateStore phys_state; +}; + +//---------------------------------------------------------------------------// +// TESTS +//---------------------------------------------------------------------------// + +TEST_F(PhysicsStepUtilsTest, calc_tabulated_physics_step) +{ + ParticleTrackView particle( + this->particles()->host_pointers(), par_state.ref(), ThreadId{0}); + PhysicsTrackView phys(this->physics()->host_pointers(), + phys_state.ref(), + ParticleId{0}, + MaterialId{0}, + ThreadId{0}); + ParticleTrackView::Initializer_t par_init; + PhysicsTrackView::Initializer_t phys_init; + + // XXX add tests for a variety of energy ranges and multiple material IDs + { + par_init.energy = MevEnergy{1e3}; + par_init.particle_id = particles()->find("celeriton"); + particle = par_init; + phys = phys_init; + + real_type step = celeritas::calc_tabulated_physics_step(particle, phys); + EXPECT_EQ(0, step); + } +} + +TEST_F(PhysicsStepUtilsTest, calc_energy_loss) {} + +TEST_F(PhysicsStepUtilsTest, select_model) {} diff --git a/test/physics/base/PhysicsTestBase.cc b/test/physics/base/PhysicsTestBase.cc new file mode 100644 index 0000000000..e0a07e1a80 --- /dev/null +++ b/test/physics/base/PhysicsTestBase.cc @@ -0,0 +1,159 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file PhysicsTestBase.cc +//---------------------------------------------------------------------------// +#include "PhysicsTestBase.hh" + +#include "physics/material/MaterialParams.hh" +#include "physics/base/ParticleParams.hh" +#include "physics/base/PhysicsParams.hh" + +using namespace celeritas; + +namespace celeritas_test +{ +//---------------------------------------------------------------------------// +//! Virtual destructor +PhysicsTestBase::~PhysicsTestBase() = default; + +//---------------------------------------------------------------------------// +void PhysicsTestBase::SetUp() +{ + materials_ = this->build_materials(); + CELER_ASSERT(materials_); + particles_ = this->build_particles(); + CELER_ASSERT(particles_); + physics_ = this->build_physics(); + CELER_ASSERT(physics_); +} + +//---------------------------------------------------------------------------// +auto PhysicsTestBase::build_materials() const -> SPConstMaterials +{ + using namespace celeritas::units; + MaterialParams::Input inp; + inp.elements + = {{1, AmuMass{1.0}, "celerogen"}, {4, AmuMass{10.0}, "celerinium"}}; + inp.materials.push_back({1e20, + 300, + MatterState::gas, + {{ElementId{0}, 1.0}}, + "lo density celerogen"}); + inp.materials.push_back({1e21, + 300, + MatterState::liquid, + {{ElementId{0}, 1.0}}, + "hi density celerogen"}); + inp.materials.push_back({1e23, + 300, + MatterState::solid, + {{ElementId{1}, 1.0}}, + "solid celerinium"}); + return std::make_shared(std::move(inp)); +} + +//---------------------------------------------------------------------------// +auto PhysicsTestBase::build_particles() const -> SPConstParticles +{ + using namespace celeritas::units; + namespace pdg = celeritas::pdg; + + constexpr auto zero = celeritas::zero_quantity(); + constexpr auto stable = celeritas::ParticleDef::stable_decay_constant(); + + ParticleParams::Input inp; + inp.push_back({"gamma", pdg::gamma(), zero, zero, stable}); + inp.push_back( + {"celeriton", PDGNumber{1337}, MevMass{1}, ElementaryCharge{1}, stable}); + inp.push_back({"anti-celeriton", + PDGNumber{-1337}, + MevMass{1}, + ElementaryCharge{-1}, + stable}); + return std::make_shared(std::move(inp)); +} + +//---------------------------------------------------------------------------// +auto PhysicsTestBase::build_physics() const -> SPConstPhysics +{ + using Barn = MockProcess::BarnMicroXs; + PhysicsParams::Input physics_inp; + physics_inp.materials = this->materials(); + physics_inp.particles = this->particles(); + + // Add a few processes + MockProcess::Input inp; + inp.materials = this->materials(); + inp.interact = this->make_model_callback(); + { + inp.label = "scattering"; + inp.applic = {make_applicability("gamma", 1e-6, 100), + make_applicability("celeriton", 1, 100)}; + inp.xs = Barn{1.0}; + inp.energy_loss = {}; + inp.range = {}; + physics_inp.processes.push_back(std::make_shared(inp)); + } + { + inp.label = "absorption"; + inp.applic = {make_applicability("gamma", 1e-6, 100)}; + inp.xs = Barn{2.0}; + inp.energy_loss = {}; + inp.range = {}; + physics_inp.processes.push_back(std::make_shared(inp)); + } + { + // Three different models for the single process + inp.label = "purrs"; + inp.applic = {make_applicability("celeriton", 1e-3, 1), + make_applicability("celeriton", 1, 10), + make_applicability("celeriton", 10, 100)}; + inp.xs = Barn{3.0}; + inp.energy_loss = 0.2; + inp.range = 2; + physics_inp.processes.push_back(std::make_shared(inp)); + } + { + // Two models for anti-celeriton + inp.label = "hisses"; + inp.applic = {make_applicability("anti-celeriton", 1e-3, 1), + make_applicability("anti-celeriton", 1, 100)}; + inp.xs = Barn{4.0}; + inp.energy_loss = 0.3; + inp.range = 3; + physics_inp.processes.push_back(std::make_shared(inp)); + } + { + inp.label = "meows"; + inp.applic = {make_applicability("celeriton", 1e-3, 10), + make_applicability("anti-celeriton", 1e-3, 10)}; + inp.xs = Barn{5.0}; + inp.energy_loss = 0.4; + inp.range = 4; + physics_inp.processes.push_back(std::make_shared(inp)); + } + return std::make_shared(std::move(physics_inp)); +} + +//---------------------------------------------------------------------------// +auto PhysicsTestBase::make_applicability(const char* name, + double lo_energy, + double hi_energy) const + -> Applicability +{ + using celeritas::units::MevEnergy; + CELER_EXPECT(particles_); + CELER_EXPECT(name); + CELER_EXPECT(lo_energy <= hi_energy); + Applicability result; + result.particle = particles_->find(name); + result.lower = MevEnergy{lo_energy}; + result.upper = MevEnergy{hi_energy}; + return result; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas_test diff --git a/test/physics/base/PhysicsTestBase.hh b/test/physics/base/PhysicsTestBase.hh new file mode 100644 index 0000000000..71aae1687c --- /dev/null +++ b/test/physics/base/PhysicsTestBase.hh @@ -0,0 +1,92 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file PhysicsTestBase.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "gtest/Test.hh" + +#include "physics/material/MaterialParams.hh" +#include "physics/base/ParticleParams.hh" +#include "MockProcess.hh" + +namespace celeritas +{ +class MaterialParams; +class ParticleParams; +class PhysicsParams; +} // namespace celeritas + +namespace celeritas_test +{ +//---------------------------------------------------------------------------// +/*! + * Construct mock materials, particles, and physics data. + * + * This class creates two elements, three single-element materials, three + * particles, and five MockProcesses, each of which emits one or more + * MockModels. + * - gamma:scattering + * - gamma:absorption + * - celeriton:scattering + * - celeriton:purrs + * - celeriton:meows + * - anti-celeriton:hisses + * - anti-celeriton:meows + */ +class PhysicsTestBase : public celeritas::Test +{ + protected: + //!@{ + //! Type aliases + using SPConstMaterials = std::shared_ptr; + using SPConstParticles = std::shared_ptr; + using SPConstPhysics = std::shared_ptr; + using Applicability = celeritas::Applicability; + using ModelId = celeritas::ModelId; + using ModelCallback = std::function; + //!@} + + protected: + virtual ~PhysicsTestBase(); + + void SetUp() override; + + virtual SPConstMaterials build_materials() const; + virtual SPConstParticles build_particles() const; + virtual SPConstPhysics build_physics() const; + + const SPConstMaterials& materials() const { return materials_; } + const SPConstParticles& particles() const { return particles_; } + const SPConstPhysics& physics() const { return physics_; } + + public: + Applicability make_applicability(const char* name, + double lo_energy, + double hi_energy) const; + + ModelCallback make_model_callback() const + { + return [this](ModelId id) { interactions_.push_back(id); }; + } + + celeritas::Span called_models() const + { + return make_span(interactions_); + } + + private: + //// DATA //// + + SPConstMaterials materials_; + SPConstParticles particles_; + SPConstPhysics physics_; + + mutable std::vector interactions_; +}; + +//---------------------------------------------------------------------------// +} // namespace celeritas_test From ea77e160039fe8cb09c5ef60a55d2d1def79c4cf Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Thu, 11 Feb 2021 21:59:39 -0500 Subject: [PATCH 5/6] Add CUDA test --- scripts/dev/celeritas-gen.py | 2 + src/physics/base/PhysicsInterface.hh | 8 +- test/CMakeLists.txt | 2 +- test/base/Pie.test.cc | 4 +- test/base/Pie.test.cu | 6 +- test/base/Pie.test.hh | 12 ++- test/physics/base/Physics.test.cc | 129 ++++++++++++++++++++++++++- test/physics/base/Physics.test.cu | 51 +++++++++++ test/physics/base/Physics.test.hh | 116 ++++++++++++++++++++++++ 9 files changed, 316 insertions(+), 14 deletions(-) create mode 100644 test/physics/base/Physics.test.cu create mode 100644 test/physics/base/Physics.test.hh diff --git a/scripts/dev/celeritas-gen.py b/scripts/dev/celeritas-gen.py index e679847e24..ff1a5cca57 100755 --- a/scripts/dev/celeritas-gen.py +++ b/scripts/dev/celeritas-gen.py @@ -179,6 +179,8 @@ class {name}Test : public celeritas::Test {lowabbr}_test_kernel<<>>( input.num_threads); + CELER_CUDA_CHECK_ERROR(); + {capabbr}TestOutput result; return result; }} diff --git a/src/physics/base/PhysicsInterface.hh b/src/physics/base/PhysicsInterface.hh index 016dcfa6ab..a94d330f61 100644 --- a/src/physics/base/PhysicsInterface.hh +++ b/src/physics/base/PhysicsInterface.hh @@ -273,10 +273,10 @@ struct PhysicsStateData * Resize a material state in host code. */ template -inline void -resize(PhysicsStateData* data, - const PhysicsParamsData& params, - size_type size) +inline void resize( + PhysicsStateData* data, + const PhysicsParamsData& params, + size_type size) { CELER_EXPECT(size > 0); CELER_EXPECT(params.max_particle_processes > 0); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c46a47dde5..50c86bc234 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -172,7 +172,7 @@ set(CELERITASTEST_LINK_LIBRARIES CeleritasPhysicsTest) celeritas_setup_tests(SERIAL PREFIX physics/base) celeritas_cudaoptional_test(physics/base/Particle) -celeritas_add_test(physics/base/Physics.test.cc) +celeritas_cudaoptional_test(physics/base/Physics) celeritas_add_test(physics/base/PhysicsStepUtils.test.cc ${_not_impl}) celeritas_setup_tests(SERIAL PREFIX physics/grid) diff --git a/test/base/Pie.test.cc b/test/base/Pie.test.cc index 82b3ef5440..c4218534f0 100644 --- a/test/base/Pie.test.cc +++ b/test/base/Pie.test.cc @@ -183,9 +183,7 @@ TEST_F(PieTest, TEST_IF_CELERITAS_CUDA(device)) kernel_input.states = device_states; kernel_input.result = device_result.device_pointers(); -#if CELERITAS_USE_CUDA - p_test(kernel_input); -#endif + pie_cuda_test(kernel_input); std::vector result(device_result.size()); device_result.copy_to_host(celeritas::make_span(result)); diff --git a/test/base/Pie.test.cu b/test/base/Pie.test.cu index b872d7a9a1..4df3cd2256 100644 --- a/test/base/Pie.test.cu +++ b/test/base/Pie.test.cu @@ -16,7 +16,7 @@ namespace //---------------------------------------------------------------------------// // KERNELS //---------------------------------------------------------------------------// -__global__ void p_test_kernel( +__global__ void pie_cuda_test_kernel( const MockParamsPies params, const MockStatePies states, const celeritas::Span results) @@ -62,11 +62,11 @@ __global__ void p_test_kernel( // TESTING INTERFACE //---------------------------------------------------------------------------// //! Run on device and return results -void p_test(PTestInput input) +void pie_cuda_test(PTestInput input) { celeritas::KernelParamCalculator calc_launch_params; auto params = calc_launch_params(input.states.size()); - p_test_kernel<<>>( + pie_cuda_test_kernel<<>>( input.params, input.states, input.result); CELER_CUDA_CHECK_ERROR(); diff --git a/test/base/Pie.test.hh b/test/base/Pie.test.hh index f14ceb9c31..92bb53de0a 100644 --- a/test/base/Pie.test.hh +++ b/test/base/Pie.test.hh @@ -5,9 +5,10 @@ //---------------------------------------------------------------------------// //! \file Pie.test.hh //---------------------------------------------------------------------------// - +#include "base/Assert.hh" #include "base/Pie.hh" #include "base/Types.hh" +#include "celeritas_config.h" namespace celeritas_test { @@ -166,7 +167,14 @@ struct PTestInput //---------------------------------------------------------------------------// //! Run on device and return results -void p_test(PTestInput); +void pie_cuda_test(PTestInput); + +#if !CELERITAS_USE_CUDA +inline void pie_cuda_test(PTestInput) +{ + CELER_NOT_CONFIGURED("CUDA"); +} +#endif //---------------------------------------------------------------------------// } // namespace celeritas_test diff --git a/test/physics/base/Physics.test.cc b/test/physics/base/Physics.test.cc index 4ac0d9f91c..e26073c8ea 100644 --- a/test/physics/base/Physics.test.cc +++ b/test/physics/base/Physics.test.cc @@ -12,10 +12,13 @@ #include "base/Range.hh" #include "base/PieStateStore.hh" #include "physics/base/ParticleParams.hh" + #include "PhysicsTestBase.hh" +#include "Physics.test.hh" using namespace celeritas; using namespace celeritas_test; +using MevEnergy = celeritas::units::MevEnergy; //---------------------------------------------------------------------------// // PHYSICS BASE CLASS @@ -98,7 +101,6 @@ class PhysicsTrackViewHostTest : public PhysicsTestBase using StateStore = PieStateStore; using ParamsHostRef = PhysicsParamsData; - using MevEnergy = celeritas::units::MevEnergy; //!@} void SetUp() override @@ -329,3 +331,128 @@ TEST_F(PhysicsTrackViewHostTest, model_finder) EXPECT_EQ(5, find_model(MevEnergy{50}).unchecked_get()); EXPECT_FALSE(find_model(MevEnergy{100.1})); } + +TEST_F(PhysicsTrackViewHostTest, cuda_surrogate) +{ + std::vector step; + for (const char* particle : {"gamma", "anti-celeriton"}) + { + PhysicsTrackView phys = this->make_track_view(particle, MaterialId{1}); + + for (real_type energy : {1e-5, 1e-3, 1., 100., 1e5}) + { + step.push_back(celeritas_test::calc_step(phys, MevEnergy{energy})); + } + } + + // PRINT_EXPECTED(step); + const double expected_step[] = {166.6666666667, + 166.6666666667, + 166.6666666667, + 166.6666666667, + inf, + 0.003, + 0.003, + 0.7573333333333, + 8.1598, + 8.1598}; + EXPECT_VEC_SOFT_EQ(expected_step, step); +} + +//---------------------------------------------------------------------------// +// PHYSICS TRACK VIEW (DEVICE) +//---------------------------------------------------------------------------// + +#define PHYS_DEVICE_TEST TEST_IF_CELERITAS_CUDA(PhysicsTrackViewDeviceTest) +class PHYS_DEVICE_TEST : public PhysicsTestBase +{ + using Base = PhysicsTestBase; + + public: + //!@{ + //! Type aliases + using StateStore = PieStateStore; + //!@} + + void SetUp() override + { + Base::SetUp(); + + CELER_ASSERT(this->physics()); + } + + StateStore states; + celeritas::StatePie inits; +}; + +TEST_F(PHYS_DEVICE_TEST, all) +{ + // Construct initial conditions + { + celeritas::StatePie + temp_inits; + + auto init_builder = make_pie_builder(&temp_inits); + PhysTestInit thread_init; + for (unsigned int matid : {0, 2}) + { + thread_init.mat = MaterialId{matid}; + for (real_type energy : {1e-5, 1e-3, 1., 100., 1e5}) + { + thread_init.energy = MevEnergy{energy}; + for (unsigned int particle : {0, 1, 2}) + { + thread_init.particle = ParticleId{particle}; + init_builder.push_back(thread_init); + } + } + } + this->inits = temp_inits; + } + + states = StateStore(*this->physics(), this->inits.size()); + celeritas::DeviceVector step(this->states.size()); + + PTestInput inp; + inp.params = this->physics()->device_pointers(); + inp.states = this->states.ref(); + inp.inits = this->inits; + inp.result = step.device_pointers(); + + phys_cuda_test(inp); + std::vector step_host(step.size()); + step.copy_to_host(make_span(step_host)); + // clang-format on + const double expected_step_host[] = {1666.666666667, + 0.002, + 0.003, + 1666.666666667, + 0.002, + 0.003, + 1666.666666667, + 0.556, + 0.7573333333333, + 1666.666666667, + 8.1598, + 8.1598, + inf, + 8.1598, + 8.1598, + 1.666666666667, + 0.002, + 0.003, + 1.666666666667, + 0.002, + 0.003, + 1.666666666667, + 0.5555555555556, + 0.5555555555556, + 1.666666666667, + 1.25, + 1.25, + inf, + 8.1598, + 8.1598}; + // clang-format off + EXPECT_VEC_SOFT_EQ(expected_step_host, step_host); +} diff --git a/test/physics/base/Physics.test.cu b/test/physics/base/Physics.test.cu new file mode 100644 index 0000000000..8e17b3aa61 --- /dev/null +++ b/test/physics/base/Physics.test.cu @@ -0,0 +1,51 @@ +//---------------------------------*-CUDA-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Physics.test.cu +//---------------------------------------------------------------------------// +#include "Physics.test.hh" + +#include "base/KernelParamCalculator.cuda.hh" + +using namespace celeritas; + +namespace celeritas_test +{ +namespace +{ +//---------------------------------------------------------------------------// +// KERNELS +//---------------------------------------------------------------------------// + +__global__ void phys_test_kernel(const PTestInput inp) +{ + auto tid = celeritas::KernelParamCalculator::thread_id(); + if (tid.get() >= inp.states.size()) + return; + + const auto& init = inp.inits[tid]; + PhysicsTrackView phys(inp.params, inp.states, init.particle, init.mat, tid); + + phys = PhysicsTrackInitializer{}; + inp.result[tid.get()] = calc_step(phys, init.energy); +} +} // namespace + +//---------------------------------------------------------------------------// +// TESTING INTERFACE +//---------------------------------------------------------------------------// +//! Run on device and return results +void phys_cuda_test(const PTestInput& input) +{ + CELER_ASSERT(input.inits.size() == input.states.size()); + KernelParamCalculator calc_launch_params; + auto params = calc_launch_params(input.states.size()); + phys_test_kernel<<>>(input); + + CELER_CUDA_CHECK_ERROR(); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas_test diff --git a/test/physics/base/Physics.test.hh b/test/physics/base/Physics.test.hh new file mode 100644 index 0000000000..3244c2ab18 --- /dev/null +++ b/test/physics/base/Physics.test.hh @@ -0,0 +1,116 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file Physics.test.hh +//---------------------------------------------------------------------------// +#include "celeritas_config.h" +#include "base/Assert.hh" +#include "base/Pie.hh" +#include "base/Range.hh" +#include "base/Span.hh" +#include "physics/base/PhysicsInterface.hh" +#include "physics/base/Units.hh" +#include "physics/base/Types.hh" +#include "physics/material/Types.hh" + +// Kernel functions +#include "base/Algorithms.hh" +#include "base/NumericLimits.hh" +#include "physics/base/PhysicsTrackView.hh" + +namespace celeritas_test +{ +using celeritas::MemSpace; +using celeritas::Ownership; +//---------------------------------------------------------------------------// +// TESTING INTERFACE +//---------------------------------------------------------------------------// +//! Input data +struct PhysTestInit +{ + celeritas::units::MevEnergy energy; + celeritas::MaterialId mat; + celeritas::ParticleId particle; +}; + +struct PTestInput +{ + celeritas::PhysicsParamsData + params; + celeritas::PhysicsStateData states; + celeritas::StatePie + inits; + + // Calculated "step" per track + celeritas::Span result; +}; + +//---------------------------------------------------------------------------// +// HELPER FUNCTIONS +//---------------------------------------------------------------------------// +inline CELER_FUNCTION celeritas::real_type + calc_step(celeritas::PhysicsTrackView& phys, celeritas::units::MevEnergy energy) +{ + using namespace celeritas; + + // Calc total macro_xs over processsess + real_type total_xs = 0; + for (auto pp_idx : range(phys.num_particle_processes())) + { + ParticleProcessId ppid{pp_idx}; + real_type process_xs = 0; + if (auto id = phys.value_grid(ValueGridType::macro_xs, ppid)) + { + auto calc_xs = phys.make_calculator(id); + process_xs = calc_xs(energy); + } + + // Zero cross section if outside of model range + auto find_model = phys.make_model_finder(ppid); + if (!find_model(energy)) + { + process_xs = 0; + } + + phys.per_process_xs(ppid) = process_xs; + total_xs += process_xs; + } + phys.interaction_mfp(1 / total_xs); + + // Calc minimum range + const auto inf = numeric_limits::infinity(); + real_type step = inf; + for (auto pp_idx : range(phys.num_particle_processes())) + { + ParticleProcessId ppid{pp_idx}; + if (auto id = phys.value_grid(ValueGridType::range, ppid)) + { + auto calc_range = phys.make_calculator(id); + step = min(step, calc_range(energy)); + } + } + if (step != inf) + { + step = phys.range_to_step(step); + } + + // Take minimum of step and half the MFP + step = min(step, 0.5 * phys.interaction_mfp()); + return step; +} + +//---------------------------------------------------------------------------// +//! Run on device and return results +void phys_cuda_test(const PTestInput&); + +#if !CELERITAS_USE_CUDA +inline void phys_cuda_test(const PTestInput&) +{ + CELER_NOT_CONFIGURED("CUDA"); +} +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas_test From 63f4ecb25816ac035b0a90a563f5e97e5f6de6d9 Mon Sep 17 00:00:00 2001 From: Seth R Johnson Date: Fri, 12 Feb 2021 07:21:32 -0500 Subject: [PATCH 6/6] Update comments --- src/physics/base/PhysicsParams.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/physics/base/PhysicsParams.hh b/src/physics/base/PhysicsParams.hh index f52b25c676..066febcf35 100644 --- a/src/physics/base/PhysicsParams.hh +++ b/src/physics/base/PhysicsParams.hh @@ -37,8 +37,8 @@ class ParticleParams; * Input options are: * - \c min_range: below this value, there is no extra transformation from * particle range to step length. - * - \c max_step_over_range: at higher energies (longer ranges), gradually - * decrease the minimum step length until it's this fraction of the tabulated + * - \c max_step_over_range: at higher energy (longer range), gradually + * decrease the maximum step length until it's this fraction of the tabulated * range. */ class PhysicsParams