From 01780d3ae202e7d64644155f48905fc5da43a707 Mon Sep 17 00:00:00 2001 From: Gerrit Leck Date: Tue, 5 Nov 2024 16:18:06 +0100 Subject: [PATCH 1/2] Get constant neutrino loss fraction --- artisoptions_kilonova_lte.h | 4 ++++ decay.cc | 14 ++++++++++---- globals.h | 2 ++ packet.cc | 2 ++ 4 files changed, 18 insertions(+), 4 deletions(-) diff --git a/artisoptions_kilonova_lte.h b/artisoptions_kilonova_lte.h index 3f35ddf43..ce75d2c0a 100644 --- a/artisoptions_kilonova_lte.h +++ b/artisoptions_kilonova_lte.h @@ -151,5 +151,9 @@ constexpr auto PARTICLE_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; constexpr auto GAMMA_THERMALISATION_SCHEME = ThermalisationScheme::DETAILED; +constexpr float BETAMINUS_ENERGY_GAMMA_SPLITUP = 0.45; + +constexpr float BETAMINUS_ENERGY_ELECTRON_SPLITUP = 0.2; + // NOLINTEND(modernize*,misc-unused-parameters) #endif // ARTISOPTIONS_H diff --git a/decay.cc b/decay.cc index 6326a684e..e0238fec4 100644 --- a/decay.cc +++ b/decay.cc @@ -232,8 +232,10 @@ void printout_nuclidemeanlife(const int z, const int a) { // contributed energy release per decay [erg] for decaytype (e.g. decaytypes::DECAYTYPE_BETAPLUS) (excludes neutrinos!) [[nodiscard]] auto nucdecayenergy(const int nucindex, const int decaytype) -> double { - const double endecay = nuclides[nucindex].endecay_gamma + nucdecayenergyparticle(nucindex, decaytype); - + double endecay = nuclides[nucindex].endecay_gamma + nucdecayenergyparticle(nucindex, decaytype); + if(globals::packet_setup_phase && decaytype == DECAYTYPE_BETAMINUS) { + endecay = (BETAMINUS_ENERGY_GAMMA_SPLITUP + BETAMINUS_ENERGY_ELECTRON_SPLITUP) * nuclides[nucindex].endecay_q[decaytype]; + } return endecay; } @@ -1408,8 +1410,12 @@ void setup_radioactive_pellet(const double e0, const int mgi, Packet &pkt) { pkt.pellet_nucindex = nucindex; pkt.pellet_decaytype = decaytype; - const auto engamma = nucdecayenergygamma(nucindex); - const auto enparticle = nucdecayenergyparticle(nucindex, decaytype); + auto engamma = nucdecayenergygamma(nucindex); + auto enparticle = nucdecayenergyparticle(nucindex, decaytype); + if(pkt.pellet_decaytype == DECAYTYPE_BETAMINUS) { + engamma = BETAMINUS_ENERGY_GAMMA_SPLITUP * nuclides[nucindex].endecay_q[DECAYTYPE_BETAMINUS]; + enparticle = BETAMINUS_ENERGY_ELECTRON_SPLITUP * nuclides[nucindex].endecay_q[DECAYTYPE_BETAMINUS]; + } pkt.originated_from_particlenotgamma = (rng_uniform() >= engamma / (engamma + enparticle)); pkt.nu_cmf = enparticle / H; // will be overwritten for gamma rays, but affects the thermalisation of particles diff --git a/globals.h b/globals.h index ee8a013e5..75a9c76f0 100644 --- a/globals.h +++ b/globals.h @@ -318,6 +318,8 @@ inline int num_grey_timesteps; inline int n_titer; inline bool lte_iteration; +inline bool packet_setup_phase{false}; + inline std::deque mutex_cellcachemacroatom; inline void setup_mpi_vars() { diff --git a/packet.cc b/packet.cc index 812da7e0b..8f49e7653 100644 --- a/packet.cc +++ b/packet.cc @@ -100,7 +100,9 @@ void packet_init(Packet *pkt) const double e0_tinf = etot_tinf / globals::npkts; printout("packet e0 (t_0 to t_inf) %g erg\n", e0_tinf); + globals::packet_setup_phase = true; decay::setup_decaypath_energy_per_mass(); + globals::packet_setup_phase = false; // Need to get a normalisation factor auto en_cumulative = std::vector(grid::ngrid); From e68491594aef5a434005f93579e70b0310690736 Mon Sep 17 00:00:00 2001 From: Gerrit Leck Date: Thu, 7 Nov 2024 16:37:23 +0100 Subject: [PATCH 2/2] update --- artisoptions_kilonova_lte.h | 2 +- decay.cc | 20 ++++++++++---------- gammapkt.cc | 5 ++++- globals.h | 2 -- packet.cc | 2 -- 5 files changed, 15 insertions(+), 16 deletions(-) diff --git a/artisoptions_kilonova_lte.h b/artisoptions_kilonova_lte.h index ce75d2c0a..9de58357f 100644 --- a/artisoptions_kilonova_lte.h +++ b/artisoptions_kilonova_lte.h @@ -6,7 +6,7 @@ #include "constants.h" -constexpr int MPKTS = 15000; +constexpr int MPKTS = 100000; constexpr auto GRID_TYPE = GridType::CARTESIAN3D; constexpr int CUBOID_NCOORDGRID_X = 50; diff --git a/decay.cc b/decay.cc index e0238fec4..8761dd0de 100644 --- a/decay.cc +++ b/decay.cc @@ -232,10 +232,8 @@ void printout_nuclidemeanlife(const int z, const int a) { // contributed energy release per decay [erg] for decaytype (e.g. decaytypes::DECAYTYPE_BETAPLUS) (excludes neutrinos!) [[nodiscard]] auto nucdecayenergy(const int nucindex, const int decaytype) -> double { - double endecay = nuclides[nucindex].endecay_gamma + nucdecayenergyparticle(nucindex, decaytype); - if(globals::packet_setup_phase && decaytype == DECAYTYPE_BETAMINUS) { - endecay = (BETAMINUS_ENERGY_GAMMA_SPLITUP + BETAMINUS_ENERGY_ELECTRON_SPLITUP) * nuclides[nucindex].endecay_q[decaytype]; - } + const double endecay = nuclides[nucindex].endecay_gamma + nucdecayenergyparticle(nucindex, decaytype); + return endecay; } @@ -1001,6 +999,12 @@ void init_nuclides(const std::vector &custom_zlist, const std::vector // Read in data for gamma ray lines and make a list of them in energy order. gammapkt::init_gamma_data(); + // manipulate betaminus decay splitup ratios here + for (int i = 0; i < static_cast(nuclides.size()); i++) { + nuclides[i].endecay_gamma = BETAMINUS_ENERGY_GAMMA_SPLITUP * nuclides[i].endecay_q[DECAYTYPE_BETAMINUS]; + nuclides[i].endecay_electron = BETAMINUS_ENERGY_ELECTRON_SPLITUP * nuclides[i].endecay_q[DECAYTYPE_BETAMINUS]; + } + // TODO: generalise this to all included nuclides printout("decayenergy(NI56), decayenergy(CO56), decayenergy_gamma(CO56): %g, %g, %g\n", nucdecayenergytotal(28, 56) / MEV, nucdecayenergytotal(27, 56) / MEV, nucdecayenergygamma(27, 56) / MEV); @@ -1410,12 +1414,8 @@ void setup_radioactive_pellet(const double e0, const int mgi, Packet &pkt) { pkt.pellet_nucindex = nucindex; pkt.pellet_decaytype = decaytype; - auto engamma = nucdecayenergygamma(nucindex); - auto enparticle = nucdecayenergyparticle(nucindex, decaytype); - if(pkt.pellet_decaytype == DECAYTYPE_BETAMINUS) { - engamma = BETAMINUS_ENERGY_GAMMA_SPLITUP * nuclides[nucindex].endecay_q[DECAYTYPE_BETAMINUS]; - enparticle = BETAMINUS_ENERGY_ELECTRON_SPLITUP * nuclides[nucindex].endecay_q[DECAYTYPE_BETAMINUS]; - } + const auto engamma = nucdecayenergygamma(nucindex); + const auto enparticle = nucdecayenergyparticle(nucindex, decaytype); pkt.originated_from_particlenotgamma = (rng_uniform() >= engamma / (engamma + enparticle)); pkt.nu_cmf = enparticle / H; // will be overwritten for gamma rays, but affects the thermalisation of particles diff --git a/gammapkt.cc b/gammapkt.cc index 56893e361..65b7a3e8a 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -236,7 +236,10 @@ void init_xcom_photoion_data() { __host__ __device__ auto choose_gamma_ray(const int nucindex) -> double { // Routine to choose which gamma ray line it'll be. - const double E_gamma = decay::nucdecayenergygamma(nucindex); // Average energy per gamma line of a decay + double E_gamma = 0.; // Average energy per gamma line of a decay + for (ptrdiff_t n = 0; n < std::ssize(gamma_spectra[nucindex]); n++) { + E_gamma += gamma_spectra[nucindex][n].probability * gamma_spectra[nucindex][n].energy; + } const double zrand = rng_uniform(); double runtot = 0.; diff --git a/globals.h b/globals.h index 75a9c76f0..ee8a013e5 100644 --- a/globals.h +++ b/globals.h @@ -318,8 +318,6 @@ inline int num_grey_timesteps; inline int n_titer; inline bool lte_iteration; -inline bool packet_setup_phase{false}; - inline std::deque mutex_cellcachemacroatom; inline void setup_mpi_vars() { diff --git a/packet.cc b/packet.cc index 8f49e7653..812da7e0b 100644 --- a/packet.cc +++ b/packet.cc @@ -100,9 +100,7 @@ void packet_init(Packet *pkt) const double e0_tinf = etot_tinf / globals::npkts; printout("packet e0 (t_0 to t_inf) %g erg\n", e0_tinf); - globals::packet_setup_phase = true; decay::setup_decaypath_energy_per_mass(); - globals::packet_setup_phase = false; // Need to get a normalisation factor auto en_cumulative = std::vector(grid::ngrid);