Skip to content

Commit

Permalink
Pre-calculate macroatom stats for n lowest levels
Browse files Browse the repository at this point in the history
  • Loading branch information
@lukeshingles
lukeshingles committed Sep 11, 2024
1 parent 94b09bf commit 4f17aa1
Show file tree
Hide file tree
Showing 3 changed files with 139 additions and 118 deletions.
237 changes: 119 additions & 118 deletions macroatom.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,123 +34,6 @@ constexpr bool LOG_MACROATOM = false;

FILE *macroatom_file{};

auto calculate_macroatom_transitionrates(const int modelgridindex, const int element, const int ion, const int level,
const double t_mid, CellCacheLevels &chlevel) {
// printout("Calculating transition rates for element %d ion %d level %d\n", element, ion, level);
auto processrates = std::array<double, MA_ACTION_COUNT>{};

const auto T_e = grid::get_Te(modelgridindex);
const auto nne = grid::get_nne(modelgridindex);
const double epsilon_current = epsilon(element, ion, level);
const double statweight = stat_weight(element, ion, level);
const auto nnlevel = get_levelpop(modelgridindex, element, ion, level);

// Downward transitions within the current ionisation stage:
// radiative/collisional deexcitation and internal downward jumps
double sum_internal_down_same = 0.;
double sum_raddeexc = 0.;
double sum_coldeexc = 0.;
const int ndowntrans = get_ndowntrans(element, ion, level);
const auto *const leveldowntranslist = get_downtranslist(element, ion, level);
auto *const arr_sum_epstrans_rad_deexc = chlevel.sum_epstrans_rad_deexc;
auto *const arr_sum_internal_down_same = chlevel.sum_internal_down_same;
for (int i = 0; i < ndowntrans; i++) {
const auto &downtrans = leveldowntranslist[i];
const int lower = downtrans.targetlevelindex;
const auto A_ul = downtrans.einstein_A;
const double epsilon_target = epsilon(element, ion, lower);
const double epsilon_trans = epsilon_current - epsilon_target;

const double R = rad_deexcitation_ratecoeff(modelgridindex, element, ion, lower, epsilon_trans, A_ul, statweight,
nnlevel, t_mid);
const double C = col_deexcitation_ratecoeff(T_e, nne, epsilon_trans, element, ion, level, downtrans);

sum_raddeexc += R * epsilon_trans;
sum_coldeexc += C * epsilon_trans;
sum_internal_down_same += (R + C) * epsilon_target;

arr_sum_epstrans_rad_deexc[i] = sum_raddeexc;
arr_sum_internal_down_same[i] = sum_internal_down_same;

// printout("checking downtrans %d to level %d: R %g, C %g, epsilon_trans %g\n",i,lower,R,C,epsilon_trans);
}
processrates[MA_ACTION_RADDEEXC] = sum_raddeexc;
processrates[MA_ACTION_COLDEEXC] = sum_coldeexc;
processrates[MA_ACTION_INTERNALDOWNSAME] = sum_internal_down_same;

// Calculate sum for upward internal transitions
// transitions within the current ionisation stage
double sum_internal_up_same = 0.;
const int nuptrans = get_nuptrans(element, ion, level);
const auto *const uptranslist = get_uptranslist(element, ion, level);
for (int i = 0; i < nuptrans; i++) {
const auto &uptrans = uptranslist[i];
const double epsilon_trans = epsilon(element, ion, uptrans.targetlevelindex) - epsilon_current;

const double R = rad_excitation_ratecoeff(modelgridindex, element, ion, level, i, epsilon_trans, nnlevel,
uptrans.lineindex, t_mid);
const double C = col_excitation_ratecoeff(T_e, nne, element, ion, level, i, epsilon_trans, statweight);
const double NT = nonthermal::nt_excitation_ratecoeff(modelgridindex, element, ion, level, i, uptrans.lineindex);

sum_internal_up_same += (R + C + NT) * epsilon_current;
chlevel.sum_internal_up_same[i] = sum_internal_up_same;
}
processrates[MA_ACTION_INTERNALUPSAME] = sum_internal_up_same;

assert_always(std::isfinite(processrates[MA_ACTION_INTERNALUPSAME]));

// Downward transitions to lower ionisation stages:
// radiative/collisional recombination and internal downward jumps
// checks only if there is a lower ion, doesn't make sure that Z(ion)=Z(ion-1)+1
double sum_internal_down_lower = 0.;
double sum_radrecomb = 0.;
double sum_colrecomb = 0.;
if (ion > 0 && level <= get_maxrecombininglevel(element, ion)) {
// nlevels = get_nlevels(element,ion-1);
const int nlevels = get_ionisinglevels(element, ion - 1);
// nlevels = get_ionisinglevels(element,ion-1);
for (int lower = 0; lower < nlevels; lower++) {
const double epsilon_target = epsilon(element, ion - 1, lower);
const double epsilon_trans = epsilon_current - epsilon_target;

const double R = rad_recombination_ratecoeff(T_e, nne, element, ion, level, lower, modelgridindex);
const double C = col_recombination_ratecoeff(modelgridindex, element, ion, level, lower, epsilon_trans);

sum_internal_down_lower += (R + C) * epsilon_target;

sum_radrecomb += R * epsilon_trans;
sum_colrecomb += C * epsilon_trans;
}
}
processrates[MA_ACTION_INTERNALDOWNLOWER] = sum_internal_down_lower;
processrates[MA_ACTION_RADRECOMB] = sum_radrecomb;
processrates[MA_ACTION_COLRECOMB] = sum_colrecomb;

// Transitions to higher ionisation stages
double sum_up_highernt = 0.;
double sum_up_higher = 0.;
const int ionisinglevels = get_ionisinglevels(element, ion);
if (ion < get_nions(element) - 1 && level < ionisinglevels) {
if (NT_ON) {
sum_up_highernt = nonthermal::nt_ionization_ratecoeff(modelgridindex, element, ion) * epsilon_current;
}

const auto nphixstargets = get_nphixstargets(element, ion, level);
for (int phixstargetindex = 0; phixstargetindex < nphixstargets; phixstargetindex++) {
const double epsilon_trans = get_phixs_threshold(element, ion, level, phixstargetindex);

const double R = get_corrphotoioncoeff(element, ion, level, phixstargetindex, modelgridindex);
const double C = col_ionization_ratecoeff(T_e, nne, element, ion, level, phixstargetindex, epsilon_trans);

sum_up_higher += (R + C) * epsilon_current;
}
}
processrates[MA_ACTION_INTERNALUPHIGHERNT] = sum_up_highernt;
processrates[MA_ACTION_INTERNALUPHIGHER] = sum_up_higher;

return processrates;
}

auto do_macroatom_internal_down_same(const int element, const int ion, const int level,
const CellCacheLevels &chlevel) -> int {
const int ndowntrans = get_ndowntrans(element, ion, level);
Expand Down Expand Up @@ -306,6 +189,124 @@ void do_macroatom_raddeexcitation(Packet &pkt, const int element, const int ion,

} // anonymous namespace

auto calculate_macroatom_transitionrates(const int modelgridindex, const int element, const int ion, const int level,
const double t_mid,
CellCacheLevels &chlevel) -> std::array<double, MA_ACTION_COUNT> {
// printout("Calculating transition rates for element %d ion %d level %d\n", element, ion, level);
auto processrates = std::array<double, MA_ACTION_COUNT>{};

const auto T_e = grid::get_Te(modelgridindex);
const auto nne = grid::get_nne(modelgridindex);
const double epsilon_current = epsilon(element, ion, level);
const double statweight = stat_weight(element, ion, level);
const auto nnlevel = get_levelpop(modelgridindex, element, ion, level);

// Downward transitions within the current ionisation stage:
// radiative/collisional deexcitation and internal downward jumps
double sum_internal_down_same = 0.;
double sum_raddeexc = 0.;
double sum_coldeexc = 0.;
const int ndowntrans = get_ndowntrans(element, ion, level);
const auto *const leveldowntranslist = get_downtranslist(element, ion, level);
auto *const arr_sum_epstrans_rad_deexc = chlevel.sum_epstrans_rad_deexc;
auto *const arr_sum_internal_down_same = chlevel.sum_internal_down_same;
for (int i = 0; i < ndowntrans; i++) {
const auto &downtrans = leveldowntranslist[i];
const int lower = downtrans.targetlevelindex;
const auto A_ul = downtrans.einstein_A;
const double epsilon_target = epsilon(element, ion, lower);
const double epsilon_trans = epsilon_current - epsilon_target;

const double R = rad_deexcitation_ratecoeff(modelgridindex, element, ion, lower, epsilon_trans, A_ul, statweight,
nnlevel, t_mid);
const double C = col_deexcitation_ratecoeff(T_e, nne, epsilon_trans, element, ion, level, downtrans);

sum_raddeexc += R * epsilon_trans;
sum_coldeexc += C * epsilon_trans;
sum_internal_down_same += (R + C) * epsilon_target;

arr_sum_epstrans_rad_deexc[i] = sum_raddeexc;
arr_sum_internal_down_same[i] = sum_internal_down_same;

// printout("checking downtrans %d to level %d: R %g, C %g, epsilon_trans %g\n",i,lower,R,C,epsilon_trans);
}
processrates[MA_ACTION_RADDEEXC] = sum_raddeexc;
processrates[MA_ACTION_COLDEEXC] = sum_coldeexc;
processrates[MA_ACTION_INTERNALDOWNSAME] = sum_internal_down_same;

// Calculate sum for upward internal transitions
// transitions within the current ionisation stage
double sum_internal_up_same = 0.;
const int nuptrans = get_nuptrans(element, ion, level);
const auto *const uptranslist = get_uptranslist(element, ion, level);
for (int i = 0; i < nuptrans; i++) {
const auto &uptrans = uptranslist[i];
const double epsilon_trans = epsilon(element, ion, uptrans.targetlevelindex) - epsilon_current;

const double R = rad_excitation_ratecoeff(modelgridindex, element, ion, level, i, epsilon_trans, nnlevel,
uptrans.lineindex, t_mid);
const double C = col_excitation_ratecoeff(T_e, nne, element, ion, level, i, epsilon_trans, statweight);
const double NT = nonthermal::nt_excitation_ratecoeff(modelgridindex, element, ion, level, i, uptrans.lineindex);

sum_internal_up_same += (R + C + NT) * epsilon_current;
chlevel.sum_internal_up_same[i] = sum_internal_up_same;
}
processrates[MA_ACTION_INTERNALUPSAME] = sum_internal_up_same;

assert_always(std::isfinite(processrates[MA_ACTION_INTERNALUPSAME]));

// Downward transitions to lower ionisation stages:
// radiative/collisional recombination and internal downward jumps
// checks only if there is a lower ion, doesn't make sure that Z(ion)=Z(ion-1)+1
double sum_internal_down_lower = 0.;
double sum_radrecomb = 0.;
double sum_colrecomb = 0.;
if (ion > 0 && level <= get_maxrecombininglevel(element, ion)) {
// nlevels = get_nlevels(element,ion-1);
const int nlevels = get_ionisinglevels(element, ion - 1);
// nlevels = get_ionisinglevels(element,ion-1);
for (int lower = 0; lower < nlevels; lower++) {
const double epsilon_target = epsilon(element, ion - 1, lower);
const double epsilon_trans = epsilon_current - epsilon_target;

const double R = rad_recombination_ratecoeff(T_e, nne, element, ion, level, lower, modelgridindex);
const double C = col_recombination_ratecoeff(modelgridindex, element, ion, level, lower, epsilon_trans);

sum_internal_down_lower += (R + C) * epsilon_target;

sum_radrecomb += R * epsilon_trans;
sum_colrecomb += C * epsilon_trans;
}
}
processrates[MA_ACTION_INTERNALDOWNLOWER] = sum_internal_down_lower;
processrates[MA_ACTION_RADRECOMB] = sum_radrecomb;
processrates[MA_ACTION_COLRECOMB] = sum_colrecomb;

// Transitions to higher ionisation stages
double sum_up_highernt = 0.;
double sum_up_higher = 0.;
const int ionisinglevels = get_ionisinglevels(element, ion);
if (ion < get_nions(element) - 1 && level < ionisinglevels) {
if (NT_ON) {
sum_up_highernt = nonthermal::nt_ionization_ratecoeff(modelgridindex, element, ion) * epsilon_current;
}

const auto nphixstargets = get_nphixstargets(element, ion, level);
for (int phixstargetindex = 0; phixstargetindex < nphixstargets; phixstargetindex++) {
const double epsilon_trans = get_phixs_threshold(element, ion, level, phixstargetindex);

const double R = get_corrphotoioncoeff(element, ion, level, phixstargetindex, modelgridindex);
const double C = col_ionization_ratecoeff(T_e, nne, element, ion, level, phixstargetindex, epsilon_trans);

sum_up_higher += (R + C) * epsilon_current;
}
}
processrates[MA_ACTION_INTERNALUPHIGHERNT] = sum_up_highernt;
processrates[MA_ACTION_INTERNALUPHIGHER] = sum_up_higher;

return processrates;
}

// handle activated macro atoms
__host__ __device__ void do_macroatom(Packet &pkt, const MacroAtomState &pktmastate) {
const int modelgridindex = grid::get_cell_modelgridindex(pkt.where);
Expand Down Expand Up @@ -367,7 +368,7 @@ __host__ __device__ void do_macroatom(Packet &pkt, const MacroAtomState &pktmast
assert_testmodeonly(globals::cellcache[cellcacheslotid].cellnumber == modelgridindex);

// If there are no precalculated rates available then calculate them
if (chlevel.processrates[MA_ACTION_INTERNALUPHIGHER] < 0) {
if (level >= PRECALCULATED_MACROATOM_LEVELS && chlevel.processrates[MA_ACTION_INTERNALUPHIGHER] < 0) {
chlevel.processrates = calculate_macroatom_transitionrates(modelgridindex, element, ion, level, t_mid, chlevel);
}
}
Expand Down
4 changes: 4 additions & 0 deletions macroatom.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ void macroatom_open_file(int my_rank);
void macroatom_close_file();

void do_macroatom(Packet &pkt, const MacroAtomState &pktmastate);
auto calculate_macroatom_transitionrates(int modelgridindex, int element, int ion, int level, double t_mid,
CellCacheLevels &chlevel) -> std::array<double, MA_ACTION_COUNT>;

[[nodiscard]] auto rad_deexcitation_ratecoeff(int modelgridindex, int element, int ion, int lower, double epsilon_trans,
float A_ul, double upperstatweight, double nnlevelupper,
Expand All @@ -31,4 +33,6 @@ void do_macroatom(Packet &pkt, const MacroAtomState &pktmastate);
[[nodiscard]] auto col_excitation_ratecoeff(float T_e, float nne, int element, int ion, int lower, int uptransindex,
double epsilon_trans, double lowerstatweight) -> double;

constexpr int PRECALCULATED_MACROATOM_LEVELS = 20;

#endif // MACROATOM_H
16 changes: 16 additions & 0 deletions update_grid.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
#include "grid.h"
#include "kpkt.h"
#include "ltepop.h"
#include "macroatom.h"
#include "nltepop.h"
#include "nonthermal.h"
#include "radfield.h"
Expand Down Expand Up @@ -1199,6 +1200,7 @@ void cellcache_change_cell(const int modelgridindex) {

cacheslot.cellnumber = modelgridindex;
cacheslot.chi_ff_nnionpart = -1.;
const double t_mid = globals::timesteps[globals::timestep].mid;

const int nelements = get_nelements();
for (int element = 0; element < nelements; element++) {
Expand Down Expand Up @@ -1240,6 +1242,20 @@ void cellcache_change_cell(const int modelgridindex) {
}

if (modelgridindex >= 0) {
for (int element = 0; element < nelements; element++) {
const int nions = get_nions(element);
for (int ion = 0; ion < nions; ion++) {
const int nlevels = get_nlevels(element, ion);
auto &chion = cacheslot.chelements[element].chions[ion];
const auto macroatom_nlevels = std::min(nlevels, PRECALCULATED_MACROATOM_LEVELS);
for (int level = 0; level < macroatom_nlevels; level++) {
auto &chlevel = chion.chlevels[level];
chion.chlevels[level].processrates =
calculate_macroatom_transitionrates(modelgridindex, element, ion, level, t_mid, chlevel);
}
}
}

std::ranges::fill(cacheslot.ch_allcont_departureratios, -1.);

const auto nnetot = grid::get_nnetot(modelgridindex);
Expand Down

0 comments on commit 4f17aa1

Please sign in to comment.