Skip to content

Improve reporting of macro bound-free heating #1186

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

Merged
jhmatthews merged 11 commits into from
Jun 18, 2026
Merged

Improve reporting of macro bound-free heating #1186

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
f1d6a23
some initial preparatory changes for macro_bf reporting
jhmatthews Oct 20, 2025
8eb4ea5
added qrecomb to communicate plasma routines
jhmatthews Oct 21, 2025
344c97f
corrected bug I introduced and added reporting of sums
jhmatthews Oct 23, 2025
5c3de8f
added qrecomb cool to windsave2table
jhmatthews Oct 23, 2025
abb7406
improved macro reporting and new cool_di_macro plasma variable
jhmatthews Oct 30, 2025
0210b6e
added macro energy flow rates
jhmatthews Oct 30, 2025
719e2e0
fixed bug in energy_flow_in
jhmatthews Dec 16, 2025
eae05c0
corrected error in macro-atom energy flow parallelisation (diagnostic…
jhmatthews Jan 7, 2026
1336ae8
fixed lte and test_cooling behaviour for new variables
jhmatthews May 26, 2026
fa68f56
rm unused variable, n
jhmatthews May 26, 2026
b44eb9f
protect macromain access with if statements
jhmatthews May 26, 2026
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
13 changes: 8 additions & 5 deletions source/communicate_macro.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -219,7 +219,7 @@ broadcast_updated_macro_atom_properties (const int n_start, const int n_stop, co

d_xsignal (files.root, "%-20s Begin macro atom updated properties communication\n", "NOK");
const int n_cells_max = get_max_cells_per_rank (NPLASMA);
const int comm_buffer_size = calculate_comm_buffer_size (1 + 3 * n_cells_max, n_cells_max * (6 * size_gamma_est + 2 * size_Jbar_est));
const int comm_buffer_size = calculate_comm_buffer_size (1 + 3 * n_cells_max, n_cells_max * (6 * size_gamma_est + 2 * size_Jbar_est + 1));

char *const comm_buffer = malloc (comm_buffer_size);
if (comm_buffer == NULL)
Expand Down Expand Up @@ -248,6 +248,7 @@ broadcast_updated_macro_atom_properties (const int n_start, const int n_stop, co
MPI_Pack (macromain[n_plasma].alpha_st_old, size_gamma_est, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&macromain[n_plasma].kpkt_rates_known, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&macromain[n_plasma].matrix_rates_known, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&macromain[n_plasma].energy_flow_in, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
}
}

Expand All @@ -271,6 +272,7 @@ broadcast_updated_macro_atom_properties (const int n_start, const int n_stop, co
MPI_Unpack (comm_buffer, comm_buffer_size, &position, macromain[n_plasma].alpha_st_old, size_gamma_est, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &macromain[n_plasma].kpkt_rates_known, 1, MPI_INT, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &macromain[n_plasma].matrix_rates_known, 1, MPI_INT, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &macromain[n_plasma].energy_flow_in, 1, MPI_DOUBLE, MPI_COMM_WORLD);
}
}
}
Expand Down Expand Up @@ -409,15 +411,15 @@ reduce_macro_atom_estimators (void)
gamma_helper = calloc (sizeof (double), NPLASMA * 4 * size_gamma_est);
alpha_helper = calloc (sizeof (double), NPLASMA * 2 * size_alpha_est);
level_helper = calloc (sizeof (double), NPLASMA * nlevels_macro);
cell_helper = calloc (sizeof (double), 8 * NPLASMA);
cell_helper = calloc (sizeof (double), 9 * NPLASMA);
cooling_bf_helper = calloc (sizeof (double), NPLASMA * 2 * nphot_total);
cooling_bb_helper = calloc (sizeof (double), NPLASMA * nlines);

jbar_helper2 = calloc (sizeof (double), NPLASMA * size_Jbar_est);
gamma_helper2 = calloc (sizeof (double), NPLASMA * 4 * size_gamma_est);
alpha_helper2 = calloc (sizeof (double), NPLASMA * 2 * size_alpha_est);
level_helper2 = calloc (sizeof (double), NPLASMA * nlevels_macro);
cell_helper2 = calloc (sizeof (double), 8 * NPLASMA);
cell_helper2 = calloc (sizeof (double), 9 * NPLASMA);
cooling_bf_helper2 = calloc (sizeof (double), NPLASMA * 2 * nphot_total);
cooling_bb_helper2 = calloc (sizeof (double), NPLASMA * nlines);

Expand All @@ -436,7 +438,7 @@ reduce_macro_atom_estimators (void)
cell_helper[mpi_i + 5 * NPLASMA] = macromain[mpi_i].cooling_ff / np_mpi_global;
cell_helper[mpi_i + 6 * NPLASMA] = macromain[mpi_i].cooling_ff_lofreq / np_mpi_global;
cell_helper[mpi_i + 7 * NPLASMA] = macromain[mpi_i].cooling_adiabatic / np_mpi_global;

cell_helper[mpi_i + 8 * NPLASMA] = macromain[mpi_i].energy_flow_out / np_mpi_global;


for (n = 0; n < nlevels_macro; n++)
Expand Down Expand Up @@ -477,7 +479,7 @@ reduce_macro_atom_estimators (void)

/* because in the above loop we have already divided by number of processes, we can now do a sum
with MPI_Reduce, passing it MPI_SUM as an argument. This will give us the mean across threads */
MPI_Allreduce (cell_helper, cell_helper2, NPLASMA * 8, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (cell_helper, cell_helper2, NPLASMA * 9, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (level_helper, level_helper2, NPLASMA * nlevels_macro, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (jbar_helper, jbar_helper2, NPLASMA * size_Jbar_est, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
MPI_Allreduce (gamma_helper, gamma_helper2, NPLASMA * 4 * size_gamma_est, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
Expand All @@ -501,6 +503,7 @@ reduce_macro_atom_estimators (void)
macromain[mpi_i].cooling_ff = cell_helper2[mpi_i + 5 * NPLASMA];
macromain[mpi_i].cooling_ff_lofreq = cell_helper2[mpi_i + 6 * NPLASMA];
macromain[mpi_i].cooling_adiabatic = cell_helper2[mpi_i + 7 * NPLASMA];
macromain[mpi_i].energy_flow_out = cell_helper2[mpi_i + 8 * NPLASMA];


for (n = 0; n < nlevels_macro; n++)
Expand Down
14 changes: 11 additions & 3 deletions source/communicate_plasma.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
/* these next two numbers are the number of basic doubles and integers
in the plasma structure to be packed (i.e. things that are not arrays or pointers and
have a size of MPI_INT or MPI_DOUBLE)*/
#define N_BASIC_DOUBLES 73
#define N_BASIC_DOUBLES 75
#define N_BASIC_INTS 22

/**********************************************************/
Expand Down Expand Up @@ -113,8 +113,10 @@ broadcast_plasma_grid (const int n_start, const int n_stop, const int n_cells_ra
MPI_Pack (&cell->heat_ind_comp, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->heat_lines_macro, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->heat_photo_macro, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->heat_qrecomb_macro, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->cool_lines_macro, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->cool_bf_macro, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->cool_di_macro, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->heat_photo, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->heat_z, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
MPI_Pack (&cell->heat_auger, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD);
Expand Down Expand Up @@ -272,8 +274,10 @@ broadcast_plasma_grid (const int n_start, const int n_stop, const int n_cells_ra
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->heat_ind_comp, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->heat_lines_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->heat_photo_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->heat_qrecomb_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->cool_lines_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->cool_bf_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->cool_di_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->heat_photo, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->heat_z, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, comm_buffer_size, &position, &cell->heat_auger, 1, MPI_DOUBLE, MPI_COMM_WORLD);
Expand Down Expand Up @@ -615,8 +619,8 @@ broadcast_updated_plasma_properties (const int n_start_rank, const int n_stop_ra
//OLD const int num_ints = 1 + n_cells_max * (20 + nphot_total + 2 * NXBANDS + 2 * N_PHOT_PROC + nions);
const int num_ints = 1 + n_cells_max * (N_BASIC_INTS + nphot_total + 2 * NXBANDS + 2 * N_PHOT_PROC + nions);
const int num_doubles =
n_cells_max * (N_BASIC_DOUBLES + 1 * 3 + 9 * 4 + 6 * NFLUX_ANGLES + 3 * NFORCE_DIRECTIONS + 9 * nions + 1 * nlte_levels + 3 * nphot_total +
1 * n_inner_tot + 9 * NXBANDS + 1 * NBINS_IN_CELL_SPEC);
n_cells_max * (N_BASIC_DOUBLES + 1 * 3 + 9 * 4 + 6 * NFLUX_ANGLES + 3 * NFORCE_DIRECTIONS + 9 * nions + 1 * nlte_levels +
3 * nphot_total + 1 * n_inner_tot + 9 * NXBANDS + 1 * NBINS_IN_CELL_SPEC);

const int size_of_comm_buffer = calculate_comm_buffer_size (num_ints, num_doubles);
char *const comm_buffer = malloc (size_of_comm_buffer);
Expand Down Expand Up @@ -669,8 +673,10 @@ broadcast_updated_plasma_properties (const int n_start_rank, const int n_stop_ra
MPI_Pack (&plasmamain[n_plasma].heat_ind_comp, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].heat_lines_macro, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].heat_photo_macro, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].heat_qrecomb_macro, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].cool_lines_macro, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].cool_bf_macro, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].cool_di_macro, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].heat_photo, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].heat_z, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n_plasma].heat_auger, 1, MPI_DOUBLE, comm_buffer, size_of_comm_buffer, &position, MPI_COMM_WORLD);
Expand Down Expand Up @@ -841,8 +847,10 @@ broadcast_updated_plasma_properties (const int n_start_rank, const int n_stop_ra
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].heat_ind_comp, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].heat_lines_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].heat_photo_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].heat_qrecomb_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].cool_lines_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].cool_bf_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].cool_di_macro, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].heat_photo, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].heat_z, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (comm_buffer, size_of_comm_buffer, &position, &plasmamain[n_plasma].heat_auger, 1, MPI_DOUBLE, MPI_COMM_WORLD);
Expand Down
18 changes: 10 additions & 8 deletions source/cooling.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,14 +147,16 @@ xtotal_emission (one, f1, f2)
if (geo.rt_mode == RT_MODE_MACRO) //Switch for macro atoms (SS)
{
xplasma->cool_bf_macro = total_fb_matoms (xplasma, t_e, f1, f2);
xplasma->cool_rr = xplasma->cool_bf_macro + total_fb (xplasma, t_e, f1, f2, FB_REDUCED, OUTER_SHELL); //outer shellrecombinations
//The first term here is the fb cooling due to macro ions and the second gives
//the fb cooling due to simple ions.
//total_fb has been modified to exclude recombinations treated using macro atoms.
//Note: This the fb_matom call makes no use of f1 or f2. They are passed for
//now in case they should be used in the future. But they could
//also be removed.
// (SS)
/* include macro collisional ionization cooling (this used to be part of total_fb_matoms) */
xplasma->cool_di_macro = cooling_di_matoms (xplasma, t_e, f1, f2);


xplasma->cool_rr = xplasma->cool_bf_macro + xplasma->cool_di_macro + total_fb (xplasma, t_e, f1, f2, FB_REDUCED, OUTER_SHELL); //outer shellrecombinations
/* The sfurst term here is the fb cooling due to macro ions and the thirdgives
the fb cooling due to simple ions. total_fb has been modified to exclude recombinations treated using macro atoms.
Note: This the fb_matom call makes no use of f1 or f2. They are passed for now in case they should be used in the future.
But they could also be removed. (SS) */

cooling = xplasma->cool_rr;
xplasma->lum_lines = total_bb_cooling (xplasma, t_e);
cooling += xplasma->lum_lines;
Expand Down
14 changes: 9 additions & 5 deletions source/emission.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,23 +189,27 @@ total_emission (xplasma, f1, f2)
{
if (geo.rt_mode == RT_MODE_MACRO)
{
xplasma->lum_rr = (total_fb_matoms (xplasma, t_e, f1, f2) + total_fb (xplasma, t_e, f1, f2, FB_FULL, OUTER_SHELL)); //outer shellrecombinations

/*
*The first term here is the fb cooling due to macro ions and the second gives
*the fb cooling due to simple ions.
*total_fb has been modified to exclude recombinations treated using macro atoms.
*/
xplasma->lum_rr = (total_fb_matoms (xplasma, t_e, f1, f2) + total_fb (xplasma, t_e, f1, f2, FB_FULL, OUTER_SHELL)); //outer shellrecombinations
/* historically we include macro collisional ionization cooling (this used to be part of total_fb_matoms) */
xplasma->lum_rr += cooling_di_matoms (xplasma, t_e, f1, f2);


xplasma->lum_tot = xplasma->cool_rr;
/* Note: This the fb_matom call makes no use of f1 or f2. They are passed for
/* Note: the fb_matom call makes no use of f1 or f2. They are passed for
* now in case they should be used in the future. But they could
* also be removed.
* (SS)
* also be removed. (SS)
*/

xplasma->lum_lines = total_bb_cooling (xplasma, t_e);
xplasma->lum_tot += xplasma->lum_lines;
/* total_bb_cooling gives the total cooling rate due to bb transisions whether they
are macro atoms or simple ions. */

xplasma->lum_ff = total_free (xplasma, t_e, f1, f2);
xplasma->lum_tot += xplasma->lum_ff;

Expand Down
Loading
Loading