diff --git a/app/celer-dump-data.cc b/app/celer-dump-data.cc index 2d9c5601e3..3ed11c94a9 100644 --- a/app/celer-dump-data.cc +++ b/app/celer-dump-data.cc @@ -123,8 +123,8 @@ void print_isotopes(std::vector& isotopes) cout << R"gfm( # Isotopes -| Isotope ID | Name | Atomic number | Atomic mass number | Binding energy (MeV) | Nuclear mass (MeV) | -| ---------- | ------ | ------------- | ------------------ | -------------------- | ------------------ | +| Isotope ID | Name | Atomic number | Atomic mass number | Binding energy (MeV) | Proton loss energy (MeV) | Neutron loss energy (MeV) | Nuclear mass (MeV) | +| ---------- | ------ | ------------- | ------------------ | -------------------- | ------------------------ | ------------------------- | ------------------ | )gfm"; for (unsigned int isotope_id : range(isotopes.size())) @@ -137,6 +137,8 @@ void print_isotopes(std::vector& isotopes) << setw(13) << isotope.atomic_number << " | " << setw(18) << isotope.atomic_mass_number << " | " << setw(20) << isotope.binding_energy << " | " + << setw(24) << isotope.proton_loss_energy << " | " + << setw(25) << isotope.neutron_loss_energy << " | " << setw(18) << isotope.nuclear_mass << " |\n"; // clang-format on } diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index 04084ddd52..0325a17b72 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -62,6 +62,7 @@ list(APPEND SOURCES mat/MaterialParams.cc mat/MaterialParamsOutput.cc mat/detail/Utils.cc + neutron/model/CascadeOptions.cc neutron/process/NeutronElasticProcess.cc neutron/process/NeutronInelasticProcess.cc optical/CerenkovParams.cc @@ -151,6 +152,7 @@ if(CELERITAS_USE_JSON) ext/GeantPhysicsOptionsIO.json.cc field/FieldDriverOptionsIO.json.cc field/RZMapFieldInputIO.json.cc + neutron/model/CascadeOptionsIO.json.cc phys/PrimaryGeneratorOptionsIO.json.cc user/RootStepWriterIO.json.cc ) diff --git a/src/celeritas/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index 4bcb2211d5..905aa35e73 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -425,6 +425,29 @@ std::vector import_isotopes() isotope.atomic_mass_number = g4isotope.GetN(); isotope.binding_energy = G4NucleiProperties::GetBindingEnergy( isotope.atomic_mass_number, isotope.atomic_number); + + // Binding energy difference for losing a nucleon + if (isotope.atomic_mass_number > 1 && isotope.atomic_number > 1 + && isotope.atomic_mass_number >= isotope.atomic_number) + { + isotope.proton_loss_energy + = G4NucleiProperties::GetBindingEnergy( + isotope.atomic_mass_number, isotope.atomic_number) + - G4NucleiProperties::GetBindingEnergy( + isotope.atomic_mass_number - 1, + isotope.atomic_number - 1); + isotope.neutron_loss_energy + = G4NucleiProperties::GetBindingEnergy( + isotope.atomic_mass_number, isotope.atomic_number) + - G4NucleiProperties::GetBindingEnergy( + isotope.atomic_mass_number - 1, isotope.atomic_number); + } + else + { + isotope.proton_loss_energy = 0; + isotope.neutron_loss_energy = 0; + } + isotope.nuclear_mass = G4NucleiProperties::GetNuclearMass( isotope.atomic_mass_number, isotope.atomic_number); } diff --git a/src/celeritas/io/ImportElement.hh b/src/celeritas/io/ImportElement.hh index a50d77be15..845c77e828 100644 --- a/src/celeritas/io/ImportElement.hh +++ b/src/celeritas/io/ImportElement.hh @@ -22,10 +22,12 @@ namespace celeritas struct ImportIsotope { std::string name; //!< Isotope label - int atomic_number; //!< Atomic number Z - int atomic_mass_number; //!< Atomic number A - double binding_energy; //!< Nuclear binding energy [MeV] - double nuclear_mass; //!< Sum of nucleons' mass + binding energy [MeV] + int atomic_number{0}; //!< Atomic number Z + int atomic_mass_number{0}; //!< Atomic number A + double binding_energy{0}; //!< Nuclear binding energy (BE) [MeV] + double proton_loss_energy{0}; //!< BE(A, Z) - BE(A-1, Z-1) [MeV] + double neutron_loss_energy{0}; //!< BE(A, Z) - BE(A-1, Z) [MeV] + double nuclear_mass{0}; //!< Sum of nucleons' mass + binding energy [MeV] }; //---------------------------------------------------------------------------// diff --git a/src/celeritas/mat/IsotopeView.hh b/src/celeritas/mat/IsotopeView.hh index 11a4eff637..5a88f82174 100644 --- a/src/celeritas/mat/IsotopeView.hh +++ b/src/celeritas/mat/IsotopeView.hh @@ -50,6 +50,12 @@ class IsotopeView // Nuclear binding energy CELER_FORCEINLINE_FUNCTION MevEnergy binding_energy() const; + // Nuclear binding energy difference for a proton loss + CELER_FORCEINLINE_FUNCTION MevEnergy proton_loss_energy() const; + + // Nuclear binding energy difference for a neutron loss + CELER_FORCEINLINE_FUNCTION MevEnergy neutron_loss_energy() const; + // Sum of nucleons + binding energy CELER_FORCEINLINE_FUNCTION MevMass nuclear_mass() const; @@ -115,6 +121,24 @@ CELER_FUNCTION units::MevEnergy IsotopeView::binding_energy() const return isotope_def().binding_energy; } +//---------------------------------------------------------------------------// +/*! + * Nuclear binding energy difference for a proton loss. + */ +CELER_FUNCTION units::MevEnergy IsotopeView::proton_loss_energy() const +{ + return isotope_def().proton_loss_energy; +} + +//---------------------------------------------------------------------------// +/*! + * Nuclear binding energy difference for a neutron loss. + */ +CELER_FUNCTION units::MevEnergy IsotopeView::neutron_loss_energy() const +{ + return isotope_def().neutron_loss_energy; +} + //---------------------------------------------------------------------------// /*! * Nuclear mass, the sum of the nucleons' mass and their binding energy. diff --git a/src/celeritas/mat/MaterialData.hh b/src/celeritas/mat/MaterialData.hh index 7f0bcc5986..1aef02833d 100644 --- a/src/celeritas/mat/MaterialData.hh +++ b/src/celeritas/mat/MaterialData.hh @@ -33,7 +33,9 @@ struct IsotopeRecord AtomicNumber atomic_number; //!< Atomic number Z AtomicMassNumber atomic_mass_number; //!< Atomic number A - units::MevEnergy binding_energy; //!< Nuclear binding energy + units::MevEnergy binding_energy; //!< Nuclear binding energy (BE) + units::MevEnergy proton_loss_energy; //!< BE(A, Z) - BE(A-1, Z-1) + units::MevEnergy neutron_loss_energy; //!< BE(A, Z) - BE(A-1, Z) units::MevMass nuclear_mass; //!< Nucleons' mass + binding energy }; diff --git a/src/celeritas/mat/MaterialParams.cc b/src/celeritas/mat/MaterialParams.cc index 289c1e72ea..1627940c5b 100644 --- a/src/celeritas/mat/MaterialParams.cc +++ b/src/celeritas/mat/MaterialParams.cc @@ -78,6 +78,10 @@ MaterialParams::from_import(ImportData const& data) = AtomicNumber{isotope.atomic_mass_number}; isotope_params.binding_energy = units::MevEnergy(isotope.binding_energy); + isotope_params.proton_loss_energy + = units::MevEnergy(isotope.proton_loss_energy); + isotope_params.neutron_loss_energy + = units::MevEnergy(isotope.neutron_loss_energy); // Convert from MeV (Geant4) to MeV/c^2 (Celeritas) isotope_params.nuclear_mass = units::MevMass(isotope.nuclear_mass); @@ -373,6 +377,8 @@ void MaterialParams::append_isotope_def(IsotopeInput const& inp, CELER_EXPECT(inp.atomic_number); CELER_EXPECT(inp.atomic_mass_number); CELER_EXPECT(inp.binding_energy >= zero_quantity()); + CELER_EXPECT(inp.proton_loss_energy >= zero_quantity()); + CELER_EXPECT(inp.neutron_loss_energy >= zero_quantity()); CELER_EXPECT(inp.nuclear_mass > zero_quantity()); IsotopeRecord result; @@ -381,6 +387,8 @@ void MaterialParams::append_isotope_def(IsotopeInput const& inp, result.atomic_number = inp.atomic_number; result.atomic_mass_number = inp.atomic_mass_number; result.binding_energy = inp.binding_energy; + result.proton_loss_energy = inp.proton_loss_energy; + result.neutron_loss_energy = inp.neutron_loss_energy; result.nuclear_mass = inp.nuclear_mass; // Add to host vector diff --git a/src/celeritas/mat/MaterialParams.hh b/src/celeritas/mat/MaterialParams.hh index 1cb4b0e557..3d91e53c11 100644 --- a/src/celeritas/mat/MaterialParams.hh +++ b/src/celeritas/mat/MaterialParams.hh @@ -53,11 +53,14 @@ class MaterialParams final : public ParamsDataInterface //!@{ //! \name Type aliases using AtomicMassNumber = AtomicNumber; + using MevEnergy = units::MevEnergy; //!@} AtomicNumber atomic_number; //!< Atomic number Z AtomicMassNumber atomic_mass_number; //!< Atomic number A - units::MevEnergy binding_energy; //!< Nuclear binding energy + MevEnergy binding_energy; //!< Nuclear binding energy (BE) + MevEnergy proton_loss_energy; //!< BE(A, Z) - BE(A-1, Z-1) + MevEnergy neutron_loss_energy; //!< BE(A, Z) - BE(A-1, Z) units::MevMass nuclear_mass; //!< Nucleons' mass + binding energy Label label; //!< Isotope name }; diff --git a/src/celeritas/mat/MaterialParamsOutput.cc b/src/celeritas/mat/MaterialParamsOutput.cc index ab24baabb0..1c23770339 100644 --- a/src/celeritas/mat/MaterialParamsOutput.cc +++ b/src/celeritas/mat/MaterialParamsOutput.cc @@ -54,6 +54,8 @@ void MaterialParamsOutput::output(JsonPimpl* j) const auto atomic_number = json::array(); auto atomic_mass_number = json::array(); auto binding_energy = json::array(); + auto proton_loss_energy = json::array(); + auto neutron_loss_energy = json::array(); auto nuclear_mass = json::array(); for (auto id : range(IsotopeId{material_->num_isotopes()})) @@ -64,6 +66,9 @@ void MaterialParamsOutput::output(JsonPimpl* j) const atomic_mass_number.push_back( iso_view.atomic_mass_number().unchecked_get()); binding_energy.push_back(iso_view.binding_energy().value()); + proton_loss_energy.push_back(iso_view.proton_loss_energy().value()); + neutron_loss_energy.push_back( + iso_view.neutron_loss_energy().value()); nuclear_mass.push_back(iso_view.nuclear_mass().value()); } obj["isotopes"] = { @@ -71,10 +76,16 @@ void MaterialParamsOutput::output(JsonPimpl* j) const {"atomic_number", std::move(atomic_number)}, {"atomic_mass_number", std::move(atomic_mass_number)}, {"binding_energy", std::move(binding_energy)}, + {"proton_loss_energy", std::move(proton_loss_energy)}, + {"neutron_loss_energy", std::move(neutron_loss_energy)}, {"nuclear_mass", std::move(nuclear_mass)}, }; units["binding_energy"] = accessor_unit_label(); + units["proton_loss_energy"] + = accessor_unit_label(); + units["neutron_loss_energy"] + = accessor_unit_label(); units["nuclear_mass"] = accessor_unit_label(); } diff --git a/src/celeritas/neutron/data/NeutronInelasticData.hh b/src/celeritas/neutron/data/NeutronInelasticData.hh index 45aedb11a8..0010dee34e 100644 --- a/src/celeritas/neutron/data/NeutronInelasticData.hh +++ b/src/celeritas/neutron/data/NeutronInelasticData.hh @@ -26,9 +26,11 @@ struct NeutronInelasticScalars // Action and particle IDs ActionId action_id; ParticleId neutron_id; + ParticleId proton_id; // Particle mass * c^2 [MeV] units::MevMass neutron_mass; + units::MevMass proton_mass; //! Number of nucleon-nucleon channels static CELER_CONSTEXPR_FUNCTION size_type num_channels() { return 3; } @@ -43,7 +45,9 @@ struct NeutronInelasticScalars //! Whether data are assigned explicit CELER_FUNCTION operator bool() const { - return action_id && neutron_id && neutron_mass > zero_quantity(); + return action_id && neutron_id && proton_id + && neutron_mass > zero_quantity() + && neutron_mass > zero_quantity(); } }; @@ -59,6 +63,69 @@ struct StepanovParameters Real3 coeffs; //!< coefficients of a second order Stepanov's function }; +//---------------------------------------------------------------------------// +/*! + * Components of nuclear zone properties of the Bertini cascade model. + */ +struct ZoneComponent +{ + using NucleonArray = Array; //!< [proton, neutron] + + real_type radius{}; //!< radius of zones in [femtometer] + real_type volume{}; //!< volume of zones in [femtometer^3] + NucleonArray density{0, 0}; //!< nucleon densities [1/femtometer^3] + NucleonArray fermi_mom{0, 0}; //!< fermi momenta in [MeV/c] + NucleonArray potential{0, 0}; //!< nucleon potentials [MeV] +}; + +//---------------------------------------------------------------------------// +/*! + * Data characterizing the nuclear zones. + */ +struct NuclearZones +{ + ItemRange zones; + + //! Whether all data are assigned and valid + explicit CELER_FUNCTION operator bool() const { return !zones.empty(); } +}; + +//---------------------------------------------------------------------------// +/*! + * Device data for nuclear zone properties + */ +template +struct NuclearZoneData +{ + template + using Items = Collection; + template + using IsotopeItems = Collection; + + //// MEMBER DATA //// + + // Nuclear zone data + Items components; + IsotopeItems zones; + + //! Whether the data are assigned + explicit CELER_FUNCTION operator bool() const + { + return !components.empty() && !zones.empty(); + } + + //! Assign from another set of data + template + NuclearZoneData& operator=(NuclearZoneData const& other) + { + CELER_EXPECT(other); + components = other.components; + zones = other.zones; + + return *this; + } +}; + //---------------------------------------------------------------------------// /*! * Device data for creating an interactor. @@ -90,11 +157,14 @@ struct NeutronInelasticData // Backend data Items reals; + // Nuclear zone data + NuclearZoneData nuclear_zones; + //! Whether the data are assigned explicit CELER_FUNCTION operator bool() const { return scalars && !micro_xs.empty() && !nucleon_xs.empty() - && !xs_params.empty() && !reals.empty(); + && !xs_params.empty() && !reals.empty() && nuclear_zones; } //! Assign from another set of data @@ -107,6 +177,7 @@ struct NeutronInelasticData nucleon_xs = other.nucleon_xs; xs_params = other.xs_params; reals = other.reals; + nuclear_zones = other.nuclear_zones; return *this; } }; diff --git a/src/celeritas/neutron/model/CascadeOptions.cc b/src/celeritas/neutron/model/CascadeOptions.cc new file mode 100644 index 0000000000..53de19fd12 --- /dev/null +++ b/src/celeritas/neutron/model/CascadeOptions.cc @@ -0,0 +1,47 @@ +//----------------------------------*-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/neutron/model/CascadeOptions.cc +//---------------------------------------------------------------------------// +#include "CascadeOptions.hh" + +#include "corecel/Assert.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Throw a runtime assertion if any of the input is invalid. + */ +void validate_input(CascadeOptions const& opts) +{ + CELER_VALIDATE(opts.prob_pion_absorption >= 0, + << "invalid prob_pion_absorption " + << opts.prob_pion_absorption); + CELER_VALIDATE(opts.radius_scale > 0, + << "invalid radius_scale " << opts.radius_scale); + CELER_VALIDATE(opts.radius_small > 0, + << "invalid radius_small " << opts.radius_small); + CELER_VALIDATE(opts.radius_alpha > 0, + << "invalid radius_alpha " << opts.radius_alpha); + CELER_VALIDATE(opts.radius_trailing >= 0, + << "invalid radius_trailing " << opts.radius_trailing); + CELER_VALIDATE(opts.fermi_scale > 0, + << "invalid fermi_scale " << opts.fermi_scale); + CELER_VALIDATE(opts.xsec_scale > 0, + << "invalid xsec_scale " << opts.xsec_scale); + CELER_VALIDATE(opts.gamma_qd_scale > 0, + << "invalid gamma_qd_scale " << opts.gamma_qd_scale); + CELER_VALIDATE(opts.dp_max_doublet > 0, + << "invalid dp_max_doublet " << opts.dp_max_doublet); + CELER_VALIDATE(opts.dp_max_triplet > 0, + << "invalid dp_max_triplet " << opts.dp_max_triplet); + CELER_VALIDATE(opts.dp_max_alpha > 0, + << "invalid dp_max_alpha " << opts.dp_max_alpha); + CELER_ENSURE(opts); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/neutron/model/CascadeOptions.hh b/src/celeritas/neutron/model/CascadeOptions.hh new file mode 100644 index 0000000000..7972d10b55 --- /dev/null +++ b/src/celeritas/neutron/model/CascadeOptions.hh @@ -0,0 +1,128 @@ +//----------------------------------*-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/neutron/model/CascadeOptions.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +//#include "celeritas/Units.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Configuration options for the Bertini cascade model. + */ +struct CascadeOptions +{ + /// Top-level Configuration flags + + //! The flag for using the PreCompound model + bool use_precompound = false; + + //! The flag for using ABLA data + bool use_abla = false; + + //! The flag for using the three body momentum + bool use_three_body_mom = false; + + //! The flag for using the phase space + bool use_phase_space = false; + + //! The flag for applying coalescence + bool do_coalescence = true; + + //! The probability for the pion-nucleon absorption + real_type prob_pion_absorption = 0; + + /// Nuclear structure parameters + + //! The flag for using two parameters for the nuclear radius + bool use_two_params = false; + + //! The scale for the nuclear radius in [femtometer] + real_type radius_scale = 2.81967; + + //! The radius scale for small A (A < 4) + real_type radius_small = 8; + + //! The radius scale for alpha (A = 4) + real_type radius_alpha = 0.7; + + //! The trailing factor for the nuclear radius + real_type radius_trailing = 0; + + //! The scale for the Fermi radius in [MeV/c] + real_type fermi_scale = 1932; + + //! The scale for cross sections + real_type xsec_scale = 1; + + //! The scale for the quasi-deuteron cross section of gamma + real_type gamma_qd_scale = 1; + + /// Final-state clustering cuts + + //! The final state clustering cut (2 clusters) + real_type dp_max_doublet = 0.09; + + //! The final state clustering cut (3 clusters) + real_type dp_max_triplet = 0.108; + + //! The final state clustering cut (4 clusters) + real_type dp_max_alpha = 0.115; + + //! Whether all data are assigned and valid + explicit CELER_FUNCTION operator bool() const + { + // clang-format off + return (prob_pion_absorption >= 0) + && (radius_scale > 0) + && (radius_small > 0) + && (radius_alpha > 0) + && (radius_trailing >= 0) + && (fermi_scale > 0) + && (xsec_scale > 0) + && (gamma_qd_scale > 0) + && (dp_max_doublet > 0) + && (dp_max_triplet > 0) + && (dp_max_alpha > 0); + // clang-format on + } +}; + +//---------------------------------------------------------------------------// +//! Equality operator +constexpr bool operator==(CascadeOptions const& a, CascadeOptions const& b) +{ + // clang-format off + return a.use_precompound == b.use_precompound + && a.use_abla == b.use_abla + && a.use_three_body_mom == b.use_three_body_mom + && a.use_phase_space == b.use_phase_space + && a.do_coalescence == b.do_coalescence + && a.prob_pion_absorption == b.prob_pion_absorption + && a.use_two_params == b.use_two_params + && a.radius_scale == b.radius_scale + && a.radius_small == b.radius_small + && a.radius_alpha == b.radius_alpha + && a.radius_trailing == b.radius_trailing + && a.fermi_scale == b.fermi_scale + && a.xsec_scale == b.xsec_scale + && a.gamma_qd_scale == b.gamma_qd_scale + && a.dp_max_doublet == b.dp_max_doublet + && a.dp_max_triplet == b.dp_max_triplet + && a.dp_max_alpha == b.dp_max_alpha; + // clang-format on +} + +//---------------------------------------------------------------------------// +// Throw a runtime assertion if any of the input is invalid +void validate_input(CascadeOptions const&); + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/neutron/model/CascadeOptionsIO.json.cc b/src/celeritas/neutron/model/CascadeOptionsIO.json.cc new file mode 100644 index 0000000000..bc9f3995fa --- /dev/null +++ b/src/celeritas/neutron/model/CascadeOptionsIO.json.cc @@ -0,0 +1,76 @@ +//----------------------------------*-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/neutron/model/CascadeOptionsIO.json.cc +//---------------------------------------------------------------------------// +#include "CascadeOptionsIO.json.hh" + +#include +#include + +#include "corecel/io/JsonUtils.json.hh" + +#include "CascadeOptions.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Read options from JSON. + */ +void from_json(nlohmann::json const& j, CascadeOptions& opts) +{ +#define FDO_INPUT(NAME) CELER_JSON_LOAD_OPTION(j, opts, NAME) + + FDO_INPUT(use_precompound); + FDO_INPUT(use_abla); + FDO_INPUT(use_three_body_mom); + FDO_INPUT(use_phase_space); + FDO_INPUT(do_coalescence); + FDO_INPUT(prob_pion_absorption); + FDO_INPUT(use_two_params); + FDO_INPUT(radius_scale); + FDO_INPUT(radius_small); + FDO_INPUT(radius_alpha); + FDO_INPUT(radius_trailing); + FDO_INPUT(fermi_scale); + FDO_INPUT(xsec_scale); + FDO_INPUT(gamma_qd_scale); + FDO_INPUT(dp_max_doublet); + FDO_INPUT(dp_max_triplet); + FDO_INPUT(dp_max_alpha); + +#undef FDO_INPUT +} + +//---------------------------------------------------------------------------// +/*! + * Write options to JSON. + */ +void to_json(nlohmann::json& j, CascadeOptions const& opts) +{ + j = nlohmann::json{ + CELER_JSON_PAIR(opts, use_precompound), + CELER_JSON_PAIR(opts, use_abla), + CELER_JSON_PAIR(opts, use_three_body_mom), + CELER_JSON_PAIR(opts, use_phase_space), + CELER_JSON_PAIR(opts, do_coalescence), + CELER_JSON_PAIR(opts, prob_pion_absorption), + CELER_JSON_PAIR(opts, use_two_params), + CELER_JSON_PAIR(opts, radius_scale), + CELER_JSON_PAIR(opts, radius_small), + CELER_JSON_PAIR(opts, radius_alpha), + CELER_JSON_PAIR(opts, radius_trailing), + CELER_JSON_PAIR(opts, fermi_scale), + CELER_JSON_PAIR(opts, xsec_scale), + CELER_JSON_PAIR(opts, gamma_qd_scale), + CELER_JSON_PAIR(opts, dp_max_doublet), + CELER_JSON_PAIR(opts, dp_max_triplet), + CELER_JSON_PAIR(opts, dp_max_alpha), + }; +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/neutron/model/CascadeOptionsIO.json.hh b/src/celeritas/neutron/model/CascadeOptionsIO.json.hh new file mode 100644 index 0000000000..1cbe0bf43f --- /dev/null +++ b/src/celeritas/neutron/model/CascadeOptionsIO.json.hh @@ -0,0 +1,24 @@ +//----------------------------------*-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/neutron/model/CascadeOptionsIO.json.hh +//---------------------------------------------------------------------------// +#pragma once + +#include + +namespace celeritas +{ +struct CascadeOptions; +//---------------------------------------------------------------------------// + +// Read options from JSON +void from_json(nlohmann::json const& j, CascadeOptions& opts); + +// Write options to JSON +void to_json(nlohmann::json& j, CascadeOptions const& opts); + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/neutron/model/NeutronInelasticModel.cc b/src/celeritas/neutron/model/NeutronInelasticModel.cc index bbe2b1fc0a..54b9cde458 100644 --- a/src/celeritas/neutron/model/NeutronInelasticModel.cc +++ b/src/celeritas/neutron/model/NeutronInelasticModel.cc @@ -18,6 +18,10 @@ #include "celeritas/phys/PDGNumber.hh" #include "celeritas/phys/ParticleParams.hh" +#include "CascadeOptions.hh" + +#include "detail/NuclearZoneBuilder.hh" + namespace celeritas { //---------------------------------------------------------------------------// @@ -27,6 +31,7 @@ namespace celeritas NeutronInelasticModel::NeutronInelasticModel(ActionId id, ParticleParams const& particles, MaterialParams const& materials, + CascadeOptions const& options, ReadData load_data) { CELER_EXPECT(id); @@ -37,13 +42,16 @@ NeutronInelasticModel::NeutronInelasticModel(ActionId id, // Save IDs data.scalars.action_id = id; data.scalars.neutron_id = particles.find(pdg::neutron()); + data.scalars.proton_id = particles.find(pdg::proton()); - CELER_VALIDATE(data.scalars.neutron_id, - << "missing neutron particles (required for " + CELER_VALIDATE(data.scalars.neutron_id && data.scalars.proton_id, + << "missing neutron and/or proton particles (required for " << this->description() << ")"); // Save particle properties data.scalars.neutron_mass = particles.get(data.scalars.neutron_id).mass(); + data.scalars.proton_mass = particles.get(data.scalars.proton_id).mass(); + CELER_EXPECT(data.scalars); // Load neutron inelastic cross section data CollectionBuilder micro_xs{&data.micro_xs}; @@ -76,6 +84,16 @@ NeutronInelasticModel::NeutronInelasticModel(ActionId id, CELER_ASSERT(data.nucleon_xs.size() == num_channels); CELER_ASSERT(data.xs_params.size() == data.nucleon_xs.size()); + // Build (A, Z)-dependent nuclear zone data + detail::NuclearZoneBuilder build_nuclear_zones( + options, data.scalars, &data.nuclear_zones); + + for (auto iso_id : range(IsotopeId{materials.num_isotopes()})) + { + build_nuclear_zones(materials.get(iso_id)); + } + CELER_ASSERT(data.nuclear_zones.zones.size() == materials.num_isotopes()); + // Move to mirrored data, copying to device data_ = CollectionMirror{std::move(data)}; CELER_ENSURE(this->data_); diff --git a/src/celeritas/neutron/model/NeutronInelasticModel.hh b/src/celeritas/neutron/model/NeutronInelasticModel.hh index fd9ab43858..a703fb47c9 100644 --- a/src/celeritas/neutron/model/NeutronInelasticModel.hh +++ b/src/celeritas/neutron/model/NeutronInelasticModel.hh @@ -19,6 +19,7 @@ namespace celeritas { +struct CascadeOptions; struct ImportPhysicsVector; class MaterialParams; class ParticleParams; @@ -43,6 +44,7 @@ class NeutronInelasticModel final : public Model NeutronInelasticModel(ActionId id, ParticleParams const& particles, MaterialParams const& materials, + CascadeOptions const& options, ReadData load_data); // Particle types and energy ranges that this model applies to diff --git a/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh new file mode 100644 index 0000000000..96b835982d --- /dev/null +++ b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh @@ -0,0 +1,405 @@ +//----------------------------------*-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/neutron/model/detail/NuclearZoneBuilder.hh +//---------------------------------------------------------------------------// +#pragma once + +#include +#include + +#include "corecel/math/SoftEqual.hh" +#include "celeritas/Types.hh" +#include "celeritas/mat/IsotopeView.hh" +#include "celeritas/neutron/data/NeutronInelasticData.hh" +#include "celeritas/phys/AtomicNumber.hh" + +#include "../CascadeOptions.hh" + +namespace celeritas +{ +namespace detail +{ +//---------------------------------------------------------------------------// +/*! + * Construct NuclearZoneData for NeutronInelasticModel. + */ +class NuclearZoneBuilder +{ + public: + //!@{ + //! \name Type aliases + using AtomicMassNumber = AtomicNumber; + using MevMass = units::MevMass; + using Energy = units::MevEnergy; + using ComponentVec = std::vector; + using Data = HostVal; + //!@} + + public: + // Construct with cascade options and shared data + inline NuclearZoneBuilder(CascadeOptions const& options, + NeutronInelasticScalars const& scalars, + Data* data); + + // Construct nuclear zone data for a target (isotope) + inline void operator()(IsotopeView const& target); + + private: + // Calculate the number of nuclear zones + inline size_type num_nuclear_zones(AtomicMassNumber a); + + // Calculate the nuclear radius + inline real_type calc_nuclear_radius(AtomicMassNumber a); + + // Calculate components of nuclear zone properties + inline void + calc_zone_component(IsotopeView const& target, ComponentVec& components); + + // Integrate the Woods-Saxon potential in [rmin, rmax] + inline real_type integrate_woods_saxon(real_type rmin, + real_type rmax, + real_type radius) const; + + // Integrate the Gaussoan potential in [rmin, rmax] + inline real_type + integrate_gaussian(real_type rmin, real_type rmax, real_type radius) const; + + // Integrate a nuclear potential by adaptive quadrature + inline real_type integrate_potential(real_type rmin, + real_type delta_r, + real_type integral1, + real_type ws_shift = 0) const; + + private: + //// DATA //// + + // Cascade model configurations and nuclear structure parameters + CascadeOptions const& options_; + + real_type skin_depth_; + MevMass neutron_mass_; + MevMass proton_mass_; + + CollectionBuilder components_; + CollectionBuilder zones_; + + //// COMMON PROPERTIES //// + + static CELER_CONSTEXPR_FUNCTION real_type half() { return 0.5; } + static CELER_CONSTEXPR_FUNCTION real_type four_thirds_pi() + { + return 4 * constants::pi / real_type{3}; + } +}; + +//---------------------------------------------------------------------------// +// INLINE DEFINITIONS +//---------------------------------------------------------------------------// +/*! + * Construct with cascade options and data. + */ +NuclearZoneBuilder::NuclearZoneBuilder(CascadeOptions const& options, + NeutronInelasticScalars const& scalars, + Data* data) + : options_(options) + , skin_depth_(0.611207 * options.radius_scale) + , neutron_mass_(scalars.neutron_mass) + , proton_mass_(scalars.proton_mass) + , components_(&data->components) + , zones_(&data->zones) +{ + CELER_EXPECT(options_); + CELER_EXPECT(data); +} + +//---------------------------------------------------------------------------// +/*! + * Construct nuclear zone data for a single isotope. + */ +void NuclearZoneBuilder::operator()(IsotopeView const& target) +{ + AtomicNumber a = target.atomic_mass_number(); + size_type num_zones = this->num_nuclear_zones(a); + + std::vector comp(num_zones); + this->calc_zone_component(target, comp); + + NuclearZones nucl_zones; + nucl_zones.zones = components_.insert_back(comp.begin(), comp.end()); + zones_.push_back(nucl_zones); +} + +//---------------------------------------------------------------------------// +/*! + * Return the number of nuclear zones. + */ +size_type NuclearZoneBuilder::num_nuclear_zones(AtomicMassNumber a) +{ + if (a < AtomicMassNumber{5}) + { + return 1; + } + if (a < AtomicMassNumber{100}) + { + return 3; + } + return 6; +} + +//---------------------------------------------------------------------------// +/*! + * Calculate components of nuclear zone data: the nuclear zone radius, volume, + * density, Fermi momentum and potential function as in G4NucleiModel and as + * documented in section 24.2.3 of the Geant4 Physics Reference (release 11.2). + */ +void NuclearZoneBuilder::calc_zone_component(IsotopeView const& target, + ComponentVec& components) +{ + AtomicNumber a = target.atomic_mass_number(); + + // Calculate nuclear radius + real_type nuclear_radius = this->calc_nuclear_radius(a); + + real_type skin_ratio = nuclear_radius / skin_depth_; + real_type skin_decay = std::exp(-skin_ratio); + + // Temporary data for the zone-by-zone density function + size_type num_of_zones = components.size(); + std::vector integral(num_of_zones); + + // Fill the nuclear radius by each zone + auto amass = a.get(); + real_type ymin = (amass < 12) ? 0 : -skin_ratio; + + if (amass < 5) + { + // Light ions treated as simple balls + components[0].radius = nuclear_radius; + integral[0] = 1; + } + else if (amass < 12) + { + // Small nuclei have a three-zone Gaussian potential + real_type gauss_radius = std::sqrt( + ipow<2>(nuclear_radius) * (1 - 1 / static_cast(amass)) + + real_type{6.4}); + // Precompute y = sqrt(-log(alpha)) where alpha[3] = {0.7, 0.3, 0.01} + constexpr Real3 y = {0.597223, 1.09726, 2.14597}; + for (auto i : range(num_of_zones)) + { + components[i].radius = gauss_radius * y[i]; + integral[i] = this->integrate_gaussian(ymin, y[i], gauss_radius); + ymin = y[i]; + } + } + else if (amass < 100) + { + // Intermediate nuclei have a three-zone Woods-Saxon potential + constexpr Real3 alpha = {0.7, 0.3, 0.01}; + for (auto i : range(num_of_zones)) + { + real_type y = std::log((1 + skin_decay) / alpha[i] - 1); + components[i].radius = nuclear_radius + skin_depth_ * y; + integral[i] = this->integrate_woods_saxon(ymin, y, nuclear_radius); + ymin = y; + } + } + else + { + // Heavy nuclei have a six-zone Woods-Saxon potential + constexpr Array alpha = {0.9, 0.6, 0.4, 0.2, 0.1, 0.05}; + for (auto i : range(num_of_zones)) + { + real_type y = std::log((1 + skin_decay) / alpha[i] - 1); + components[i].radius = nuclear_radius + skin_depth_ * y; + integral[i] = this->integrate_woods_saxon(ymin, y, nuclear_radius); + ymin = y; + } + } + real_type total_integral + = std::accumulate(integral.begin(), integral.end(), real_type{0}); + + // Fill the nuclear volume by each zone + for (auto i : range(num_of_zones)) + { + components[i].volume = ipow<3>(components[i].radius); + if (i > 0) + { + components[i].volume -= ipow<3>(components[i - 1].radius); + } + components[i].volume *= this->four_thirds_pi(); + } + + // Fill the nuclear zone density, fermi momentum, and potential + int num_protons = target.atomic_number().get(); + int num_nucleons[] = {num_protons, a.get() - num_protons}; + real_type mass[] = {proton_mass_.value(), neutron_mass_.value()}; + real_type dm[] = {target.proton_loss_energy().value(), + target.neutron_loss_energy().value()}; + static_assert(std::size(mass) == ZoneComponent::NucleonArray::size()); + static_assert(std::size(dm) == std::size(mass)); + + for (auto i : range(num_of_zones)) + { + for (auto ptype : range(ZoneComponent::NucleonArray::size())) + { + components[i].density[ptype] + = static_cast(num_nucleons[ptype]) * integral[i] + / (total_integral * components[i].volume); + components[i].fermi_mom[ptype] + = options_.fermi_scale + * std::cbrt(components[i].density[ptype]); + components[i].potential[ptype] + = this->half() * ipow<2>(components[i].fermi_mom[ptype]) + / mass[ptype] + + dm[ptype]; + } + } +} + +//---------------------------------------------------------------------------// +/*! + * Calculate the nuclear radius (R) computed from the atomic mass number (A). + * + * For \f$ A > 4 \f$, the nuclear radius with two parameters takes the form, + * \f[ + R = [ 1.16 * A^{1/3} - 1.3456 / A^{1/3} ] \cdot R_{scale} + \f] + * where \f$ R_{scale} \f$ is a configurable parameter in [femtometer], while + * \f$ R = 1.2 A^{1/3} \cdot R_{scale} \f$ (default) with a single parameter. + */ +real_type NuclearZoneBuilder::calc_nuclear_radius(AtomicMassNumber a) +{ + // Nuclear radius computed from A + real_type nuclear_radius{0}; + auto amass = a.get(); + + if (amass > 4) + { + real_type cbrt_a = std::cbrt(static_cast(amass)); + real_type par_a = (options_.use_two_params ? 1.16 : 1.2); + real_type par_b = (options_.use_two_params ? -1.3456 : 0); + nuclear_radius = options_.radius_scale + * (par_a * cbrt_a + par_b / cbrt_a); + } + else + { + nuclear_radius = options_.radius_small + * (amass == 4 ? options_.radius_alpha : 1); + } + + return nuclear_radius; +} + +//---------------------------------------------------------------------------// +/*! + * Integrate the Woods-Saxon potential, \f$ V(r) \f$, + * + * \f[ + V(r) = frac{V_{o}}{1 + e^{\frac{r - R}{a}}} + \f] + * numerically over the volume from \f$ r_{min} \f$ to \f$ r_{rmax} \f$, where + * \f$ V_{o}, R, a\f$ are the potential well depth, nuclear radius, and + * surface thickness (skin depth), respectively. + */ +real_type +NuclearZoneBuilder::integrate_woods_saxon(real_type rmin, + real_type rmax, + real_type nuclear_radius) const +{ + real_type skin_ratio = nuclear_radius / skin_depth_; + + real_type ws_shift = 2 * skin_ratio; + real_type delta_r = rmax - rmin; + + // Average value of \f$ r^{2} V(r) dr \f$ at boundaries + real_type integral = half() * delta_r + * (rmin * (rmin + ws_shift) / (1 + std::exp(rmin)) + + rmax * (rmax + ws_shift) / (1 + std::exp(rmax))); + + real_type result = integrate_potential(rmin, delta_r, integral, ws_shift); + + return ipow<3>(skin_depth_) + * (result + + ipow<2>(skin_ratio) + * std::log((1 + std::exp(-rmin)) / (1 + std::exp(-rmax)))); +} + +//---------------------------------------------------------------------------// +/*! + * Integrate the Gaussian potential function from \f$ r_1 \f$ to \f$ r_2 \f$. + */ +real_type NuclearZoneBuilder::integrate_gaussian(real_type rmin, + real_type rmax, + real_type gauss_radius) const +{ + real_type delta_r = rmax - rmin; + real_type rmin_sq = ipow<2>(rmin); + real_type rmax_sq = ipow<2>(rmax); + + real_type integral + = this->half() * delta_r + * (rmin_sq * std::exp(-rmin_sq) + rmax_sq * std::exp(-rmax_sq)); + + real_type result = integrate_potential(rmin, delta_r, integral); + + return ipow<3>(gauss_radius) * result; +} + +//---------------------------------------------------------------------------// +/*! + * Integrate a nuclear potential numerically by adaptive quadrature. + */ +real_type NuclearZoneBuilder::integrate_potential(real_type rmin, + real_type delta_r, + real_type integral, + real_type ws_shift) const +{ + constexpr int max_trials = 1000; + constexpr real_type epsilon = 1e-3; + + int depth = 1; + real_type interval = delta_r; + real_type result = 0; + + bool succeeded = false; + int remaining_trials = max_trials; + SoftEqual const soft_eq{epsilon}; + + do + { + delta_r *= this->half(); + + real_type r = rmin - delta_r; + real_type fi = 0; + + for (int i = 0; i < depth; ++i) + { + r += interval; + fi += (ws_shift > 0) ? r * (r + ws_shift) / (1 + std::exp(r)) + : ipow<2>(r) * std::exp(-ipow<2>(r)); + } + + result = this->half() * integral + fi * delta_r; + + if (soft_eq(1, integral / result)) + { + succeeded = true; + } + else + { + depth *= 2; + interval = delta_r; + integral = result; + } + } while (!succeeded && --remaining_trials > 0); + + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/src/celeritas/neutron/process/NeutronInelasticProcess.cc b/src/celeritas/neutron/process/NeutronInelasticProcess.cc index f687812202..09e1ee3851 100644 --- a/src/celeritas/neutron/process/NeutronInelasticProcess.cc +++ b/src/celeritas/neutron/process/NeutronInelasticProcess.cc @@ -9,6 +9,7 @@ #include "corecel/Assert.hh" #include "celeritas/grid/ValueGridBuilder.hh" +#include "celeritas/neutron/model/CascadeOptions.hh" #include "celeritas/neutron/model/NeutronInelasticModel.hh" #include "celeritas/phys/PDGNumber.hh" @@ -38,8 +39,9 @@ NeutronInelasticProcess::NeutronInelasticProcess(SPConstParticles particles, */ auto NeutronInelasticProcess::build_models(ActionIdIter id) const -> VecModel { + CascadeOptions options; // TODO: options as an argument of this process return {std::make_shared( - *id++, *particles_, *materials_, load_data_)}; + *id++, *particles_, *materials_, options, load_data_)}; } //---------------------------------------------------------------------------// diff --git a/test/celeritas/em/CoulombScattering.test.cc b/test/celeritas/em/CoulombScattering.test.cc index db30d37f70..2651c66e5f 100644 --- a/test/celeritas/em/CoulombScattering.test.cc +++ b/test/celeritas/em/CoulombScattering.test.cc @@ -60,11 +60,15 @@ class CoulombScatteringTest : public InteractorHostTestBase mat_inp.isotopes = {{AtomicNumber{29}, AtomicNumber{63}, MevEnergy{551.384}, + MevEnergy{6.122}, + MevEnergy{10.864}, MevMass{58618.5}, "63Cu"}, {AtomicNumber{29}, AtomicNumber{65}, MevEnergy{569.211}, + MevEnergy{7.454}, + MevEnergy{9.911}, MevMass{60479.8}, "65Cu"}}; mat_inp.elements = {{AtomicNumber{29}, diff --git a/test/celeritas/mat/Material.test.cc b/test/celeritas/mat/Material.test.cc index bcef56e9e8..e0ea9ef5cd 100644 --- a/test/celeritas/mat/Material.test.cc +++ b/test/celeritas/mat/Material.test.cc @@ -320,7 +320,7 @@ TEST_F(MaterialTest, output) if (CELERITAS_USE_JSON && CELERITAS_UNITS == CELERITAS_UNITS_CGS) { EXPECT_JSON_EQ( - R"json({"_category":"internal","_label":"material","_units":{"atomic_mass":"amu","mean_excitation_energy":"MeV","binding_energy":"MeV","nuclear_mass":"MeV/c^2"},"elements":{"atomic_mass":[1.008,26.9815385,22.98976928,126.90447],"atomic_number":[1,13,11,53],"coulomb_correction":[6.400821803338426e-05,0.010734632775699565,0.00770256745342534,0.15954439947436763],"isotope_fractions":[[0.9,0.1],[0.7,0.3],[1.0],[0.05,0.15,0.8]],"isotope_ids":[[0,1],[2,3],[4],[5,6,7]],"label":["H","Al","Na","I"],"mass_radiation_coeff":[0.0158611264432063,0.04164723292591279,0.03605392839455309,0.11791841505608874]},"isotopes":{"atomic_mass_number":[1,2,27,28,23,125,126,127],"atomic_number":[1,1,13,13,11,53,53,53],"binding_energy":[0.0,2.22457,224.952,232.677,186.564,1056.29,1063.43,1072.58],"label":["1H","2H","27Al","28Al","23Na","125I","126I","127I"],"nuclear_mass":[938.272,1875.61,25126.5,26058.3,21409.2,116321.0,117253.0,118184.0]},"materials":{"density":[3.6700020622594716,0.0,0.00017976000000000003,0.00017943386624303615],"electron_density":[9.4365282069664e+23,0.0,1.073948435904467e+20,1.072e+20],"element_frac":[[0.5,0.5],[],[1.0],[1.0]],"element_id":[[2,3],[],[0],[0]],"label":["NaI","hard vacuum","H2@1","H2@2"],"matter_state":["solid","unspecified","gas","gas"],"mean_excitation_energy":[0.00040000760709482647,0.0,1.9199999999999986e-05,1.9199999999999986e-05],"number_density":[2.948915064677e+22,0.0,1.073948435904467e+20,1.072e+20],"radiation_length":[3.5393292693170424,null,350729.99844063615,351367.4750467326],"temperature":[293.0,0.0,100.0,110.0],"zeff":[32.0,0.0,1.0,1.0]}})json", + R"json({"_category":"internal","_label":"material","_units":{"atomic_mass":"amu","binding_energy":"MeV","mean_excitation_energy":"MeV","neutron_loss_energy":"MeV","nuclear_mass":"MeV/c^2","proton_loss_energy":"MeV"},"elements":{"atomic_mass":[1.008,26.9815385,22.98976928,126.90447],"atomic_number":[1,13,11,53],"coulomb_correction":[6.400821803338426e-05,0.010734632775699565,0.00770256745342534,0.15954439947436763],"isotope_fractions":[[0.9,0.1],[0.7,0.3],[1.0],[0.05,0.15,0.8]],"isotope_ids":[[0,1],[2,3],[4],[5,6,7]],"label":["H","Al","Na","I"],"mass_radiation_coeff":[0.0158611264432063,0.04164723292591279,0.03605392839455309,0.11791841505608874]},"isotopes":{"atomic_mass_number":[1,2,27,28,23,125,126,127],"atomic_number":[1,1,13,13,11,53,53,53],"binding_energy":[0.0,2.22457,224.952,232.677,186.564,1056.29,1063.43,1072.58],"label":["1H","2H","27Al","28Al","23Na","125I","126I","127I"],"neutron_loss_energy":[0.0,2.22457,13.058,7.725,12.42,9.543,7.145,9.144],"nuclear_mass":[938.272,1875.61,25126.5,26058.3,21409.2,116321.0,117253.0,118184.0],"proton_loss_energy":[0.0,2.22457,8.271,9.553,8.794,5.601,6.177,6.208]},"materials":{"density":[3.6700020622594716,0.0,0.00017976000000000003,0.00017943386624303615],"electron_density":[9.4365282069664e+23,0.0,1.073948435904467e+20,1.072e+20],"element_frac":[[0.5,0.5],[],[1.0],[1.0]],"element_id":[[2,3],[],[0],[0]],"label":["NaI","hard vacuum","H2@1","H2@2"],"matter_state":["solid","unspecified","gas","gas"],"mean_excitation_energy":[0.00040000760709482647,0.0,1.9199999999999986e-05,1.9199999999999986e-05],"number_density":[2.948915064677e+22,0.0,1.073948435904467e+20,1.072e+20],"radiation_length":[3.5393292693170424,null,350729.99844063615,351367.4750467326],"temperature":[293.0,0.0,100.0,110.0],"zeff":[32.0,0.0,1.0,1.0]}})json", to_string(out)); } } diff --git a/test/celeritas/mat/MaterialTestBase.hh b/test/celeritas/mat/MaterialTestBase.hh index b0598bc4d9..1b2b6774c8 100644 --- a/test/celeritas/mat/MaterialTestBase.hh +++ b/test/celeritas/mat/MaterialTestBase.hh @@ -31,44 +31,60 @@ class MaterialTestBase {AtomicNumber{1}, AtomicNumber{1}, MevEnergy{0}, + MevEnergy{0}, + MevEnergy{0}, MevMass{938.272}, "1H"}, {AtomicNumber{1}, AtomicNumber{2}, MevEnergy{2.22457}, + MevEnergy{2.22457}, + MevEnergy{2.22457}, MevMass{1875.61}, "2H"}, // Al {AtomicNumber{13}, AtomicNumber{27}, MevEnergy{224.952}, + MevEnergy{8.271}, + MevEnergy{13.058}, MevMass{25126.5}, "27Al"}, {AtomicNumber{13}, AtomicNumber{28}, MevEnergy{232.677}, + MevEnergy{9.553}, + MevEnergy{7.725}, MevMass{26058.3}, "28Al"}, // Na {AtomicNumber{11}, AtomicNumber{23}, MevEnergy{186.564}, + MevEnergy{8.794}, + MevEnergy{12.420}, MevMass{21409.2}, "23Na"}, // I {AtomicNumber{53}, AtomicNumber{125}, MevEnergy{1056.29}, + MevEnergy{5.601}, + MevEnergy{9.543}, MevMass{116321}, "125I"}, {AtomicNumber{53}, AtomicNumber{126}, MevEnergy{1063.43}, + MevEnergy{6.177}, + MevEnergy{7.145}, MevMass{117253}, "126I"}, {AtomicNumber{53}, AtomicNumber{127}, MevEnergy{1072.58}, + MevEnergy{6.208}, + MevEnergy{9.144}, MevMass{118184}, "127I"}}; diff --git a/test/celeritas/neutron/NeutronElastic.test.cc b/test/celeritas/neutron/NeutronElastic.test.cc index 877039300a..3ef3763a4d 100644 --- a/test/celeritas/neutron/NeutronElastic.test.cc +++ b/test/celeritas/neutron/NeutronElastic.test.cc @@ -249,15 +249,15 @@ TEST_F(NeutronElasticTest, extended) 999.98224958135802, 9999.9796839871542}; std::vector const expected_angle = {0.73642517015232023, - -0.93575721388581834, + -0.93575721388581845, 0.39720848171626949, - 0.63289309823793904, + 0.63289309823793916, -0.98646693371059246, -0.15892409193659782, 0.99930375891282708, 0.97355882843189279, 0.99963860089317569, - 0.99998997488585095}; + 0.99998997488584618}; real_type energy = 1e-5; real_type const factor = 1e+1; @@ -298,8 +298,8 @@ TEST_F(NeutronElasticTest, stress_test) 99.915949160022208, 9999.9519447032781}; std::vector cos_theta; - real_type const expected_angle[] = {-0.0062993051693808555, - 0.0036712145032219813, + real_type const expected_angle[] = {-0.0062993051693811574, + 0.0036712145032221713, -0.29680452456419526, 0.97426025128623828, 0.9999755334707946}; diff --git a/test/celeritas/neutron/NeutronInelastic.test.cc b/test/celeritas/neutron/NeutronInelastic.test.cc index 8005ca7da9..ad6a2fe785 100644 --- a/test/celeritas/neutron/NeutronInelastic.test.cc +++ b/test/celeritas/neutron/NeutronInelastic.test.cc @@ -16,6 +16,7 @@ #include "celeritas/mat/MaterialTrackView.hh" #include "celeritas/neutron/NeutronTestBase.hh" #include "celeritas/neutron/interactor/NeutronInelasticInteractor.hh" +#include "celeritas/neutron/model/CascadeOptions.hh" #include "celeritas/neutron/model/NeutronInelasticModel.hh" #include "celeritas/neutron/xs/NeutronInelasticMicroXsCalculator.hh" #include "celeritas/neutron/xs/NucleonNucleonXsCalculator.hh" @@ -47,10 +48,15 @@ class NeutronInelasticTest : public NeutronTestBase this->set_inc_particle(pdg::neutron(), MevEnergy{100}); this->set_inc_direction({0, 0, 1}); - // Set up the default material + // Build the model with the default material this->set_material("HeCu"); - model_ = std::make_shared( - ActionId{0}, particles, *this->material_params(), read_el_data); + CascadeOptions options; + model_ + = std::make_shared(ActionId{0}, + particles, + *this->material_params(), + options, + read_el_data); } protected: @@ -174,6 +180,202 @@ TEST_F(NeutronInelasticTest, nucleon_xs) EXPECT_VEC_SOFT_EQ(expected_xs, xs); } +TEST_F(NeutronInelasticTest, model_data) +{ + // Test neutron inelastic interactions + NeutronInelasticRef shared = model_->host_ref(); + + // Set the target to (light) isotope He3 + IsotopeId iso_id{0}; + + // Check the size of the number of nuclear zones + NuclearZones he3_nuclear_zones = shared.nuclear_zones.zones[iso_id]; + EXPECT_EQ(he3_nuclear_zones.zones.size(), 1); + + // Check zone data + for (auto sid : he3_nuclear_zones.zones) + { + ZoneComponent components = shared.nuclear_zones.components[sid]; + EXPECT_SOFT_EQ(8, components.radius); + EXPECT_SOFT_EQ(2144.6605848506319, components.volume); + // proton + EXPECT_SOFT_EQ(0.00093254849467907434, components.density[0]); + EXPECT_SOFT_EQ(188.75462299392046, components.fermi_mom[0]); + EXPECT_SOFT_EQ(24.476129399543886, components.potential[0]); + // neutron + EXPECT_SOFT_EQ(0.00046627424733953717, components.density[1]); + EXPECT_SOFT_EQ(149.81464355220513, components.fermi_mom[1]); + EXPECT_SOFT_EQ(55.944047271729367, components.potential[1]); + } + + // Set the target to (heavy) isotope Cu63 + IsotopeId iso_cu63{2}; + + // Check the size of the number of nuclear zones + NuclearZones cu63_nuclear_zones = shared.nuclear_zones.zones[iso_cu63]; + EXPECT_EQ(cu63_nuclear_zones.zones.size(), 3); + + // Check zone data + std::vector radii; + std::vector volumes; + std::vector densities; + std::vector fermi_moms; + std::vector potentials; + + for (auto sid : cu63_nuclear_zones.zones) + { + ZoneComponent components = shared.nuclear_zones.components[sid]; + radii.push_back(components.radius); + volumes.push_back(components.volume); + for (auto nucleon_index : range(2)) + { + densities.push_back(components.density[nucleon_index]); + fermi_moms.push_back(components.fermi_mom[nucleon_index]); + potentials.push_back(components.potential[nucleon_index]); + } + } + + real_type const expected_cu_radii[] + = {12.0056427171327, 14.924785000235, 21.383497297248}; + + real_type const expected_cu_volumes[] + = {7248.44509638664, 6677.12103595803, 27031.1207141316}; + + // clang-format off + real_type const expected_cu_densities[] + = {0.00218857360950474, 0.00256591388700555, 0.0011959208904212, + 0.00140211414739038, 0.000190555762441557, 0.000223410204241825}; + + real_type const expected_cu_fermi_moms[] + = {250.838488281962, 264.497241970299, 205.072742538403, + 216.239442265022, 111.176880545669, 117.23072673814}; + + real_type const expected_cu_potentials[] + = {39.6516941248423, 48.0933349774238, 28.5327876764636, + 35.7475768798981, 12.7087352945686, 18.1775106385426}; + // clang-format on + + EXPECT_VEC_SOFT_EQ(expected_cu_radii, radii); + EXPECT_VEC_SOFT_EQ(expected_cu_volumes, volumes); + EXPECT_VEC_SOFT_EQ(expected_cu_densities, densities); + EXPECT_VEC_SOFT_EQ(expected_cu_fermi_moms, fermi_moms); + EXPECT_VEC_SOFT_EQ(expected_cu_potentials, potentials); + + // Set the target to (very heavy, A > 100) isotope Pb208 + IsotopeId iso_pb208{7}; + + // Check the size of the number of nuclear zones + NuclearZones pb208_nuclear_zones = shared.nuclear_zones.zones[iso_pb208]; + EXPECT_EQ(pb208_nuclear_zones.zones.size(), 6); + + // Check zone data + radii.clear(); + volumes.clear(); + densities.clear(); + fermi_moms.clear(); + potentials.clear(); + + for (auto sid : pb208_nuclear_zones.zones) + { + ZoneComponent components = shared.nuclear_zones.components[sid]; + radii.push_back(components.radius); + volumes.push_back(components.volume); + for (auto nucleon_index : range(2)) + { + densities.push_back(components.density[nucleon_index]); + fermi_moms.push_back(components.fermi_mom[nucleon_index]); + potentials.push_back(components.potential[nucleon_index]); + } + } + + // clang-format off + real_type const expected_pb_radii[] + = {16.2612786504728, 19.3490859199713, 20.7466319700195, + 22.4369887370502, 23.834545403954, 25.122295335588}; + + real_type const expected_pb_volumes[] + = {18011.6162284334, 12332.183871369, 7061.35140601041, + 9908.04814765151, 9403.27500414658, 9698.58388620079}; + + real_type const expected_pb_densities[] + = {0.00225402783765291, 0.00346350618956423, 0.00177111372940623, + 0.00272146743786811, 0.00115221167124763, 0.0017704715924049, + 0.000673003348019808, 0.00103412709573775, 0.000333946126232526, + 0.000513136730552418, 0.000166530199539091, 0.000255887867584456}; + + real_type const expected_pb_fermi_moms[] + = {253.314595582345, 292.311416929582, 233.752324669988, + 269.737608596091, 202.5432986892, 233.72407141927, + 169.308799502358, 195.373247117507, 134.039436095239, + 154.674298965551, 106.293030558995, 122.656439519448}; + + real_type const expected_pb_potentials[] + = {42.1989261226909, 52.8390035394296, 37.1214353127717, + 46.0871655105693, 29.8653511196685, 36.4383237834904, + 23.279671229538, 27.6809580702491, 17.5782866566401, + 20.0994918268739, 14.0247531445444, 15.3741494083458}; + // clang-format on + + EXPECT_VEC_SOFT_EQ(expected_pb_radii, radii); + EXPECT_VEC_SOFT_EQ(expected_pb_volumes, volumes); + EXPECT_VEC_SOFT_EQ(expected_pb_densities, densities); + EXPECT_VEC_SOFT_EQ(expected_pb_fermi_moms, fermi_moms); + EXPECT_VEC_SOFT_EQ(expected_pb_potentials, potentials); + + // Set the target to isotope B11 and validate the Gaussian potential + IsotopeId iso_b11{9}; + + // Check the size of the number of nuclear zones + NuclearZones b11_nuclear_zones = shared.nuclear_zones.zones[iso_b11]; + EXPECT_EQ(b11_nuclear_zones.zones.size(), 3); + + // Clear zone data + radii.clear(); + volumes.clear(); + densities.clear(); + fermi_moms.clear(); + potentials.clear(); + + for (auto sid : b11_nuclear_zones.zones) + { + ZoneComponent components = shared.nuclear_zones.components[sid]; + radii.push_back(components.radius); + volumes.push_back(components.volume); + for (auto nucleon_index : range(2)) + { + densities.push_back(components.density[nucleon_index]); + fermi_moms.push_back(components.fermi_mom[nucleon_index]); + potentials.push_back(components.potential[nucleon_index]); + } + } + + real_type const expected_b_radii[] + = {4.54355900675881, 8.34774540792328, 16.3261316488719}; + + real_type const expected_b_volumes[] + = {392.895565424687, 2043.77151146382, 15791.3107610447}; + + // clang-format off + real_type const expected_b_densities[] + = {0.00169984536247725, 0.0020398144349727, 0.000949798766722706, + 0.00113975852006725, 0.00015141027051547, 0.000181692324618564}; + + real_type const expected_b_fermi_moms[] + = {230.573961004183, 245.021395491472, 189.911379816935, + 201.810955147759, 102.973522012059, 109.425695565028}; + + real_type const expected_b_potentials[] + = {39.5599907769562, 43.4025388663509, 30.4485502393036, + 33.1276701038221, 16.8795715233179, 17.8260857964719}; + // clang-format on + + EXPECT_VEC_SOFT_EQ(expected_b_radii, radii); + EXPECT_VEC_SOFT_EQ(expected_b_volumes, volumes); + EXPECT_VEC_SOFT_EQ(expected_b_densities, densities); + EXPECT_VEC_SOFT_EQ(expected_b_fermi_moms, fermi_moms); + EXPECT_VEC_SOFT_EQ(expected_b_potentials, potentials); +} + //---------------------------------------------------------------------------// } // namespace test } // namespace celeritas diff --git a/test/celeritas/neutron/NeutronTestBase.cc b/test/celeritas/neutron/NeutronTestBase.cc index 68eddbb5c2..7816942c7f 100644 --- a/test/celeritas/neutron/NeutronTestBase.cc +++ b/test/celeritas/neutron/NeutronTestBase.cc @@ -27,11 +27,13 @@ NeutronTestBase::NeutronTestBase() using namespace units; constexpr auto zero = zero_quantity(); - constexpr MevMass neutronmass{939.5654133}; + constexpr MevMass neutron_mass{939.5654133}; + constexpr MevMass proton_mass{938.272013}; // Setup default particle params ParticleParams::Input par_inp = { - {"neutron", pdg::neutron(), neutronmass, zero, stable_decay_constant}}; + {"neutron", pdg::neutron(), neutron_mass, zero, stable_decay_constant}, + {"proton", pdg::proton(), proton_mass, zero, stable_decay_constant}}; this->set_particle_params(std::move(par_inp)); // Setup default material params @@ -41,23 +43,73 @@ NeutronTestBase::NeutronTestBase() mat_inp.isotopes = {{AtomicNumber{2}, AtomicNumber{3}, units::MevEnergy{7.71804}, + units::MevEnergy{5.49}, + units::MevEnergy{44}, units::MevMass{3016.0}, "3He"}, {AtomicNumber{2}, AtomicNumber{4}, units::MevEnergy{28.2957}, + units::MevEnergy{19.814}, + units::MevEnergy{20.578}, units::MevMass{4002.6}, "4He"}, {AtomicNumber{29}, AtomicNumber{63}, units::MevEnergy{551.384}, + units::MevEnergy{6.122}, + units::MevEnergy{10.864}, units::MevMass{58618.5}, "63Cu"}, {AtomicNumber{29}, AtomicNumber{65}, units::MevEnergy{569.211}, + units::MevEnergy{7.454}, + units::MevEnergy{9.911}, units::MevMass{60479.8}, - "65Cu"}}; + "65Cu"}, + {AtomicNumber{82}, + AtomicNumber{202}, + units::MevEnergy{1607.5056}, + units::MevEnergy{6.637}, + units::MevEnergy{8.395}, + units::MevMass{189958.3388}, + "202Pb"}, + {AtomicNumber{82}, + AtomicNumber{206}, + units::MevEnergy{1622.3240}, + units::MevEnergy{7.254}, + units::MevEnergy{8.087}, + units::MevMass{191822.6512}, + "206Pb"}, + {AtomicNumber{82}, + AtomicNumber{207}, + units::MevEnergy{1629.0618}, + units::MevEnergy{7.488}, + units::MevEnergy{6.738}, + units::MevMass{192755.4787}, + "207Pb"}, + {AtomicNumber{82}, + AtomicNumber{208}, + units::MevEnergy{1636.4296}, + units::MevEnergy{8.004}, + units::MevEnergy{7.368}, + units::MevMass{193687.6762}, + "208Pb"}, + {AtomicNumber{5}, + AtomicNumber{10}, + units::MevEnergy{64.751}, + units::MevEnergy{6.587}, + units::MevEnergy{8.437}, + units::MevMass{9324.437}, + "10B"}, + {AtomicNumber{5}, + AtomicNumber{11}, + units::MevEnergy{76.205}, + units::MevEnergy{11.229}, + units::MevEnergy{11.454}, + units::MevMass{10252.548}, + "11B"}}; // Elements mat_inp.elements = {{AtomicNumber{2},