diff --git a/source/communicate_macro.c b/source/communicate_macro.c index 573404861..0b49aaf6c 100644 --- a/source/communicate_macro.c +++ b/source/communicate_macro.c @@ -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) @@ -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 (¯omain[n_plasma].kpkt_rates_known, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD); MPI_Pack (¯omain[n_plasma].matrix_rates_known, 1, MPI_INT, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD); + MPI_Pack (¯omain[n_plasma].energy_flow_in, 1, MPI_DOUBLE, comm_buffer, comm_buffer_size, &position, MPI_COMM_WORLD); } } @@ -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, ¯omain[n_plasma].kpkt_rates_known, 1, MPI_INT, MPI_COMM_WORLD); MPI_Unpack (comm_buffer, comm_buffer_size, &position, ¯omain[n_plasma].matrix_rates_known, 1, MPI_INT, MPI_COMM_WORLD); + MPI_Unpack (comm_buffer, comm_buffer_size, &position, ¯omain[n_plasma].energy_flow_in, 1, MPI_DOUBLE, MPI_COMM_WORLD); } } } @@ -409,7 +411,7 @@ 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); @@ -417,7 +419,7 @@ reduce_macro_atom_estimators (void) 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); @@ -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++) @@ -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); @@ -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++) diff --git a/source/communicate_plasma.c b/source/communicate_plasma.c index a23a4f5ba..3f46c9e88 100644 --- a/source/communicate_plasma.c +++ b/source/communicate_plasma.c @@ -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 /**********************************************************/ @@ -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); @@ -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); @@ -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); @@ -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); @@ -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); diff --git a/source/cooling.c b/source/cooling.c index b05797b6e..78a099f01 100644 --- a/source/cooling.c +++ b/source/cooling.c @@ -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; diff --git a/source/emission.c b/source/emission.c index 36a502510..855fb8608 100644 --- a/source/emission.c +++ b/source/emission.c @@ -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; diff --git a/source/estimators_macro.c b/source/estimators_macro.c index 9a4554865..21a5d02c3 100644 --- a/source/estimators_macro.c +++ b/source/estimators_macro.c @@ -536,13 +536,19 @@ normalise_macro_estimators (PlasmaPtr xplasma) xplasma->heat_lines += heat_contribution = macro_bb_heating (xplasma, xplasma->t_e); xplasma->heat_lines_macro = heat_contribution; - xplasma->heat_tot += heat_contribution; /* Get the bf heating contributions here too. (SS June 04) */ + /* JM Oct 2025 -- we now separate out three body recombination and photoionization heating */ + /* so that we can track them separately */ - xplasma->heat_photo += heat_contribution = macro_bf_heating (xplasma, xplasma->t_e); + xplasma->heat_photo += heat_contribution = macro_photo_heating (xplasma, xplasma->t_e); xplasma->heat_photo_macro = heat_contribution; - xplasma->heat_tot += heat_contribution; + + xplasma->heat_photo += heat_contribution = macro_qrecomb_heating (xplasma, xplasma->t_e); + xplasma->heat_qrecomb_macro = heat_contribution; + + /* ensure the total heating is incremented too, for line and bound-free */ + xplasma->heat_tot += xplasma->heat_lines_macro + xplasma->heat_photo_macro + xplasma->heat_qrecomb_macro; /* finally, check if we have any places where stimulated recombination wins over photoionization */ @@ -557,6 +563,13 @@ normalise_macro_estimators (PlasmaPtr xplasma) mplasma->kpkt_rates_known = FALSE; mplasma->matrix_rates_known = FALSE; + /* record a total energy flow into macro-atoms, by summing over all matom_abs contributions */ + mplasma->energy_flow_in = 0.0; + for (i = 0; i < nlevels_macro; i++) + { + mplasma->energy_flow_in += mplasma->matom_abs[i]; + } + return (0); } @@ -597,7 +610,6 @@ total_fb_matoms (xplasma, t_e, f1, f2) struct topbase_phot *cont_ptr; double total, density; int i, j; - double q_ioniz (); MacroPtr mplasma; mplasma = ¯omain[xplasma->nplasma]; @@ -630,23 +642,69 @@ total_fb_matoms (xplasma, t_e, f1, f2) - mplasma->alpha_st_old[xconfig[i].bfu_indx_first + j] - alpha_sp (cont_ptr, xplasma, 0)) * PLANCK * phot_top[xconfig[i].bfu_jump[j]].freq[0] * density * xplasma->ne * xplasma->vol; + total += cool_contribution; + } + } + + xplasma->t_e = t_e_store; //restore the original value + + } + + return (total); +} + + +/**********************************************************/ +/** + * @brief computes the cooling rate due to collisional ionization in macro-atoms + * + * @param [in] PlasmaPtr xplasma + * @param [in] double t_e electron temperature + * @param [in] double f1 lower frequency + * @param [in] double f2 upper frequency + * @return total + * + * @details + * computes the cooling rate due to collisional ionization in macro-atoms + * + **********************************************************/ + +double +cooling_di_matoms (xplasma, t_e, f1, f2) + PlasmaPtr xplasma; + double t_e; + double f1, f2; +{ + double cool_contribution; + struct topbase_phot *cont_ptr; + double total, density; + int i, j; + + total = 0; + + if (geo.macro_simple == FALSE) //allow for "only-simple" calculations (SS May04) + { + for (i = 0; i < nlte_levels; i++) + { + for (j = 0; j < xconfig[i].n_bfu_jump; j++) + { + /* Need the density for the lower level in the ionization process. */ + cont_ptr = &phot_top[xconfig[i].bfu_jump[j]]; + /* Now add the collisional ionization term. */ density = den_config (xplasma, cont_ptr->nlev); - cool_contribution += + cool_contribution = q_ioniz (cont_ptr, t_e) * density * xplasma->ne * PLANCK * phot_top[xconfig[i].bfu_jump[j]].freq[0] * xplasma->vol; - /* That's the bf cooling contribution. */ total += cool_contribution; } } - - xplasma->t_e = t_e_store; //restore the original value - } return (total); } + /**********************************************************/ /** * @brief computes the total cooling in matom bb transitions @@ -763,27 +821,25 @@ macro_bb_heating (xplasma, t_e) /**********************************************************/ /** - * @brief computes the total heating due to bf transitions for macro atoms. + * @brief computes the total heating due to photoionization for macro atoms. * * @param [in] PlasmaPtr xplasma * @param [in] double t_e electron temperature * @return total * * @details - * computes the total heating due to bf transitions for macro atoms. - * The heating in simple ions - * is taken care of elsewhere. + * computes the total heating due to bf transitions for macro atoms, not including the three body recombination part. + * The heating in simple ions is taken care of elsewhere. * It is used by the heating/cooling calculation to get the temperature. - * **********************************************************/ double -macro_bf_heating (xplasma, t_e) +macro_photo_heating (xplasma, t_e) PlasmaPtr xplasma; double t_e; { double heat_contribution; - double total, upper_density, lower_density; + double total, lower_density; int i, j; double q_recomb (); MacroPtr mplasma; @@ -804,6 +860,44 @@ macro_bf_heating (xplasma, t_e) mplasma->gamma_old[xconfig[i].bfu_indx_first + j]) * PLANCK * phot_top[xconfig[i].bfu_jump[j]].freq[0] * lower_density * xplasma->vol; + total += heat_contribution; + } + } + + return (total); +} + + +/**********************************************************/ +/** + * @brief computes the total heating due to threebody recombination bf transitions for macro atoms. + * + * @param [in] PlasmaPtr xplasma + * @param [in] double t_e electron temperature + * @return total + * + * @details + * computes the total heating due to the three body recombination part for macro atoms. + * The heating in simple ions is taken care of elsewhere. + * It is used by the heating/cooling calculation to get the temperature. + **********************************************************/ + +double +macro_qrecomb_heating (xplasma, t_e) + PlasmaPtr xplasma; + double t_e; +{ + double heat_contribution; + double total, upper_density; + int i, j; + + total = 0; // initialise + + for (i = 0; i < nlte_levels; i++) + { + for (j = 0; j < xconfig[i].n_bfu_jump; j++) + { + heat_contribution = 0.0; /* Three body recombination part. */ upper_density = den_config (xplasma, phot_top[xconfig[i].bfu_jump[j]].uplev); heat_contribution += diff --git a/source/ionization.c b/source/ionization.c index 5d4ed9691..6f1ee0b24 100644 --- a/source/ionization.c +++ b/source/ionization.c @@ -474,6 +474,7 @@ static double lte_heat_photo_orig, lte_heat_ff_orig; static double lte_heat_comp_orig, lte_heat_ind_comp_orig; static double lte_heat_lines_orig, lte_heat_auger_orig; static double lte_heat_lines_macro_orig, lte_heat_photo_macro_orig; +static double lte_heat_qrecomb_macro_orig; static double lte_heat_ch_ex_orig; @@ -641,15 +642,21 @@ calc_te (PlasmaPtr xplasma, double tmin, double tmax) xplasma->heat_tot += xplasma->heat_lines_macro; xplasma->heat_lines += xplasma->heat_lines_macro; - /* Similaryly for macro_atom_bf_heating */ + /* Similarly for macro_atom_bf_heating, for both photoionization + and three body recombination components */ xplasma->heat_tot -= xplasma->heat_photo_macro; xplasma->heat_photo -= xplasma->heat_photo_macro; + xplasma->heat_tot -= xplasma->heat_qrecomb_macro; + xplasma->heat_photo -= xplasma->heat_qrecomb_macro; - xplasma->heat_photo_macro = macro_bf_heating (xplasma, xplasma->t_e); + xplasma->heat_photo_macro = macro_photo_heating (xplasma, xplasma->t_e); + xplasma->heat_qrecomb_macro = macro_qrecomb_heating (xplasma, xplasma->t_e); xplasma->heat_tot += xplasma->heat_photo_macro; xplasma->heat_photo += xplasma->heat_photo_macro; + xplasma->heat_tot += xplasma->heat_qrecomb_macro; + xplasma->heat_photo += xplasma->heat_qrecomb_macro; return (xplasma->t_e); @@ -702,7 +709,7 @@ zero_emit (double t) { double difference; - /*Original method */ + /* Original method */ xxxplasma->t_e = t; @@ -715,23 +722,29 @@ zero_emit (double t) xxxplasma->heat_tot += xxxplasma->heat_lines_macro; xxxplasma->heat_lines += xxxplasma->heat_lines_macro; + /* Similarly for macro_atom_bf_heating, recompute for both photoionization + and three body recombination components */ + xxxplasma->heat_tot -= xxxplasma->heat_photo_macro; xxxplasma->heat_photo -= xxxplasma->heat_photo_macro; + xxxplasma->heat_tot -= xxxplasma->heat_qrecomb_macro; + xxxplasma->heat_photo -= xxxplasma->heat_qrecomb_macro; - xxxplasma->heat_photo_macro = macro_bf_heating (xxxplasma, t); + xxxplasma->heat_photo_macro = macro_photo_heating (xxxplasma, t); + xxxplasma->heat_qrecomb_macro = macro_qrecomb_heating (xxxplasma, t); xxxplasma->heat_tot += xxxplasma->heat_photo_macro; xxxplasma->heat_photo += xxxplasma->heat_photo_macro; + xxxplasma->heat_tot += xxxplasma->heat_qrecomb_macro; + xxxplasma->heat_photo += xxxplasma->heat_qrecomb_macro; - /* Finished macro atom corrections */ + /* Finished macro atom corrections. Now compute the total cooling and heating - cooling */ cooling (xxxplasma, t); difference = xxxplasma->heat_tot + xxxplasma->heat_shock - xxxplasma->cool_tot; - - return (difference); } @@ -848,12 +861,16 @@ zero_emit_lte (double t) /* Now subtract the original macro contributions (which are included in the scaled heat_photo and heat_lines above) and replace with freshly computed macro heating at the new temperature and densities */ + /* note that we have to include collisional recombination (qrecomb) */ xxxplasma->heat_photo -= lte_heat_photo_macro_orig; + xxxplasma->heat_photo -= lte_heat_qrecomb_macro_orig; - xxxplasma->heat_photo_macro = macro_bf_heating (xxxplasma, t); - + xxxplasma->heat_photo_macro = macro_photo_heating (xxxplasma, t); xxxplasma->heat_photo += xxxplasma->heat_photo_macro; + xxxplasma->heat_qrecomb_macro = macro_qrecomb_heating (xxxplasma, t); + xxxplasma->heat_photo += xxxplasma->heat_qrecomb_macro; + xxxplasma->heat_lines_macro = macro_bb_heating (xxxplasma, t); xxxplasma->heat_lines += xxxplasma->heat_lines_macro; @@ -939,6 +956,7 @@ calc_te_lte (PlasmaPtr xplasma, double tmin, double tmax) lte_heat_auger_orig = xplasma->heat_auger; lte_heat_lines_macro_orig = xplasma->heat_lines_macro; lte_heat_photo_macro_orig = xplasma->heat_photo_macro; + lte_heat_qrecomb_macro_orig = xplasma->heat_qrecomb_macro; lte_heat_ch_ex_orig = xplasma->heat_ch_ex; /* Evaluate heating-cooling difference at bracket endpoints */ diff --git a/source/log.h b/source/log.h index 241077e0d..05784060b 100644 --- a/source/log.h +++ b/source/log.h @@ -3,6 +3,7 @@ int opar(char filename[]); int add_par(char filename[]); int cpar(char filename[]); int rdpar_init(void); +char *check_and_fix_string(char *s); int string_process(char question[], char dummy[]); int string_process_from_command_line(char question[], char dummy[]); int string_process_from_file(char question[], char dummy[]); @@ -23,7 +24,6 @@ int rdpar_set_mpi_rank(int rank); int rdpar_set_verbose(int vlevel); int rdpar_check(void); /* xlog.c */ -void print_memory_usage(const char *label); int Log_init(char *filename); int Log_append(char *filename); void Log_close(void); diff --git a/source/macro_gov.c b/source/macro_gov.c index 5a98740cd..947ca81d6 100644 --- a/source/macro_gov.c +++ b/source/macro_gov.c @@ -185,6 +185,7 @@ macro_gov (p, nres, matom_or_kpkt, which_out) { emit_matom (wmain, p, nres, new_uplvl, 0, VERY_BIG); *which_out = MATOM; + mplasma->energy_flow_out += p->w; } if (p->origin < 10) @@ -214,6 +215,7 @@ macro_gov (p, nres, matom_or_kpkt, which_out) /* It escapes as a r-packet that was created by de-activation of a macro atom. */ *which_out = MATOM; + mplasma->energy_flow_out += p->w; /* Update the the photon origin to indicate the packet has been processed by a macro atom */ @@ -233,6 +235,7 @@ macro_gov (p, nres, matom_or_kpkt, which_out) /* It escapes as a r-packet that was created by de-activation of a macro atom. */ *which_out = MATOM; + mplasma->energy_flow_out += p->w; /* Update the the photon origin to indicate the packet has been processed by a macro atom */ if (p->origin < 10) @@ -266,6 +269,7 @@ macro_gov (p, nres, matom_or_kpkt, which_out) destroyed by collisionally exciting a macro atom so... */ matom_or_kpkt = MATOM; + mplasma->energy_flow_out += p->w; } else { diff --git a/source/sirocco.h b/source/sirocco.h index d52d35751..61230f460 100644 --- a/source/sirocco.h +++ b/source/sirocco.h @@ -938,10 +938,11 @@ typedef struct plasma double heat_lines, heat_ff; double heat_comp; /**< The compton heating for the cell */ double heat_ind_comp; /**< The induced compton heatingfor the cell */ - double heat_lines_macro, heat_photo_macro; /**< bb and bf heating due to macro atoms. Subset of heat_lines - and heat_photo. SS June 04. */ - double cool_lines_macro, cool_bf_macro; /**< bb and bf cooling due to macro atoms. Subset of heat_lines - and heat_photo. SS June 04. */ + double heat_lines_macro, heat_photo_macro; /**< bb and bf heating due to macro atoms. Subset of heat_lines */ + + double cool_lines_macro, cool_bf_macro; /* bb and bf cooling due to macro atoms. */ + double heat_qrecomb_macro; /**< The heating due to macro atom 3-body recombinations, added Oct 2025 */ + double cool_di_macro; /**< direct/collisional ionization cooling due to macro atoms */ double heat_photo, heat_z; /**< photoionization heating total and of metals */ double heat_auger; /**< photoionization heating due to inner shell ionizations */ double heat_ch_ex; @@ -1210,6 +1211,9 @@ typedef struct macro double *matom_abs; /**< This is the energy absorbed by the macro atom levels - recorded during the ionization cycles and used to get matom_emiss (SS) */ + + double energy_flow_out, energy_flow_in; /**< Total energy flow in/out of macro atom levels in cell during ionization cycles */ + /* This portion of the macro structure is not written out by windsave */ int kpkt_rates_known; diff --git a/source/templates.h b/source/templates.h index a6c3d5c6f..0cb9c26d3 100644 --- a/source/templates.h +++ b/source/templates.h @@ -186,9 +186,11 @@ int bf_estimators_increment(WindPtr one, PhotPtr p, double ds); int bb_estimators_increment(WindPtr one, PhotPtr p, double tau_sobolev, double dvds, int nn); int normalise_macro_estimators(PlasmaPtr xplasma); double total_fb_matoms(PlasmaPtr xplasma, double t_e, double f1, double f2); +double cooling_di_matoms(PlasmaPtr xplasma, double t_e, double f1, double f2); double total_bb_cooling(PlasmaPtr xplasma, double t_e); double macro_bb_heating(PlasmaPtr xplasma, double t_e); -double macro_bf_heating(PlasmaPtr xplasma, double t_e); +double macro_photo_heating(PlasmaPtr xplasma, double t_e); +double macro_qrecomb_heating(PlasmaPtr xplasma, double t_e); int bb_simple_heat(PlasmaPtr xplasma, PhotPtr p, double tau_sobolev, int nn); int check_stimulated_recomb(PlasmaPtr xplasma); int get_dilute_estimators(PlasmaPtr xplasma); @@ -459,6 +461,7 @@ int opar(char filename[]); int add_par(char filename[]); int cpar(char filename[]); int rdpar_init(void); +char *check_and_fix_string(char *s); int string_process(char question[], char dummy[]); int string_process_from_command_line(char question[], char dummy[]); int string_process_from_file(char question[], char dummy[]); @@ -711,7 +714,6 @@ int create_spec_table(int ndom, char rootname[]); int create_detailed_cell_spec_table(int ncell, char rootname[]); int create_big_detailed_spec_table(int ndom, char *rootname); /* xlog.c */ -void print_memory_usage(const char *label); int Log_init(char *filename); int Log_append(char *filename); void Log_close(void); diff --git a/source/test_cooling.c b/source/test_cooling.c index 4ec85a367..41dbe25c4 100644 --- a/source/test_cooling.c +++ b/source/test_cooling.c @@ -227,11 +227,16 @@ xcalc_te (PlasmaPtr xplasma, double tmin, double tmax) xplasma->heat_tot -= xplasma->heat_photo_macro; xplasma->heat_photo -= xplasma->heat_photo_macro; + xplasma->heat_tot -= xplasma->heat_qrecomb_macro; + xplasma->heat_photo -= xplasma->heat_qrecomb_macro; - xplasma->heat_photo_macro = macro_bf_heating (xplasma, xplasma->t_e); + xplasma->heat_photo_macro = macro_photo_heating (xplasma, xplasma->t_e); + xplasma->heat_qrecomb_macro = macro_qrecomb_heating (xplasma, xplasma->t_e); xplasma->heat_tot += xplasma->heat_photo_macro; xplasma->heat_photo += xplasma->heat_photo_macro; + xplasma->heat_tot += xplasma->heat_qrecomb_macro; + xplasma->heat_photo += xplasma->heat_qrecomb_macro; return (xplasma->t_e); diff --git a/source/wind_updates2d.c b/source/wind_updates2d.c index 2e38cae9d..5317a091b 100644 --- a/source/wind_updates2d.c +++ b/source/wind_updates2d.c @@ -55,6 +55,9 @@ wind_update (WindPtr w) double xsum, psum, fsum, lsum, csum, icsum, ausum, chexsum; double cool_sum, lum_sum, radiated_luminosity_sum; //1706 - the total cooling and luminosity of the wind double apsum, aausum, abstot; //Absorbed photon energy from PI and auger + double heat_macro_photo_sum, heat_macro_qrecomb_sum, heat_macro_lines_sum; + double cool_macro_photo_sum, cool_macro_di_sum, cool_macro_lines_sum; + double macro_energy_in, macro_energy_out; double flux_persist_scale; double volume; double dt_r, dt_e; @@ -149,7 +152,6 @@ wind_update (WindPtr w) } /* Calculate the densities in various ways depending on the ioniz_mode */ - Log ("XXXXX updating %3d %3d\n", geo.wcycle, n_plasma); ion_abundances (&plasmamain[n_plasma], geo.ioniz_mode); } @@ -200,7 +202,6 @@ wind_update (WindPtr w) } } /* Finished updating region outside of wind */ - /* Check the balance between the absorbed and the emitted flux */ /* NSH 0717 - ensure the cooling and luminosities reflect the current temperature */ @@ -221,6 +222,9 @@ wind_update (WindPtr w) aausum = 0.0; abstot = 0.0; chexsum = 0.0; + heat_macro_photo_sum = heat_macro_qrecomb_sum = heat_macro_lines_sum = 0.0; + cool_macro_photo_sum = cool_macro_di_sum = cool_macro_lines_sum = 0.0; + macro_energy_in = macro_energy_out = 0.0; /* Each rank now has updated plasma cells (temperature, ion abundances, heat/cool rates, etc.), so we can now find * out what the max d_t is in the wind and also sum up properties to find the total/global values */ @@ -269,6 +273,19 @@ wind_update (WindPtr w) apsum += plasmamain[n_plasma].abs_photo; aausum += plasmamain[n_plasma].abs_auger; chexsum += plasmamain[n_plasma].heat_ch_ex; + + heat_macro_photo_sum += plasmamain[n_plasma].heat_photo_macro; + heat_macro_qrecomb_sum += plasmamain[n_plasma].heat_qrecomb_macro; + heat_macro_lines_sum += plasmamain[n_plasma].heat_lines_macro; + cool_macro_photo_sum += plasmamain[n_plasma].cool_bf_macro; + cool_macro_lines_sum += plasmamain[n_plasma].cool_lines_macro; + cool_macro_di_sum += plasmamain[n_plasma].cool_di_macro; + + if (geo.rt_mode == RT_MODE_MACRO) + { + macro_energy_in += macromain[n_plasma].energy_flow_in; + macro_energy_out += macromain[n_plasma].energy_flow_out; + } } /* We can now calculate the average of the t */ @@ -332,11 +349,15 @@ wind_update (WindPtr w) ("!!wind_update: Wind cooling %8.2e (recomb %8.2e ff %8.2e compton %8.2e DR %8.2e DI %8.2e lines %8.2e adiabatic %8.2e) after update\n", cool_sum, geo.cool_rr, geo.lum_ff, geo.cool_comp, geo.cool_dr, geo.cool_di, geo.lum_lines, geo.cool_adiabatic); - if (modes.use_upweighting_of_simple_macro_atoms) + if (geo.rt_mode == RT_MODE_MACRO) { - /* If we have "indivisible packet" mode on but are using the - upweighting scheme for simple atoms then we report the flows into and out of the ion pool */ - if (geo.rt_mode == RT_MODE_MACRO) + Log ("!!wind_update: macro-atom heating: photoionization %8.2e three body recomb %8.2e lines %8.2e\n", heat_macro_photo_sum, + heat_macro_qrecomb_sum, heat_macro_lines_sum); + Log ("!!wind_update: macro-atom cooling: photoionization %8.2e collisional ionization %8.2e lines %8.2e\n", cool_macro_photo_sum, + cool_macro_di_sum, cool_macro_lines_sum); + Log ("!!wind_update: macro-atom energy flow: in %8.2e out %8.2e\n", macro_energy_in, macro_energy_out); + + if (modes.use_upweighting_of_simple_macro_atoms) { report_bf_simple_ionpool (); } @@ -717,6 +738,13 @@ init_macro_rad_properties (void) plasmamain[n_plasma].kpkt_emiss = 0.0; plasmamain[n_plasma].kpkt_abs = 0.0; + + if (geo.rt_mode == RT_MODE_MACRO) /* macromain is only allocated if geo.rt_mode == RT_MODE_MACRO */ + { + macromain[n_plasma].energy_flow_out = 0.0; + macromain[n_plasma].energy_flow_in = 0.0; + } + for (macro_level = 0; macro_level < nlevels_macro; ++macro_level) { macromain[n_plasma].matom_abs[macro_level] = 0.0; @@ -738,7 +766,6 @@ init_macro_rad_properties (void) /* calculating recomb_sp and recomb_simple is expensive due to calls to * `alpha_sp()` , so we do this part of the initialisation in parallel */ - #ifdef MPI_ON n_cells = get_parallel_nrange (rank_global, NPLASMA, np_mpi_global, &n_start, &n_stop); #else diff --git a/source/windsave2table_sub.c b/source/windsave2table_sub.c index b499f4134..ae7572584 100644 --- a/source/windsave2table_sub.c +++ b/source/windsave2table_sub.c @@ -521,15 +521,26 @@ create_heat_table (ndom, rootname) c[27] = get_one (ndom, "heat_photo_macro"); strcpy (column_name[27], "ht_ph_macro"); + c[28] = get_one (ndom, "heat_qrecomb_macro"); + strcpy (column_name[28], "ht_qr_macro"); - c[28] = get_one (ndom, "cool_lines_macro"); - strcpy (column_name[28], "cl_ln_macro"); + c[29] = get_one (ndom, "cool_lines_macro"); + strcpy (column_name[29], "cl_ln_macro"); - c[29] = get_one (ndom, "cool_bf_macro"); - strcpy (column_name[29], "cl_ph_macro"); + c[30] = get_one (ndom, "cool_bf_macro"); + strcpy (column_name[30], "cl_ph_macro"); + + c[31] = get_one (ndom, "cool_di_macro"); + strcpy (column_name[31], "cl_di_macro"); + + c[32] = get_one (ndom, "energy_in_macro"); + strcpy (column_name[32], "e_in_macro"); + + c[33] = get_one (ndom, "energy_out_macro"); + strcpy (column_name[33], "e_out_macro"); /* This should be the maximum number above +1 */ - ncols = 30; + ncols = 34; converge = get_one (ndom, "converge"); @@ -1536,6 +1547,10 @@ get_one (ndom, variable_name) { x[n] = plasmamain[nplasma].heat_photo_macro; } + else if (strcmp (variable_name, "heat_qrecomb_macro") == 0) + { + x[n] = plasmamain[nplasma].heat_qrecomb_macro; + } else if (strcmp (variable_name, "cool_lines_macro") == 0) { x[n] = plasmamain[nplasma].cool_lines_macro; @@ -1544,7 +1559,18 @@ get_one (ndom, variable_name) { x[n] = plasmamain[nplasma].cool_bf_macro; } - + else if (strcmp (variable_name, "cool_di_macro") == 0) + { + x[n] = plasmamain[nplasma].cool_di_macro; + } + else if (strcmp (variable_name, "energy_in_macro") == 0) + { + x[n] = macromain[nplasma].energy_flow_in; + } + else if (strcmp (variable_name, "energy_out_macro") == 0) + { + x[n] = macromain[nplasma].energy_flow_out; + } else if (strcmp (variable_name, "gain") == 0) { x[n] = plasmamain[nplasma].gain;