From ff78c9bdaaa33832be13e29cd6ed2f752bdcb7e3 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Wed, 12 Jun 2024 07:56:34 -0500 Subject: [PATCH 01/14] Import nucleon loss energies in IsotopeView --- app/celer-dump-data.cc | 6 ++- src/celeritas/ext/GeantImporter.cc | 10 +++++ src/celeritas/io/ImportElement.hh | 4 +- src/celeritas/mat/IsotopeView.hh | 24 +++++++++++ src/celeritas/mat/MaterialData.hh | 4 +- src/celeritas/mat/MaterialParams.cc | 8 ++++ src/celeritas/mat/MaterialParams.hh | 5 ++- src/celeritas/mat/MaterialParamsOutput.cc | 11 ++++++ test/celeritas/em/CoulombScattering.test.cc | 4 ++ test/celeritas/mat/Material.test.cc | 2 +- test/celeritas/mat/MaterialTestBase.hh | 16 ++++++++ test/celeritas/neutron/NeutronTestBase.cc | 44 +++++++++++++++++++-- 12 files changed, 129 insertions(+), 9 deletions(-) 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/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index 4bcb2211d5..2c2850833d 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -425,6 +425,16 @@ std::vector import_isotopes() isotope.atomic_mass_number = g4isotope.GetN(); isotope.binding_energy = G4NucleiProperties::GetBindingEnergy( 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); 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..08204c2d80 100644 --- a/src/celeritas/io/ImportElement.hh +++ b/src/celeritas/io/ImportElement.hh @@ -24,7 +24,9 @@ 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 binding_energy; //!< Nuclear binding energy (BE) [MeV] + double proton_loss_energy; //!< BE(A, Z) - BE(A-1, Z-1) [MeV] + double neutron_loss_energy; //!< BE(A, Z) - BE(A-1, Z) [MeV] double nuclear_mass; //!< 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/test/celeritas/em/CoulombScattering.test.cc b/test/celeritas/em/CoulombScattering.test.cc index f4d2a106e1..ce9face8dd 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/NeutronTestBase.cc b/test/celeritas/neutron/NeutronTestBase.cc index 68eddbb5c2..289a1b5271 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,59 @@ 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"}}; // Elements mat_inp.elements = {{AtomicNumber{2}, From 8004439d73811b552d37002a04220a8c505b6006 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Wed, 12 Jun 2024 07:57:18 -0500 Subject: [PATCH 02/14] Add options for the Bertini cascade model --- src/celeritas/CMakeLists.txt | 2 + src/celeritas/neutron/model/CascadeOptions.cc | 47 +++++++ src/celeritas/neutron/model/CascadeOptions.hh | 128 ++++++++++++++++++ .../neutron/model/CascadeOptionsIO.json.cc | 81 +++++++++++ .../neutron/model/CascadeOptionsIO.json.hh | 24 ++++ .../process/NeutronInelasticProcess.cc | 4 +- 6 files changed, 285 insertions(+), 1 deletion(-) create mode 100644 src/celeritas/neutron/model/CascadeOptions.cc create mode 100644 src/celeritas/neutron/model/CascadeOptions.hh create mode 100644 src/celeritas/neutron/model/CascadeOptionsIO.json.cc create mode 100644 src/celeritas/neutron/model/CascadeOptionsIO.json.hh 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/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..5b2f1b4220 --- /dev/null +++ b/src/celeritas/neutron/model/CascadeOptionsIO.json.cc @@ -0,0 +1,81 @@ +//----------------------------------*-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) \ + do \ + { \ + if (j.contains(#NAME)) \ + j.at(#NAME).get_to(opts.NAME); \ + } while (0) + + 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/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_)}; } //---------------------------------------------------------------------------// From 67c7137bf2fa10ab9c9ad1198538b4c8cc7826f7 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Wed, 12 Jun 2024 07:58:26 -0500 Subject: [PATCH 03/14] Adjust test results due to updates for nuclear mass --- test/celeritas/neutron/NeutronElastic.test.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) 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}; From 5896714cca7b453ad140f54d8f3c9d2cd7f4a3e3 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Wed, 12 Jun 2024 07:59:24 -0500 Subject: [PATCH 04/14] Add nuclear zone data for the Bertini cascade model --- .../neutron/data/NeutronInelasticData.hh | 79 +++- .../neutron/model/NeutronInelasticModel.cc | 22 +- .../neutron/model/NeutronInelasticModel.hh | 2 + .../model/detail/NuclearZoneBuilder.hh | 410 ++++++++++++++++++ .../neutron/NeutronInelastic.test.cc | 155 ++++++- 5 files changed, 661 insertions(+), 7 deletions(-) create mode 100644 src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh diff --git a/src/celeritas/neutron/data/NeutronInelasticData.hh b/src/celeritas/neutron/data/NeutronInelasticData.hh index 45aedb11a8..6be72169bd 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,73 @@ 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 +{ + size_type num_zones{}; + ItemRange zones; + + //! Whether all data are assigned and valid + explicit CELER_FUNCTION operator bool() const + { + return num_zones > 0 && !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 +161,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 +181,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/NeutronInelasticModel.cc b/src/celeritas/neutron/model/NeutronInelasticModel.cc index bbe2b1fc0a..1ba144a954 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 zone_builder( + options, data.scalars, &data.nuclear_zones); + + for (auto iso_id : range(IsotopeId{materials.num_isotopes()})) + { + zone_builder(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..ad83be0306 --- /dev/null +++ b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh @@ -0,0 +1,410 @@ +//----------------------------------*-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 "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 + explicit 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. + */ +CELER_FUNCTION +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.num_zones = num_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) +{ + return (a.get() < 5) ? 1 : (a.get() < 100) ? 3 : 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). + */ +CELER_FUNCTION +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 ur(num_of_zones + 1); + std::vector integral(num_of_zones); + real_type total_integral{0}; + + // Fill the nuclear radius by each zone + auto amass = a.get(); + ur[0] = (amass < 12) ? 0 : -skin_ratio; + + if (amass < 5) + { + // Light ions treated as simple balls + components[0].radius = nuclear_radius; + integral[0] = 1; + total_integral = integral[0]; + } + else if (amass < 100) + { + constexpr Real3 alpha = {0.7, 0.3, 0.01}; + 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}); + for (auto i : range(num_of_zones)) + { + real_type y = std::sqrt(-std::log(alpha[i])); + components[i].radius = gauss_radius * y; + ur[i + 1] = y; + integral[i] + = this->integrate_gaussian(ur[i], ur[i + 1], gauss_radius); + total_integral += integral[i]; + } + } + else + { + // Intermediate nuclei have a three-zone Woods-Saxon potential + 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; + ur[i + 1] = y; + integral[i] = this->integrate_woods_saxon( + ur[i], ur[i + 1], nuclear_radius); + total_integral += integral[i]; + } + } + } + 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; + ur[i + 1] = y; + integral[i] = this->integrate_woods_saxon( + ur[i], ur[i + 1], nuclear_radius); + total_integral += integral[i]; + } + } + + ur.clear(); + + // 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 const dim = ZoneComponent::NucleonArray::size(); + int num_protons = target.atomic_number().get(); + int num_nucleons[dim] = {num_protons, a.get() - num_protons}; + real_type mass[dim] = {proton_mass_.value(), neutron_mass_.value()}; + real_type dm[dim] = {target.proton_loss_energy().value(), + target.neutron_loss_energy().value()}; + + for (auto i : range(num_of_zones)) + { + for (auto ptype : range(dim)) + { + 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]; + } + } + + integral.clear(); +} + +//---------------------------------------------------------------------------// +/*! + * 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(false, 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 = 1e+3; + 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; + 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 (std::fabs(1 - integral / result) < epsilon) + { + succeeded = true; + } + else + { + depth *= 2; + interval = delta_r; + integral = result; + } + } while (!succeeded && --remaining_trials > 0); + + return result; +} + +//---------------------------------------------------------------------------// +} // namespace detail +} // namespace celeritas diff --git a/test/celeritas/neutron/NeutronInelastic.test.cc b/test/celeritas/neutron/NeutronInelastic.test.cc index 8005ca7da9..0b10dbce7a 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,149 @@ 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.num_zones, 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.num_zones, 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.num_zones, 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); +} + //---------------------------------------------------------------------------// } // namespace test } // namespace celeritas From b5d56652ee044ba80e6f72a1bbe8074d435976bd Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Wed, 12 Jun 2024 09:16:21 -0500 Subject: [PATCH 05/14] Fix an error from window build --- src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh index ad83be0306..804e87cc43 100644 --- a/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh +++ b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh @@ -365,7 +365,7 @@ real_type NuclearZoneBuilder::integrate_potential(real_type rmin, real_type integral, real_type ws_shift) const { - constexpr int max_trials = 1e+3; + constexpr int max_trials = 1000; constexpr real_type epsilon = 1e-3; int depth = 1; From c8bda1a5a7121afe7971fc22fa5f17adabf0a79c Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Wed, 12 Jun 2024 11:11:30 -0500 Subject: [PATCH 06/14] Guard for A and Z when calculating the binding energy difference for lossing a nucleon --- src/celeritas/ext/GeantImporter.cc | 33 +++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 10 deletions(-) diff --git a/src/celeritas/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index 2c2850833d..b8ea92e104 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -425,16 +425,29 @@ std::vector import_isotopes() isotope.atomic_mass_number = g4isotope.GetN(); isotope.binding_energy = G4NucleiProperties::GetBindingEnergy( 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); + + // Binding energy difference for lossing 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); } From 84d053dca36f85c19ffbe868394ba678b89ada80 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 06:23:31 -0500 Subject: [PATCH 07/14] Fix a typo --- src/celeritas/ext/GeantImporter.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/celeritas/ext/GeantImporter.cc b/src/celeritas/ext/GeantImporter.cc index b8ea92e104..905aa35e73 100644 --- a/src/celeritas/ext/GeantImporter.cc +++ b/src/celeritas/ext/GeantImporter.cc @@ -426,7 +426,7 @@ std::vector import_isotopes() isotope.binding_energy = G4NucleiProperties::GetBindingEnergy( isotope.atomic_mass_number, isotope.atomic_number); - // Binding energy difference for lossing a nucleon + // Binding energy difference for losing a nucleon if (isotope.atomic_mass_number > 1 && isotope.atomic_number > 1 && isotope.atomic_mass_number >= isotope.atomic_number) { From 3dc753f074c78804106444c90f8295c524d03518 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 06:25:23 -0500 Subject: [PATCH 08/14] Add default initializers in ImportIsotope --- src/celeritas/io/ImportElement.hh | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/celeritas/io/ImportElement.hh b/src/celeritas/io/ImportElement.hh index 08204c2d80..845c77e828 100644 --- a/src/celeritas/io/ImportElement.hh +++ b/src/celeritas/io/ImportElement.hh @@ -22,12 +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 (BE) [MeV] - double proton_loss_energy; //!< BE(A, Z) - BE(A-1, Z-1) [MeV] - double neutron_loss_energy; //!< BE(A, Z) - BE(A-1, Z) [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] }; //---------------------------------------------------------------------------// From 3bfad58b09f17b4fd15a5490ea8d4f000885359a Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 06:29:19 -0500 Subject: [PATCH 09/14] Use CELER_JSON_LOAD_OPTION from JasonUtils.json --- src/celeritas/neutron/model/CascadeOptionsIO.json.cc | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/celeritas/neutron/model/CascadeOptionsIO.json.cc b/src/celeritas/neutron/model/CascadeOptionsIO.json.cc index 5b2f1b4220..bc9f3995fa 100644 --- a/src/celeritas/neutron/model/CascadeOptionsIO.json.cc +++ b/src/celeritas/neutron/model/CascadeOptionsIO.json.cc @@ -22,12 +22,7 @@ namespace celeritas */ void from_json(nlohmann::json const& j, CascadeOptions& opts) { -#define FDO_INPUT(NAME) \ - do \ - { \ - if (j.contains(#NAME)) \ - j.at(#NAME).get_to(opts.NAME); \ - } while (0) +#define FDO_INPUT(NAME) CELER_JSON_LOAD_OPTION(j, opts, NAME) FDO_INPUT(use_precompound); FDO_INPUT(use_abla); From 78cd8988ee6a1ec34085f596d7747778ebc7107d Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 06:31:05 -0500 Subject: [PATCH 10/14] Rename zone_builder to build_nuclear_zones --- src/celeritas/neutron/model/NeutronInelasticModel.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/celeritas/neutron/model/NeutronInelasticModel.cc b/src/celeritas/neutron/model/NeutronInelasticModel.cc index 1ba144a954..54b9cde458 100644 --- a/src/celeritas/neutron/model/NeutronInelasticModel.cc +++ b/src/celeritas/neutron/model/NeutronInelasticModel.cc @@ -85,12 +85,12 @@ NeutronInelasticModel::NeutronInelasticModel(ActionId id, CELER_ASSERT(data.xs_params.size() == data.nucleon_xs.size()); // Build (A, Z)-dependent nuclear zone data - detail::NuclearZoneBuilder zone_builder( + detail::NuclearZoneBuilder build_nuclear_zones( options, data.scalars, &data.nuclear_zones); for (auto iso_id : range(IsotopeId{materials.num_isotopes()})) { - zone_builder(materials.get(iso_id)); + build_nuclear_zones(materials.get(iso_id)); } CELER_ASSERT(data.nuclear_zones.zones.size() == materials.num_isotopes()); From 7a200a03df3c77c962aade3558d0573d24e86733 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 07:14:15 -0500 Subject: [PATCH 11/14] Reflect review comments and suggestions by sethrj --- .../model/detail/NuclearZoneBuilder.hh | 90 +++++++++---------- 1 file changed, 44 insertions(+), 46 deletions(-) diff --git a/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh index 804e87cc43..2814934e89 100644 --- a/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh +++ b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh @@ -8,7 +8,9 @@ #pragma once #include +#include +#include "corecel/math/SoftEqual.hh" #include "celeritas/Types.hh" #include "celeritas/mat/IsotopeView.hh" #include "celeritas/neutron/data/NeutronInelasticData.hh" @@ -38,9 +40,9 @@ class NuclearZoneBuilder public: // Construct with cascade options and shared data - explicit inline NuclearZoneBuilder(CascadeOptions const& options, - NeutronInelasticScalars const& scalars, - Data* 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); @@ -138,7 +140,15 @@ void NuclearZoneBuilder::operator()(IsotopeView const& target) */ size_type NuclearZoneBuilder::num_nuclear_zones(AtomicMassNumber a) { - return (a.get() < 5) ? 1 : (a.get() < 100) ? 3 : 6; + if (a < AtomicMassNumber{5}) + { + return 1; + } + if (a < AtomicMassNumber{100}) + { + return 3; + } + return 6; } //---------------------------------------------------------------------------// @@ -161,53 +171,43 @@ void NuclearZoneBuilder::calc_zone_component(IsotopeView const& target, // Temporary data for the zone-by-zone density function size_type num_of_zones = components.size(); - std::vector ur(num_of_zones + 1); std::vector integral(num_of_zones); - real_type total_integral{0}; // Fill the nuclear radius by each zone auto amass = a.get(); - ur[0] = (amass < 12) ? 0 : -skin_ratio; + 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; - total_integral = integral[0]; } - else if (amass < 100) + else if (amass < 12) { - constexpr Real3 alpha = {0.7, 0.3, 0.01}; - 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)) { - // 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}); - for (auto i : range(num_of_zones)) - { - real_type y = std::sqrt(-std::log(alpha[i])); - components[i].radius = gauss_radius * y; - ur[i + 1] = y; - integral[i] - = this->integrate_gaussian(ur[i], ur[i + 1], gauss_radius); - total_integral += integral[i]; - } + components[i].radius = gauss_radius * y[i]; + integral[i] = this->integrate_gaussian(ymin, y[i], gauss_radius); + ymin = y[i]; } - else + } + 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)) { - // Intermediate nuclei have a three-zone Woods-Saxon potential - 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; - ur[i + 1] = y; - integral[i] = this->integrate_woods_saxon( - ur[i], ur[i + 1], nuclear_radius); - total_integral += integral[i]; - } + 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 @@ -219,14 +219,12 @@ void NuclearZoneBuilder::calc_zone_component(IsotopeView const& target, { real_type y = std::log((1 + skin_decay) / alpha[i] - 1); components[i].radius = nuclear_radius + skin_depth_ * y; - ur[i + 1] = y; - integral[i] = this->integrate_woods_saxon( - ur[i], ur[i + 1], nuclear_radius); - total_integral += integral[i]; + integral[i] = this->integrate_woods_saxon(ymin, y, nuclear_radius); + ymin = y; } } - - ur.clear(); + 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)) @@ -263,8 +261,6 @@ void NuclearZoneBuilder::calc_zone_component(IsotopeView const& target, + dm[ptype]; } } - - integral.clear(); } //---------------------------------------------------------------------------// @@ -351,7 +347,7 @@ real_type NuclearZoneBuilder::integrate_gaussian(real_type rmin, = this->half() * delta_r * (rmin_sq * std::exp(-rmin_sq) + rmax_sq * std::exp(-rmax_sq)); - real_type result = integrate_potential(false, rmin, delta_r, integral); + real_type result = integrate_potential(rmin, delta_r, integral); return ipow<3>(gauss_radius) * result; } @@ -374,6 +370,8 @@ real_type NuclearZoneBuilder::integrate_potential(real_type rmin, bool succeeded = false; int remaining_trials = max_trials; + SoftEqual const soft_eq{epsilon}; + do { delta_r *= this->half(); @@ -390,7 +388,7 @@ real_type NuclearZoneBuilder::integrate_potential(real_type rmin, result = this->half() * integral + fi * delta_r; - if (std::fabs(1 - integral / result) < epsilon) + if (soft_eq(1, integral / result)) { succeeded = true; } From 68253c41dfc8f87f4f9ec086863f40942e4b3f38 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 07:17:45 -0500 Subject: [PATCH 12/14] Add a test case for the Gaussian ptential with isotope A in [5,11] --- .../neutron/NeutronInelastic.test.cc | 53 +++++++++++++++++++ test/celeritas/neutron/NeutronTestBase.cc | 16 +++++- 2 files changed, 68 insertions(+), 1 deletion(-) diff --git a/test/celeritas/neutron/NeutronInelastic.test.cc b/test/celeritas/neutron/NeutronInelastic.test.cc index 0b10dbce7a..12b1dba8eb 100644 --- a/test/celeritas/neutron/NeutronInelastic.test.cc +++ b/test/celeritas/neutron/NeutronInelastic.test.cc @@ -321,6 +321,59 @@ TEST_F(NeutronInelasticTest, model_data) 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.num_zones, 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); } //---------------------------------------------------------------------------// diff --git a/test/celeritas/neutron/NeutronTestBase.cc b/test/celeritas/neutron/NeutronTestBase.cc index 289a1b5271..7816942c7f 100644 --- a/test/celeritas/neutron/NeutronTestBase.cc +++ b/test/celeritas/neutron/NeutronTestBase.cc @@ -95,7 +95,21 @@ NeutronTestBase::NeutronTestBase() units::MevEnergy{8.004}, units::MevEnergy{7.368}, units::MevMass{193687.6762}, - "208Pb"}}; + "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}, From 77fd2d7646506c0886e8c31041a05253aba79c00 Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 12:23:11 -0500 Subject: [PATCH 13/14] Remove num_zones from NuclearZones and use zones.size() instead --- src/celeritas/neutron/data/NeutronInelasticData.hh | 6 +----- test/celeritas/neutron/NeutronInelastic.test.cc | 8 ++++---- 2 files changed, 5 insertions(+), 9 deletions(-) diff --git a/src/celeritas/neutron/data/NeutronInelasticData.hh b/src/celeritas/neutron/data/NeutronInelasticData.hh index 6be72169bd..0010dee34e 100644 --- a/src/celeritas/neutron/data/NeutronInelasticData.hh +++ b/src/celeritas/neutron/data/NeutronInelasticData.hh @@ -84,14 +84,10 @@ struct ZoneComponent */ struct NuclearZones { - size_type num_zones{}; ItemRange zones; //! Whether all data are assigned and valid - explicit CELER_FUNCTION operator bool() const - { - return num_zones > 0 && !zones.empty(); - } + explicit CELER_FUNCTION operator bool() const { return !zones.empty(); } }; //---------------------------------------------------------------------------// diff --git a/test/celeritas/neutron/NeutronInelastic.test.cc b/test/celeritas/neutron/NeutronInelastic.test.cc index 12b1dba8eb..ad6a2fe785 100644 --- a/test/celeritas/neutron/NeutronInelastic.test.cc +++ b/test/celeritas/neutron/NeutronInelastic.test.cc @@ -190,7 +190,7 @@ TEST_F(NeutronInelasticTest, model_data) // Check the size of the number of nuclear zones NuclearZones he3_nuclear_zones = shared.nuclear_zones.zones[iso_id]; - EXPECT_EQ(he3_nuclear_zones.num_zones, 1); + EXPECT_EQ(he3_nuclear_zones.zones.size(), 1); // Check zone data for (auto sid : he3_nuclear_zones.zones) @@ -213,7 +213,7 @@ TEST_F(NeutronInelasticTest, model_data) // Check the size of the number of nuclear zones NuclearZones cu63_nuclear_zones = shared.nuclear_zones.zones[iso_cu63]; - EXPECT_EQ(cu63_nuclear_zones.num_zones, 3); + EXPECT_EQ(cu63_nuclear_zones.zones.size(), 3); // Check zone data std::vector radii; @@ -266,7 +266,7 @@ TEST_F(NeutronInelasticTest, model_data) // Check the size of the number of nuclear zones NuclearZones pb208_nuclear_zones = shared.nuclear_zones.zones[iso_pb208]; - EXPECT_EQ(pb208_nuclear_zones.num_zones, 6); + EXPECT_EQ(pb208_nuclear_zones.zones.size(), 6); // Check zone data radii.clear(); @@ -327,7 +327,7 @@ TEST_F(NeutronInelasticTest, model_data) // Check the size of the number of nuclear zones NuclearZones b11_nuclear_zones = shared.nuclear_zones.zones[iso_b11]; - EXPECT_EQ(b11_nuclear_zones.num_zones, 3); + EXPECT_EQ(b11_nuclear_zones.zones.size(), 3); // Clear zone data radii.clear(); From 37dc2c7dd6025a768629ad05c5d20634a6beaded Mon Sep 17 00:00:00 2001 From: Soon Yung Jun Date: Thu, 13 Jun 2024 12:51:52 -0500 Subject: [PATCH 14/14] Reflect more review comments by amandalund and sethrj --- .../neutron/model/detail/NuclearZoneBuilder.hh | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh index 2814934e89..96b835982d 100644 --- a/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh +++ b/src/celeritas/neutron/model/detail/NuclearZoneBuilder.hh @@ -101,7 +101,6 @@ class NuclearZoneBuilder /*! * Construct with cascade options and data. */ -CELER_FUNCTION NuclearZoneBuilder::NuclearZoneBuilder(CascadeOptions const& options, NeutronInelasticScalars const& scalars, Data* data) @@ -129,7 +128,6 @@ void NuclearZoneBuilder::operator()(IsotopeView const& target) this->calc_zone_component(target, comp); NuclearZones nucl_zones; - nucl_zones.num_zones = num_zones; nucl_zones.zones = components_.insert_back(comp.begin(), comp.end()); zones_.push_back(nucl_zones); } @@ -157,7 +155,6 @@ size_type NuclearZoneBuilder::num_nuclear_zones(AtomicMassNumber a) * 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). */ -CELER_FUNCTION void NuclearZoneBuilder::calc_zone_component(IsotopeView const& target, ComponentVec& components) { @@ -214,7 +211,6 @@ void NuclearZoneBuilder::calc_zone_component(IsotopeView const& target, { // 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); @@ -238,16 +234,17 @@ void NuclearZoneBuilder::calc_zone_component(IsotopeView const& target, } // Fill the nuclear zone density, fermi momentum, and potential - int const dim = ZoneComponent::NucleonArray::size(); int num_protons = target.atomic_number().get(); - int num_nucleons[dim] = {num_protons, a.get() - num_protons}; - real_type mass[dim] = {proton_mass_.value(), neutron_mass_.value()}; - real_type dm[dim] = {target.proton_loss_energy().value(), - target.neutron_loss_energy().value()}; + 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(dim)) + for (auto ptype : range(ZoneComponent::NucleonArray::size())) { components[i].density[ptype] = static_cast(num_nucleons[ptype]) * integral[i]