Skip to content

Commit

Permalink
Add diagnostics for ion pool for simple atoms in macro mode(#1034)
Browse files Browse the repository at this point in the history 
* Provides additional diagnostics of energy flow into and out of the ion pool (see #1024)
* Version changed to 87f due to the fact that there are additional variables in the plasma structure
* There should be no actual changes to the way photons are processed.
  • Loading branch information
@kslong
kslong authored Oct 16, 2023
1 parent 60be552 commit 219057a
Show file tree
Hide file tree
Showing 5 changed files with 113 additions and 32 deletions.
4 changes: 2 additions & 2 deletions source/Makefile
View file
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,7 @@ LDFLAGS+= -L$(LIB) -lm -lgsl -lgslcblas $(CUDA_LIBS)
#Note that version should be a single string without spaces.


VERSION = 87e
VERSION = 87f


startup:
Expand Down Expand Up @@ -448,4 +448,4 @@ libatomic.a: atomicdata.o atomicdata_init.o atomic.o

clean :
rm -f *.o *~
cd tests; make clean
cd tests; make clean
11 changes: 7 additions & 4 deletions source/python.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1036,10 +1036,13 @@ typedef struct plasma
and freqmax. This will depend on the last call to total_emission */
double heat_shock; /**< An extra heating term added to allow for shock heating of the plasma (Implementef for FU Ori Project */

/* these variables are for the BF_SIMPLE_EMISSIVITY_APPROACH
they allow one to inspect the net flow of energy into and from the simple ion
ionization pool */
double bf_simple_ionpool_in, bf_simple_ionpool_out;
double bf_simple_ionpool_in, bf_simple_ionpool_out; /**<Varibles to track net flow of energy
into and from ionization pool
in BF_SIMPLE_EMISSIVITY_APPROACH
*/
#define N_PHOT_PROC 500
int n_bf_in[N_PHOT_PROC],n_bf_out[N_PHOT_PROC];/**<Counters to track bf excitations and de-exitations.
*/

double comp_nujnu; /**< The integral of alpha(nu)nuj(nu) used to
compute compton cooling- only needs computing once per cycle
Expand Down
27 changes: 21 additions & 6 deletions source/resonate.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -587,20 +587,21 @@ kappa_bf (xplasma, freq, macro_all)
*
* @details
*
* This routine is now called before various cylces of the Monte Carlo calculation.
* This routine is now called before various cycles of the Monte Carlo calculation.
* (in run.c) It determines which
* bf processes are worth considering during the calculation that follows.
*
* To do this
* it uses the edge position (must be inside, or close to, the spectral region of
* interest) and the edge opacity (must be greater than a threshold value, currently
* To do this, it uses a typical distance a photon can travel within the cell,
* the threshold x-section, and the density of the ion of interest. It retains
* x-sections if tau at the threshold frequency is greater than
* a value set to
* 10^-6.
*
* The need for the routine is just to prevent wasting time on unimportant bf
* transitions. It is called (currently) from resonate.
* transitions.
*
* ### Notes ###
* The dimensionalty of kbf_use is currently set to NTOP_PHOT, which is large enough
* The dimensionalty of kbf_use is set to nphot_total, which is large enough
* to store every transition if necessary.
*
**********************************************************/
Expand Down Expand Up @@ -1080,6 +1081,15 @@ scatter (p, nres, nnscat)

prob_kpkt = 1. - (phot_top[*nres - NLINES - 1].freq[0] / freq_comoving);

if (*nres - NLINES - 1 >= 0)
{
xplasma->n_bf_in[*nres - NLINES - 1] += 1;


// XXXXXXXXXXXXXXXXXX 117 had and inordinate

}

if (prob_kpkt < 0)
{
/* only report an error for a negative prob_kpkt if it's large-ish in magnitude. see #436 discussion */
Expand Down Expand Up @@ -1131,6 +1141,11 @@ scatter (p, nres, nnscat)
}
}


if (*nres - NLINES - 1 >= 0)
{
xplasma->n_bf_out[*nres - NLINES - 1] += 1;
}
}
else
{
Expand Down
64 changes: 45 additions & 19 deletions source/run.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,8 +49,8 @@ calculate_ionization (restart_stat)
int restart_stat;
{
int n, nn;
double zz, z_abs_all, z_abs[N_ISTAT], z_else, ztot;
double radiated[20];
double zz, z_abs_all, z_abs_all_orig, z_orig[N_ISTAT], z_abs[N_ISTAT], z_else, z_else_orig, ztot;
double radiated[20], radiated_orig[20];
int nphot_istat[N_ISTAT];
WindPtr w;
PhotPtr p;
Expand Down Expand Up @@ -227,61 +227,87 @@ calculate_ionization (restart_stat)
/* Determine how much energy was absorbed in the wind. first zero counters.
There are counters for total energy absorbed and for each entry in the istat enum,
The second loop is for the energy radiated (i.e. that actually escapes) */
z_abs_all = z_else = 0.0;
z_abs_all = z_else = z_abs_all_orig = z_else_orig = 0.0;
for (nn = 0; nn < N_ISTAT; nn++)
{
z_abs[nn] = 0.0;
z_orig[nn] = 0.0;
nphot_istat[nn] = 0.0;
}
for (nn = 0; nn < 20; nn++)
{
/* 20 entries to match declaration above */
radiated[nn] = 0.0;
radiated_orig[nn] = 0.0;
}

/* loop over the different photon istats to determine where the luminosity went */
for (nn = 0; nn < NPHOT; nn++)
{

z_abs_all += p[nn].w;
z_abs_all_orig += p[nn].w_orig;

/* we want the istat to be >1 (not P_SCAT or P_INWIND) */
if (p[nn].istat < N_ISTAT && p[nn].istat > 1)
// if (p[nn].istat < N_ISTAT && p[nn].istat > 1)
if (p[nn].istat < N_ISTAT)
{
z_abs[p[nn].istat] += p[nn].w;
z_orig[p[nn].istat] += p[nn].w_orig;
nphot_istat[p[nn].istat]++;
}
if (p[nn].istat == P_ESCAPE)
{
radiated[p[nn].origin] += p[nn].w;
radiated_orig[p[nn].origin] += p[nn].w_orig;
}
else
{
z_else += p[nn].w;
z_else_orig += p[nn].w_orig;
}
}

Log
("!!python: Total photon luminosity after transphot %18.12e (absorbed or lost %18.12e). Radiated luminosity %18.12e\n",
z_abs_all, z_abs_all - zz, z_abs[P_ESCAPE]);
if (geo.rt_mode == RT_MODE_MACRO)

for (nn = 0; nn < N_ISTAT; nn++)
{
Log ("!!python: luminosity lost by adiabatic kpkt destruction %18.12e number of packets %d\n", z_abs[P_ADIABATIC],
nphot_istat[P_ADIABATIC]);
Log ("!!python: luminosity lost to low-frequency free-free %18.12e number of packets %d\n", z_abs[P_LOFREQ_FF],
nphot_istat[P_LOFREQ_FF]);
Log ("XXX stat %8d %8d %12.3e %12.3e\n", nn, nphot_istat[nn], z_abs[nn], z_orig[nn]);
}
for (nn = 0; nn < 20; nn++)
{
Log ("XXX rad %8d %12.3e %12.3e\n", nn, radiated[nn], radiated_orig[nn]);
}
Log ("XXX rad abs_all %12.3e %12.3e\n", z_abs_all, z_abs_all_orig);
Log ("XXX rad else l %12.3e %12.3e\n", z_else, z_else_orig);





Log
("!!python: luminosity (radiated or lost) after transphot %18.12e (absorbed or lost %18.12e %18.12e). \n",
z_abs_all, z_abs_all - zz, z_abs_all - z_abs_all_orig);
Log ("\n");
Log ("!!python: stellar photon luminosity escaping %18.12e \n", radiated[PTYPE_STAR]);
Log ("!!python: boundary layer photon luminosity escaping %18.12e \n", radiated[PTYPE_BL]);
Log ("!!python: disk photon luminosity escaping %18.12e \n", radiated[PTYPE_DISK]);
Log ("!!python: wind photon luminosity escaping %18.12e \n", radiated[PTYPE_WIND]);
Log ("!!python: agn photon luminosity escaping %18.12e \n", radiated[PTYPE_AGN]);
Log ("!!python: luminosity escaping %18.12e\n", z_abs[P_ESCAPE]);
Log ("!!python: stellar photon luminosity escaping %18.12e \n", radiated[PTYPE_STAR] + radiated[PTYPE_STAR_MATOM]);
Log ("!!python: boundary layer photon luminosity escaping %18.12e \n", radiated[PTYPE_BL] + radiated[PTYPE_BL_MATOM]);
Log ("!!python: disk photon luminosity escaping %18.12e \n", radiated[PTYPE_DISK] + radiated[PTYPE_DISK_MATOM]);
Log ("!!python: wind photon luminosity escaping %18.12e \n", radiated[PTYPE_WIND] + radiated[PTYPE_WIND_MATOM]);
Log ("!!python: agn photon luminosity escaping %18.12e \n", radiated[PTYPE_AGN] + radiated[PTYPE_AGN_MATOM]);
Log ("!!python: luminosity lost by any process %18.12e \n", z_else);
Log ("\n");
Log ("!!python: luminosity lost by being completely absorbed %18.12e \n", z_abs[P_ABSORB]);
Log ("!!python: luminosity lost by too many scatters %18.12e \n", z_abs[P_TOO_MANY_SCATTERS]);
Log ("!!python: luminosity lost by hitting the central object %18.12e \n", z_abs[P_HIT_STAR]);
Log ("!!python: luminosity lost by hitting the disk %18.12e \n", z_abs[P_HIT_DISK]);
if (geo.rt_mode == RT_MODE_MACRO)
{
Log ("!!python: luminosity lost by adiabatic kpkt destruction %18.12e number of packets %d\n", z_abs[P_ADIABATIC],
nphot_istat[P_ADIABATIC]);
Log ("!!python: luminosity lost to low-frequency free-free %18.12e number of packets %d\n", z_abs[P_LOFREQ_FF],
nphot_istat[P_LOFREQ_FF]);
}
Log ("!!python: luminosity lost by errors %18.12e \n",
z_abs[P_ERROR] + z_abs[P_ERROR_MATOM] + z_abs[P_REPOSITION_ERROR]);
Log ("!!python: luminosity lost by the unknown %18.12e \n", z_else);
if (geo.binary == TRUE)
Log ("!!python: luminosity lost by hitting the secondary %18.12e \n", z_abs[P_SEC]);

Expand Down
39 changes: 38 additions & 1 deletion source/wind_updates2d.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,6 +110,8 @@ WindPtr (w);
8 * (n_inner_tot + 10 * nions + nlte_levels + 3 * nphot_total + 15 * NXBANDS + 133 * NFLUX_ANGLES) * (floor (NPLASMA / np_mpi_global) +
1);

size_of_commbuffer += 2 * 4 * NPLASMA * N_PHOT_PROC; /* Add space for storing number of bf transitions, up and down */

size_of_specbuffer = 8 * NPLASMA * NBINS_IN_CELL_SPEC;

size_of_commbuffer += size_of_specbuffer;
Expand Down Expand Up @@ -378,6 +380,11 @@ WindPtr (w);
MPI_Pack (&plasmamain[n].xi, 1, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n].bf_simple_ionpool_in, 1, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n].bf_simple_ionpool_out, 1, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);

MPI_Pack (&plasmamain[n].n_bf_in, N_PHOT_PROC, MPI_INT, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (&plasmamain[n].n_bf_out, N_PHOT_PROC, MPI_INT, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);


MPI_Pack (plasmamain[n].cell_spec_flux, NBINS_IN_CELL_SPEC, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (plasmamain[n].F_UV_ang_x, NFLUX_ANGLES, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
MPI_Pack (plasmamain[n].F_UV_ang_y, NFLUX_ANGLES, MPI_DOUBLE, commbuffer, size_of_commbuffer, &position, MPI_COMM_WORLD);
Expand Down Expand Up @@ -527,6 +534,8 @@ WindPtr (w);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].xi, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].bf_simple_ionpool_in, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].bf_simple_ionpool_out, 1, MPI_DOUBLE, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].n_bf_in, N_PHOT_PROC, MPI_INT, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, &plasmamain[n].n_bf_out, N_PHOT_PROC, MPI_INT, MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, plasmamain[n].cell_spec_flux, NBINS_IN_CELL_SPEC, MPI_DOUBLE,
MPI_COMM_WORLD);
MPI_Unpack (commbuffer, size_of_commbuffer, &position, plasmamain[n].F_UV_ang_x, NFLUX_ANGLES, MPI_DOUBLE, MPI_COMM_WORLD);
Expand Down Expand Up @@ -1245,6 +1254,13 @@ wind_rad_init ()
plasmamain[n].bf_simple_ionpool_out = 0.0;
plasmamain[n].bf_simple_ionpool_in = 0.0;

for (i = 0; i < nphot_total; i++)
{
plasmamain[n].n_bf_in[i] = 0;
plasmamain[n].n_bf_out[i] = 0;
}


for (i = 0; i < 3; i++)
plasmamain[n].dmo_dt[i] = 0.0;
for (i = 0; i < 4; i++)
Expand Down Expand Up @@ -1399,7 +1415,8 @@ wind_rad_init ()
int
report_bf_simple_ionpool ()
{
int n;
int n, m;
int in_tot, out_tot;
double total_in = 0.0;
double total_out = 0.0;

Expand All @@ -1417,5 +1434,25 @@ report_bf_simple_ionpool ()

Log ("!! report_bf_simple_ionpool: Total flow into: %8.4e and out of: %8.4e bf_simple ion pool\n", total_in, total_out);

total_in = total_out = 0;
for (m = 0; m < nphot_total; m++)
{
in_tot = out_tot = 0;
for (n = 0; n < NPLASMA; n++)
{
in_tot += plasmamain[n].n_bf_in[m];
out_tot += plasmamain[n].n_bf_out[m];
}


Log ("!! report_bf: %3d %3d %3d %7d %7d\n", m, phot_top[m].z, phot_top[m].istate, in_tot, out_tot);

total_in += in_tot;
total_out += out_tot;
}

Log ("!! report_bf tots: %10.0f %10.0f\n", total_in, total_out);


return (0);
}

0 comments on commit 219057a

Please sign in to comment.