diff --git a/source/Makefile b/source/Makefile index eb94a0212..06b605073 100644 --- a/source/Makefile +++ b/source/Makefile @@ -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: @@ -448,4 +448,4 @@ libatomic.a: atomicdata.o atomicdata_init.o atomic.o clean : rm -f *.o *~ - cd tests; make clean \ No newline at end of file + cd tests; make clean diff --git a/source/python.h b/source/python.h index 85017402b..9473ae5d4 100644 --- a/source/python.h +++ b/source/python.h @@ -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; /**= 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 */ @@ -1131,6 +1141,11 @@ scatter (p, nres, nnscat) } } + + if (*nres - NLINES - 1 >= 0) + { + xplasma->n_bf_out[*nres - NLINES - 1] += 1; + } } else { diff --git a/source/run.c b/source/run.c index 68e436f3b..38176bc8a 100644 --- a/source/run.c +++ b/source/run.c @@ -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; @@ -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]); diff --git a/source/wind_updates2d.c b/source/wind_updates2d.c index 8b9bd0c9b..72c9708c5 100644 --- a/source/wind_updates2d.c +++ b/source/wind_updates2d.c @@ -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; @@ -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); @@ -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); @@ -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++) @@ -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; @@ -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); }