From f474fd8f13cfde960d127dcc1c8e914bbf09a181 Mon Sep 17 00:00:00 2001 From: Luke Shingles Date: Thu, 12 Sep 2024 10:08:00 +0100 Subject: [PATCH] Add kappa_gamma estimators --- gammapkt.cc | 17 +++++++++++++++++ globals.h | 5 +++++ sn3d.cc | 20 ++++++++++++++++++++ 3 files changed, 42 insertions(+) diff --git a/gammapkt.cc b/gammapkt.cc index 720815db4..11159ebab 100644 --- a/gammapkt.cc +++ b/gammapkt.cc @@ -647,6 +647,16 @@ constexpr auto meanf_sigma(const double x) -> double { return tot; } +// get the gamma-ray opacity (with the expected energy loss per interaction factor included) +auto get_kappa(const Packet &pkt) -> double { + const double xx = H * pkt.nu_cmf / ME / CLIGHT / CLIGHT; + const int mgi = grid::get_cell_modelgridindex(pkt.where); + const double chi = ((meanf_sigma(xx) * grid::get_nnetot(mgi)) + get_chi_photo_electric_rf(pkt) + + (sigma_pair_prod_rf(pkt) * (1. - (2.46636e+20 / pkt.nu_cmf)))); + const double kappa = 1 / grid::get_rho(mgi) * chi; + return kappa; +} + // Subroutine to record the heating rate in a cell due to gamma rays. // By heating rate I mean, for now, really the rate at which the code is making // k-packets in that cell which will then convert into r-packets. This is (going @@ -677,6 +687,9 @@ void update_gamma_dep(const Packet &pkt, const double dist, const int mgi, const assert_testmodeonly(heating_cont >= 0.); assert_testmodeonly(std::isfinite(heating_cont)); atomicadd(globals::dep_estimator_gamma[nonemptymgi], heating_cont); + + const double kappa_E_cont = heating_cont * get_kappa(pkt); + atomicadd(globals::estimator_gamma_kappa_integrated, kappa_E_cont); } // handle physical pair production event @@ -1079,6 +1092,9 @@ __host__ __device__ void pellet_gamma_decay(Packet &pkt) { // printout("initialise pol state of packet %g, %g, %g, %g, // %g\n",pkt.stokes_qu[0],pkt.stokes_qu[1],pkt.pol_dir[0],pkt.pol_dir[1],pkt.pol_dir[2]); // printout("pkt direction %g, %g, %g\n",pkt.dir[0],pkt.dir[1],pkt.dir[2]); + + atomicadd(globals::estimator_gamma_kappa_decayspec, get_kappa(pkt) * pkt.e_cmf); + atomicadd(globals::estimator_gamma_nu_cmf_decayspec, pkt.nu_cmf * pkt.e_cmf); } __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) { @@ -1096,6 +1112,7 @@ __host__ __device__ void do_gamma(Packet &pkt, const int nts, const double t2) { if (pkt.type != TYPE_GAMMA && pkt.type != TYPE_ESCAPE) { atomicadd(globals::timesteps[nts].gamma_dep_discrete, pkt.e_cmf); + atomicadd(globals::estimator_gamma_kappa_absorbedspec, get_kappa(pkt) * pkt.e_cmf); if constexpr (GAMMA_THERMALISATION_SCHEME != ThermalisationScheme::DETAILED) { // no transport, so the path-based gamma deposition estimator won't get updated unless we do it here diff --git a/globals.h b/globals.h index d46106e78..2b1a2eb7c 100644 --- a/globals.h +++ b/globals.h @@ -212,6 +212,11 @@ inline std::vector dep_estimator_positron; inline std::vector dep_estimator_electron; inline std::vector dep_estimator_alpha; +inline double estimator_gamma_nu_cmf_decayspec; +inline double estimator_gamma_kappa_integrated; +inline double estimator_gamma_kappa_decayspec; +inline double estimator_gamma_kappa_absorbedspec; + inline int bfestimcount{0}; // for USE_LUT_PHOTOION = true diff --git a/sn3d.cc b/sn3d.cc index 568fe2733..c936a9776 100644 --- a/sn3d.cc +++ b/sn3d.cc @@ -617,6 +617,11 @@ void zero_estimators() { std::ranges::fill(globals::dep_estimator_electron, 0.); std::ranges::fill(globals::dep_estimator_alpha, 0.); + globals::estimator_gamma_nu_cmf_decayspec = 0.; + globals::estimator_gamma_kappa_integrated = 0.; + globals::estimator_gamma_kappa_decayspec = 0.; + globals::estimator_gamma_kappa_absorbedspec = 0.; + if constexpr (USE_LUT_PHOTOION) { if (globals::nbfcontinua_ground > 0) { std::ranges::fill(globals::gammaestimator, 0.); @@ -643,6 +648,21 @@ void normalise_deposition_estimators(int nts) { globals::timesteps[nts].electron_dep = 0.; globals::timesteps[nts].alpha_dep = 0.; + double gamma_dep_erg = 0.; + for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { + gamma_dep_erg += (globals::dep_estimator_gamma[nonemptymgi] / nprocs); + } + globals::estimator_gamma_kappa_integrated /= gamma_dep_erg; + globals::estimator_gamma_kappa_decayspec /= globals::timesteps[nts].gamma_emission; + globals::estimator_gamma_kappa_absorbedspec /= globals::timesteps[nts].gamma_dep_discrete; + globals::estimator_gamma_nu_cmf_decayspec /= globals::timesteps[nts].gamma_emission; + const double mean_en_mev = H * globals::estimator_gamma_nu_cmf_decayspec / MEV; + + printout( + "timestep %d kappa_eff [cm2/g]: path-integrated %g decay_spec %g deposit_spec %g decay gamma mean E[MeV] %g\n", + nts, globals::estimator_gamma_kappa_integrated, globals::estimator_gamma_kappa_decayspec, + globals::estimator_gamma_kappa_absorbedspec, mean_en_mev); + for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) { const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);