Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Calculate gamma kappa_eff and <E_gamma> #125

Draft
wants to merge 1 commit into
base: develop
Choose a base branch
from
Draft

Calculate gamma kappa_eff and <E_gamma> #125

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 15 additions & 0 deletions gammapkt.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -1079,6 +1089,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) {
Expand All @@ -1096,6 +1109,8 @@ __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);
atomicadd(globals::estimator_gamma_nu_cmf_absorbedspec, pkt.nu_cmf * 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
Expand Down
5 changes: 5 additions & 0 deletions globals.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,6 +212,11 @@ inline std::vector<double> dep_estimator_positron;
inline std::vector<double> dep_estimator_electron;
inline std::vector<double> dep_estimator_alpha;

inline double estimator_gamma_kappa_decayspec = 0.;
inline double estimator_gamma_nu_cmf_decayspec = 0.;
inline double estimator_gamma_kappa_absorbedspec = 0.;
inline double estimator_gamma_nu_cmf_absorbedspec = 0.;

inline int bfestimcount{0};

// for USE_LUT_PHOTOION = true
Expand Down
19 changes: 19 additions & 0 deletions sn3d.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_nu_cmf_absorbedspec = 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.);
Expand All @@ -643,6 +648,20 @@ void normalise_deposition_estimators(int nts) {
globals::timesteps[nts].electron_dep = 0.;
globals::timesteps[nts].alpha_dep = 0.;

globals::estimator_gamma_kappa_decayspec /= globals::timesteps[nts].gamma_emission;
globals::estimator_gamma_nu_cmf_decayspec /= globals::timesteps[nts].gamma_emission;
const double mean_en_decay_mev = H * globals::estimator_gamma_nu_cmf_decayspec / MEV;

globals::estimator_gamma_kappa_absorbedspec /= globals::timesteps[nts].gamma_dep_discrete;
globals::estimator_gamma_nu_cmf_absorbedspec /= globals::timesteps[nts].gamma_dep_discrete;
const double mean_en_absorb_mev = H * globals::estimator_gamma_nu_cmf_absorbedspec / MEV;

printout(
"timestep %d %.2f days kappa_eff [cm2/g]: decay_spec %.3f deposit_spec %.3f. mean gamma "
"E[MeV]: decay %.2f deposit %.2f\n",
nts, globals::timesteps[nts].mid / DAY, globals::estimator_gamma_kappa_decayspec,
globals::estimator_gamma_kappa_absorbedspec, mean_en_decay_mev, mean_en_absorb_mev);

for (int nonemptymgi = 0; nonemptymgi < grid::get_nonempty_npts_model(); nonemptymgi++) {
const int mgi = grid::get_mgi_of_nonemptymgi(nonemptymgi);

Expand Down