diff --git a/README.md b/README.md index 05f3310..32c1f80 100644 --- a/README.md +++ b/README.md @@ -58,6 +58,10 @@ equations of state and opacities and emissivities. - An HDF5 reader exists for opacities defined in `nubhlight` format. See `core/opac_emis_hdf.c` for the reader and for information for how to define your own such table. +- The HDF5 opacity reader does not exactly support NuLib opacities, + however, [a fork of + NuLib](https://github.com/Yurlungur/NuLib/pull/1) exists that + contains the relevant modifications to run with nublhight. # CLANG FORMAT diff --git a/core/constants.h b/core/constants.h index bed0f3e..39f32cf 100644 --- a/core/constants.h +++ b/core/constants.h @@ -28,6 +28,10 @@ #define S4THW (S2THW * S2THW) #define NUSIGMA0 (1.7611737037e-44) // Fundamental neutrino cross section +// Frequency scale of neutrino oscillations +#define ROOT2 (1.4142135623730951) +#define NUFERM (ROOT2*HBAR*HBAR*CL*CL*CL*GFERM) + // Unit conversions #define EV (1.60217653e-12) // Electron-volt #define MEV (1.0e6 * EV) // Mega-Electron-Volt diff --git a/core/coord.c b/core/coord.c index aea5f19..7031388 100644 --- a/core/coord.c +++ b/core/coord.c @@ -402,7 +402,12 @@ void set_points() { stopx_rad[1] = startx_rad[1] + N1TOT * dx[1]; stopx_rad[2] = startx_rad[2] + N2TOT * dx[2]; stopx_rad[3] = startx_rad[3] + N3TOT * dx[3]; -#endif +#if LOCAL_ANGULAR_DISTRIBUTIONS + local_dx1_rad = (stopx_rad[1] - startx_rad[1]) / (LOCAL_ANGLES_NX1); + local_dx2_rad = (stopx_rad[2] - startx_rad[2]) / (LOCAL_ANGLES_NX2); + local_dx_costh = 2. / (LOCAL_ANGLES_NMU); +#endif // LOCAL_ANGULAR_DISTRIBUTIONS +#endif // RADIATION #elif METRIC == MKS // Calculate some radii determined by the geometry Reh = 1. + sqrt(1. - a * a); @@ -435,6 +440,11 @@ void set_points() { stopx_rad[1] = log(Rout_rad); stopx_rad[2] = startx[2] + N2TOT * dx[2]; stopx_rad[3] = startx[3] + N3TOT * dx[3]; +#if LOCAL_ANGULAR_DISTRIBUTIONS + local_dx1_rad = (stopx_rad[1] - startx_rad[1]) / (LOCAL_ANGLES_NX1); + local_dx2_rad = (stopx_rad[2] - startx_rad[2]) / (LOCAL_ANGLES_NX2); + local_dx_costh = 2. / (LOCAL_ANGLES_NMU); +#endif // LOCAL_ANGULAR_DISTRIBUTIONS #endif poly_norm = 0.5 * M_PI * 1. / (1. + 1. / (poly_alpha + 1.) * 1. / pow(poly_xt, poly_alpha)); diff --git a/core/decs.h b/core/decs.h index e55e35d..5555f66 100644 --- a/core/decs.h +++ b/core/decs.h @@ -165,10 +165,10 @@ #define RAD_TYPE_START (0) #define TYPE_TRACER (-1) #if RADIATION == RADTYPE_NEUTRINOS -#define RAD_NUM_TYPES (3) #define NU_ELECTRON (0) #define ANTINU_ELECTRON (1) #define NU_HEAVY (2) +#define ANTINU_HEAVY (3) #if MULTISCATT_TEST #define RAD_SCATT_TYPES (3) #else @@ -183,12 +183,10 @@ #define RADG_YE (4) #define RADG_YE_EM (5) #elif RADIATION == RADTYPE_LIGHT -#define RAD_NUM_TYPES (1) #define RAD_SCATT_TYPES (1) #define NRADCOMP (0) #define PHOTON (0) #else -#define RAD_NUM_TYPES (0) #define RAD_SCATT_TYPES (0) #define NRADCOMP (0) #endif @@ -368,6 +366,30 @@ extern grid_radtype_type Nem_phys, Nabs_phys, radtype_buf; extern grid_int_type Nsuper; extern grid_double_type Esuper; extern grid_prim_type psupersave; + +#if LOCAL_ANGULAR_DISTRIBUTIONS +#define LOCAL_NUM_BASES (2) +#define MOMENTS_A (0) +#define MOMENTS_B (1) +#define MOMENTS_DIFF (2) +typedef double grid_local_angles_type[LOCAL_NUM_BASES][LOCAL_ANGLES_NX1] + [LOCAL_ANGLES_NX2][RAD_NUM_TYPES] + [LOCAL_ANGLES_NMU]; +extern grid_local_angles_type local_angles; +extern double local_dx1_rad, local_dx2_rad, local_dx_costh; + +#if RAD_NUM_TYPES >= 4 +typedef double grid_Gnu_type[LOCAL_NUM_BASES][LOCAL_ANGLES_NX1] + [LOCAL_ANGLES_NX2][LOCAL_ANGLES_NMU]; +typedef double grid_local_moment_type[LOCAL_NUM_BASES][3][LOCAL_ANGLES_NX1] + [LOCAL_ANGLES_NX2]; +typedef int grid_local_basis_idx_type[LOCAL_ANGLES_NX1][LOCAL_ANGLES_NX2]; +extern grid_Gnu_type Gnu, local_stddev; +extern grid_local_moment_type local_moments; +extern grid_local_basis_idx_type local_b_osc; +#endif // #if RAD_NUM_TYPES >= 4 +#endif // LOCAL_ANGULAR_DISTRIBUTIONS + #endif // RADIATION // Default initialization is 0, which in this case is @@ -671,6 +693,12 @@ extern int global_stop[NDIM]; #define JRADLOOP for (int n = 0; n < MAXNSCATT + 2; n++) #define NULOOP for (int inu = 0; inu < NU_BINS + 1; inu++) +#define LOCALXLOOP \ + for (int i = 0; i < LOCAL_ANGLES_NX1; ++i) \ + for (int j = 0; j < LOCAL_ANGLES_NX2; ++j) +#define LOCALMULOOP for (int imu = 0; imu < LOCAL_ANGLES_NMU; ++imu) +#define LOCALXMULOOP LOCALXLOOP LOCALMULOOP + #define MY_MIN(fval1, fval2) (((fval1) < (fval2)) ? (fval1) : (fval2)) #define MY_MAX(fval1, fval2) (((fval1) > (fval2)) ? (fval1) : (fval2)) #define MY_SIGN(fval) (((fval) < 0.) ? -1. : 1.) @@ -981,6 +1009,20 @@ double Jnu_hdf(double nu, int type, const struct of_microphysics *m); double int_jnudnudOmega_hdf(const struct of_microphysics *m); double alpha_nu_hdf(double nu, int type, const struct of_microphysics *m); #endif // HDF opacities + +// oscillations.c +#if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS +double get_dt_oscillations(); +void get_local_angle_bins( + struct of_photon *ph, int *pi, int *pj, int *pmu1, int *pmu2); +void accumulate_local_angles(); +#if RAD_NUM_TYPES >= 4 +void compute_local_gnu(grid_local_angles_type local_angles, + grid_Gnu_type local_stddev, grid_Gnu_type gnu); +void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments); +void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu); +#endif // RAD_NUM_TYPES >= 4 +#endif // LOCAL_ANGULAR_DISTRIBUTIONS #endif // RADIATION // passive.c @@ -1065,6 +1107,7 @@ void set_cooling_time( void record_lepton_flux(const struct of_photon *ph); void check_nu_type(const char *location); int get_lepton_sign(const struct of_photon *ph); +int nu_is_heavy(const int radtype); #endif // NEUTRINOS #endif // RADIATION diff --git a/core/defs.h b/core/defs.h index a26d977..63c0df2 100644 --- a/core/defs.h +++ b/core/defs.h @@ -49,6 +49,16 @@ grid_radtype_type Nem_phys, Nabs_phys, radtype_buf; grid_int_type Nsuper; grid_double_type Esuper; grid_prim_type psupersave; + +#if LOCAL_ANGULAR_DISTRIBUTIONS +grid_local_angles_type local_angles; +double local_dx1_rad, local_dx2_rad, local_dx_costh; +#if RAD_NUM_TYPES >= 4 +grid_Gnu_type Gnu, local_stddev; +grid_local_moment_type local_moments; +grid_local_basis_idx_type local_b_osc; +#endif // RAD_NUM_TYPES +#endif // LOCAL_ANGULAR_DISTRIBUTIONS #endif // RADIATION #if ELECTRONS diff --git a/core/diag.c b/core/diag.c index 370a744..e326bdd 100644 --- a/core/diag.c +++ b/core/diag.c @@ -583,8 +583,12 @@ void print_rad_types() { rad_type_counts[NU_ELECTRON]); fprintf(stdout, " ANTI %.2f %%\n", rad_type_counts[ANTINU_ELECTRON]); - fprintf(stdout, " HEAVY %.2f %%\n", + fprintf(stdout, " X %.2f %%\n", rad_type_counts[NU_HEAVY]); +#if RAD_NUM_TYPES > 3 + fprintf(stdout, " ANTIX %.2f %%\n", + rad_type_counts[ANTINU_HEAVY]); +#endif fprintf(stdout, "*********************************\n\n"); } timer_stop(TIMER_DIAG); diff --git a/core/io.c b/core/io.c index 07ca050..d111939 100644 --- a/core/io.c +++ b/core/io.c @@ -359,7 +359,7 @@ void track_ph() { hsize_t dims[1] = {num_tracked_buf}; hid_t space = H5Screate_simple(1, dims, NULL); ph_dsets[n] = H5Dcreate(file_id, dsetnam, trackphfiletype, space, - H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); H5Sclose(space); } if (num_tracked > 0) @@ -637,6 +637,92 @@ void dump_grid() { free(Lambda_h2cart_cov); } +#if RADIATION +#if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS + { + double *local_angles_Xharm = safe_malloc( + (NDIM - 1) * LOCAL_ANGLES_NX1 * LOCAL_ANGLES_NX2 * sizeof(double)); +#if METRIC == MKS + double *local_angles_Xbl = safe_malloc( + (NDIM - 1) * LOCAL_ANGLES_NX1 * LOCAL_ANGLES_NX2 * sizeof(double)); +#endif // MKS + double *local_angles_Xcart = safe_malloc( + (NDIM - 1) * LOCAL_ANGLES_NX1 * LOCAL_ANGLES_NX2 * sizeof(double)); + int n = 0; + double X[NDIM], Xcart[NDIM], r, th; + for (int i = 0; i < LOCAL_ANGLES_NX1; ++i) { + X[1] = startx_rad[1] + (i + 0.5) * local_dx1_rad; // startx is a face + for (int j = 0; j < LOCAL_ANGLES_NX2; ++j) { + X[2] = startx_rad[2] + (j + 0.5) * local_dx2_rad; + cart_coord(X, Xcart); + + local_angles_Xharm[n + 0] = 0; + local_angles_Xharm[n + 1] = X[1]; + local_angles_Xharm[n + 2] = X[2]; + +#if METRIC == MKS + bl_coord(X, &r, &th); + local_angles_Xbl[n + 0] = 0; + local_angles_Xbl[n + 1] = r; + local_angles_Xbl[n + 2] = th; +#endif // MKS + + local_angles_Xcart[n + 0] = 0; + local_angles_Xcart[n + 1] = Xcart[1]; + local_angles_Xcart[n + 2] = Xcart[2]; + + n += (NDIM - 1); + } + } +#define RANK (3) + hsize_t fdims[RANK] = {LOCAL_ANGLES_NX1, LOCAL_ANGLES_NX2, NDIM - 1}; + hsize_t fstart[RANK] = {0, 0, 0}; + hsize_t fcount[RANK] = {LOCAL_ANGLES_NX1, LOCAL_ANGLES_NX2, NDIM - 1}; + hsize_t mdims[RANK] = {LOCAL_ANGLES_NX1, LOCAL_ANGLES_NX2, NDIM - 1}; + hsize_t mstart[RANK] = {0, 0, 0}; + if (!mpi_io_proc()) { + fcount[0] = 0; + fcount[1] = 0; + fcount[2] = 0; + } + WRITE_ARRAY(local_angles_Xharm, RANK, fdims, fstart, fcount, mdims, mstart, + TYPE_DBL); + free(local_angles_Xharm); + + WRITE_ARRAY(local_angles_Xcart, RANK, fdims, fstart, fcount, mdims, mstart, + TYPE_DBL); + free(local_angles_Xcart); + +#if METRIC == MKS + WRITE_ARRAY( + local_angles_Xbl, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); + free(local_angles_Xbl); +#endif // MKS +#undef RANK + } + + { + double *local_angles_mu = safe_malloc(LOCAL_ANGLES_NMU * sizeof(double)); + for (int i = 0; i < LOCAL_ANGLES_NMU; ++i) { + local_angles_mu[i] = -1 + (i + 0.5) * local_dx_costh; + } +#define RANK (1) + hsize_t fdims[RANK] = {LOCAL_ANGLES_NMU}; + hsize_t fstart[RANK] = {0}; + hsize_t fcount[RANK] = {LOCAL_ANGLES_NMU}; + hsize_t mdims[RANK] = {LOCAL_ANGLES_NMU}; + hsize_t mstart[RANK] = {0}; + if (!mpi_io_proc()) { + fcount[0] = 0; + } + WRITE_ARRAY( + local_angles_mu, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); +#undef RANK + free(local_angles_mu); + } +#endif // RADIATION +#endif // LOCAL_ANGULAR_DISTRIBUTIONS + H5Fflush(file_id, H5F_SCOPE_GLOBAL); H5Fclose(file_id); } @@ -816,7 +902,7 @@ void dump() { hsize_t str_dims[1] = {NVAR}; hid_t prim_space = H5Screate_simple(1, str_dims, NULL); hid_t str_attr = H5Acreate( - prim_dset, "vnams", strtype, prim_space, H5P_DEFAULT, H5P_DEFAULT); + prim_dset, "vnams", strtype, prim_space, H5P_DEFAULT, H5P_DEFAULT); H5Awrite(str_attr, strtype, vnams); H5Aclose(str_attr); H5Sclose(prim_space); @@ -975,6 +1061,76 @@ void dump() { WRITE_ARRAY(nuLnu, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); #undef RANK } + + // local angle histograms +#if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS + accumulate_local_angles(); + { +#define RANK (5) + hsize_t fdims[RANK] = {LOCAL_NUM_BASES, LOCAL_ANGLES_NX1, + LOCAL_ANGLES_NX2, RAD_NUM_TYPES, LOCAL_ANGLES_NMU}; + hsize_t fstart[RANK] = {0, 0, 0, 0, 0}; + hsize_t fcount[RANK] = {LOCAL_NUM_BASES, LOCAL_ANGLES_NX1, + LOCAL_ANGLES_NX2, RAD_NUM_TYPES, LOCAL_ANGLES_NMU}; + hsize_t mdims[RANK] = {LOCAL_NUM_BASES, LOCAL_ANGLES_NX1, + LOCAL_ANGLES_NX2, RAD_NUM_TYPES, LOCAL_ANGLES_NMU}; + hsize_t mstart[RANK] = {0, 0, 0, 0, 0}; + if (!mpi_io_proc()) { + fcount[0] = 0; + fcount[1] = 0; + fcount[2] = 0; + fcount[3] = 0; + fcount[4] = 0; + } + WRITE_ARRAY( + local_angles, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); +#undef RANK + } + + { +#define RANK (4) + hsize_t fdims[RANK] = {LOCAL_NUM_BASES, LOCAL_ANGLES_NX1, + LOCAL_ANGLES_NX2, LOCAL_ANGLES_NMU}; + hsize_t fstart[RANK] = {0, 0, 0, 0}; + hsize_t fcount[RANK] = {LOCAL_NUM_BASES, LOCAL_ANGLES_NX1, + LOCAL_ANGLES_NX2, LOCAL_ANGLES_NMU}; + hsize_t mdims[RANK] = {LOCAL_NUM_BASES, LOCAL_ANGLES_NX1, + LOCAL_ANGLES_NX2, LOCAL_ANGLES_NMU}; + hsize_t mstart[RANK] = {0, 0, 0, 0}; + if (!mpi_io_proc()) { + fcount[0] = 0; + fcount[1] = 0; + fcount[2] = 0; + fcount[3] = 0; + } + WRITE_ARRAY(Gnu, + RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); + WRITE_ARRAY(local_stddev, + RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); +#undef RANK + } + + { +#define RANK (4) + hsize_t fdims[RANK] = { + LOCAL_NUM_BASES, 2, LOCAL_ANGLES_NX1, LOCAL_ANGLES_NX2}; + hsize_t fstart[RANK] = {0, 0, 0, 0}; + hsize_t fcount[RANK] = { + LOCAL_NUM_BASES, 2, LOCAL_ANGLES_NX1, LOCAL_ANGLES_NX2}; + hsize_t mdims[RANK] = { + LOCAL_NUM_BASES, 2, LOCAL_ANGLES_NX1, LOCAL_ANGLES_NX2}; + hsize_t mstart[RANK] = {0, 0, 0, 0}; + if (!mpi_io_proc()) { + fcount[0] = 0; + fcount[1] = 0; + fcount[2] = 0; + fcount[3] = 0; + } + WRITE_ARRAY( + local_moments, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); +#undef RANK + } +#endif // LOCAL_ANGULAR_DISTRIBUTIONS #endif // RADIATION } @@ -1271,8 +1427,8 @@ void restart_write(int restart_type) { hsize_t fcount[RANK] = {1}; hsize_t mdims[RANK] = {1}; hsize_t mstart[RANK] = {0}; - int *particle_offsets = safe_malloc(1 * sizeof(int)); - int *particle_counts = safe_malloc(1 * sizeof(int)); + int * particle_offsets = safe_malloc(1 * sizeof(int)); + int * particle_counts = safe_malloc(1 * sizeof(int)); particle_offsets[0] = npart_offset; particle_counts[0] = npart_local; WRITE_ARRAY( @@ -1525,8 +1681,8 @@ void restart_read(char *fname) { hsize_t fcount[RANK] = {1}; hsize_t mdims[RANK] = {1}; hsize_t mstart[RANK] = {0}; - int *particle_offsets = safe_malloc(1 * sizeof(int)); - int *particle_counts = safe_malloc(1 * sizeof(int)); + int * particle_offsets = safe_malloc(1 * sizeof(int)); + int * particle_counts = safe_malloc(1 * sizeof(int)); READ_ARRAY( particle_offsets, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_INT); READ_ARRAY( @@ -1774,9 +1930,9 @@ void dump_tracers() { WRITE_HDR(ntracers, TYPE_INT); // arrays to fil - int *id = safe_malloc(ntracers_local * sizeof(int)); - int *it = safe_malloc(ntracers_local * sizeof(int)); - int *active = safe_malloc(ntracers_local * sizeof(int)); + int * id = safe_malloc(ntracers_local * sizeof(int)); + int * it = safe_malloc(ntracers_local * sizeof(int)); + int * active = safe_malloc(ntracers_local * sizeof(int)); double *time = safe_malloc(ntracers_local * sizeof(double)); double *mass = safe_malloc(ntracers_local * sizeof(double)); double *Xharm = safe_malloc(3 * ntracers_local * sizeof(double)); @@ -2072,7 +2228,7 @@ void write_array(void *data, const char *name, hsize_t rank, hsize_t *fdims, H5Tset_size(string_type, strlen(data)); plist_id = H5Pcreate(H5P_DATASET_CREATE); dset_id = H5Dcreate(file_id, name, string_type, filespace, H5P_DEFAULT, - plist_id, H5P_DEFAULT); + plist_id, H5P_DEFAULT); H5Pclose(plist_id); plist_id = H5Pcreate(H5P_DATASET_XFER); H5Dwrite(dset_id, string_type, memspace, filespace, plist_id, data); diff --git a/core/oscillations.c b/core/oscillations.c new file mode 100644 index 0000000..10a0f7c --- /dev/null +++ b/core/oscillations.c @@ -0,0 +1,229 @@ +/****************************************************************************** + * * + * OSCILLATIONS.C * + * * + * Routines for neutrino oscillations * + * * + ******************************************************************************/ +#include "decs.h" + +#if RADIATION == RADTYPE_NEUTRINOS +#if LOCAL_ANGULAR_DISTRIBUTIONS + +double get_dt_oscillations() { + set_Rmunu(); // So we have Nsph and nph + double nph_max = 0; +#pragma omp parallel for reduction(max : nph_max) collapse(3) + ZLOOP { + nph_max = MY_MAX(nph_max, nph[i][j][k]); // 1/cm^3 + } + // seconds + double dt_osc = 1. / (NUFERM * nph_max + SMALL); + dt_osc /= T_unit; // code units + dt_osc = mpi_min(dt_osc); + return dt_osc; +} + +void accumulate_local_angles() { + static const int LOCAL_ANGLES_SIZE = LOCAL_NUM_BASES * LOCAL_ANGLES_NX1 * + LOCAL_ANGLES_NX2 * LOCAL_ANGLES_NMU * + RAD_NUM_TYPES; + static const int LOCAL_STDDEV_SIZE = + LOCAL_NUM_BASES * LOCAL_ANGLES_NX1 * LOCAL_ANGLES_NX2 * LOCAL_ANGLES_NMU; + + memset(local_angles, 0, LOCAL_ANGLES_SIZE * sizeof(double)); + memset(local_stddev, 0, LOCAL_STDDEV_SIZE * sizeof(double)); + +#pragma omp parallel + { + struct of_photon *ph = photon_lists[omp_get_thread_num()]; + while (ph != NULL) { + if (ph->type != TYPE_TRACER) { + int ix1, ix2, icosth[LOCAL_NUM_BASES]; + get_local_angle_bins(ph, &ix1, &ix2, &icosth[0], &icosth[1]); + for (int b = 0; b < LOCAL_NUM_BASES; ++b) { +#pragma omp atomic + local_angles[b][ix1][ix2][ph->type][icosth[b]] += ph->w; +#pragma omp atomic + local_stddev[b][ix1][ix2][icosth[b]] += 1.; + } + } + ph = ph->next; + } + } // omp parallel + + mpi_dbl_allreduce_array((double *)local_angles, LOCAL_ANGLES_SIZE); + mpi_dbl_allreduce_array((double *)local_stddev, LOCAL_STDDEV_SIZE); + +#pragma omp parallel for collapse(4) + for (int b = 0; b < LOCAL_NUM_BASES; ++b) { + LOCALXMULOOP { + local_stddev[b][i][j][imu] = sqrt(fabs(local_stddev[b][i][j][imu])); + } + } + + // Gnu, local_moments are global +#if RAD_NUM_TYPES >= 4 + compute_local_gnu(local_angles, local_stddev, Gnu); + compute_local_moments(Gnu, local_moments); +#endif // RAD_NUM_TYPES >= 4 +} + +void get_local_angle_bins( + struct of_photon *ph, int *pi, int *pj, int *pmu1, int *pmu2) { + double X[NDIM]; + double Kcov[NDIM]; + double Kcon[NDIM]; + int k; + get_X_K_interp(ph, t, P, X, Kcov, Kcon); + Xtoijk(X, pi, pj, &k); + + // JMM: If we use more complicated bases this is more complicated + // change this for other basis vectors + double X1norm = sqrt(fabs(ggeom[*pi][*pj][CENT].gcov[1][1])); + double X2norm = sqrt(fabs(ggeom[*pi][*pj][CENT].gcov[2][2])); + double X1vec[NDIM] = {0, 1. / (X1norm + SMALL), 0, 0}; + double X2vec[NDIM] = {0, 0, 1. / (X2norm + SMALL), 0}; + + double costh1 = 0; + double costh2 = 0; + SDLOOP { + costh1 += X1vec[mu] * Kcov[mu]; // X1^a K_a + costh2 += X2vec[mu] * Kcov[mu]; // X2^a K_a + } + // cos(th) = e^a K_a / |e| |K| + // |e| = 1 by construction, but K must be normalized + // note we want to normalize the SPATIAL part of K + double knorm = 0; + for (int mu = 1; mu < NDIM; ++mu) { + for (int nu = 1; nu < NDIM; ++nu) { + knorm += fabs((ggeom[*pi][*pj][CENT].gcov[mu][nu]) * Kcon[mu] * Kcon[mu]); + } + } + // sqrt the inner product and lets go + knorm = sqrt(knorm); + knorm = 1. / (fabs(knorm) + SMALL); + costh1 *= knorm; + costh2 *= knorm; + + *pi = MY_MAX( + 0, MY_MIN(LOCAL_ANGLES_NX1 - 1, (X[1] - startx_rad[1]) / local_dx1_rad)); + *pj = MY_MAX( + 0, MY_MIN(LOCAL_ANGLES_NX2 - 1, (X[2] - startx_rad[2]) / local_dx2_rad)); + *pmu1 = + MY_MAX(0, MY_MIN(LOCAL_ANGLES_NMU - 1, (costh1 + 1) / local_dx_costh)); + *pmu2 = + MY_MAX(0, MY_MIN(LOCAL_ANGLES_NMU - 1, (costh2 + 1) / local_dx_costh)); +} + +#if RAD_NUM_TYPES >= 4 +void compute_local_gnu( + grid_local_angles_type f, grid_Gnu_type stddev, grid_Gnu_type gnu) { +#pragma omp parallel for collapse(4) + for (int b = 0; b < LOCAL_NUM_BASES; ++b) { + LOCALXMULOOP { + // TODO(JMM): Generalize this for six species? + double ELN = + (f[b][i][j][NU_ELECTRON][imu] - f[b][i][j][ANTINU_ELECTRON][imu]); + double XLN = (f[b][i][j][NU_HEAVY][imu] - f[b][i][j][ANTINU_HEAVY][imu]); + + double tot = 0; + TYPELOOP tot += f[b][i][j][itp][imu]; + double ebar = tot / (stddev[b][i][j][imu] + SMALL); + + double g_temp = ELN - 0.5*XLN; + gnu[b][i][j][imu] = (fabs(g_temp) > ebar) * g_temp; + } + } +} + +// JMM: We can also compute, e.g., the average bin momentum if we need +// to, e.g., compute higher moment integrands +void compute_local_moments(grid_Gnu_type gnu, grid_local_moment_type moments) { + // We are reducing over mu, but if we just parallel loop over b,i,j, + // there is no danger of index collisions. + LOCALMULOOP { +#pragma omp parallel for collapse(3) + for (int b = 0; b < LOCAL_NUM_BASES; ++b) { + LOCALXLOOP { + if (gnu[b][i][j][imu] < 0) { + // TODO(JMM): Pretty sure this atomic isn't needed +#pragma omp atomic + moments[b][MOMENTS_A][i][j] += fabs(gnu[b][i][j][imu]); + } else if (gnu[b][i][j][imu] > 0) { + // TODO(JMM): Pretty sure this atomic isn't needed +#pragma omp atomic + moments[b][MOMENTS_B][i][j] += fabs(gnu[b][i][j][imu]); + } // else nada. We don't care about == 0. + } + } + } // MU LOOP + +#pragma omp parallel for collapse(2) + LOCALXLOOP { + for (int b = 0; b < LOCAL_NUM_BASES; ++b) { + moments[b][MOMENTS_DIFF][i][j] = + fabs(moments[b][MOMENTS_B][i][j] - moments[b][MOMENTS_A][i][j]); + } + local_b_osc[i][j] = (local_moments[1][MOMENTS_DIFF][i][j] > + local_moments[0][MOMENTS_DIFF][i][j]); + } +} + +void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) { +#pragma omp parallel + { + struct of_photon *ph = photon_lists[omp_get_thread_num()]; + while (ph != NULL) { + if (ph->type != TYPE_TRACER) { + int ix1, ix2, icosth[LOCAL_NUM_BASES]; + get_local_angle_bins(ph, &ix1, &ix2, &icosth[0], &icosth[1]); + + int b_osc = local_b_osc[ix1][ix2]; + int imu = icosth[b_osc]; + double A = local_moments[b_osc][MOMENTS_A][ix1][ix2]; + double B = local_moments[b_osc][MOMENTS_B][ix1][ix2]; + double G = gnu[b_osc][ix1][ix2][imu]; + + // gnu == 0 when we activated stddev trigger. Don't oscillate. + if (((G != 0) || FORCE_EQUIPARTITION) && (A != 0) && (B != 0)) { + // if ((A != 0) && (B != 0)) { + // If A == B then which region we treat as shallow is + // unimportant. Psurvive = 1/3 for both regions. + int A_is_shallow = A < B; + int B_is_shallow = !(A_is_shallow); + double shallow = A_is_shallow ? A : B; + double deep = A_is_shallow ? B : A; + + int g_in_A = G < 0; + int g_in_B = G > 0; + + int in_shallow = (A_is_shallow && g_in_A) || (B_is_shallow && g_in_B); + int in_deep = !(in_shallow); + + double peq = nu_is_heavy(ph->type) ? (2./3.) : (1./3.); +#if FORCE_EQUIPARTITION + double p_survival = peq; +#else + double p_survival = + in_shallow ? peq : (1 - (1 - peq) * B / (A + SMALL)); +#endif // FORCE_EQUIPARTITION + double p_osc = 1. - p_survival; + if (get_rand() < p_osc) { + // JMM: + // Type order is NUE, NUEBAR, NUX, NUXBAR + // adding 2 on ring 0, 1, 2, 3 + // moves through without changing to antiparticle. + ph->type = (ph->type + (RAD_NUM_TYPES / 2)) % RAD_NUM_TYPES; + } + } + } + ph = ph->next; + } + } // omp parallel +} + +#endif // RAD_NUM_TYPES >= 4 + +#endif // LOCAL_ANGULAR_DISTRIBUTIONS +#endif // RADIATION diff --git a/core/rad_utils.c b/core/rad_utils.c index 19bda36..5c1d877 100644 --- a/core/rad_utils.c +++ b/core/rad_utils.c @@ -758,6 +758,9 @@ int get_lepton_sign(const struct of_photon *ph) { return -1; return 0; } +int nu_is_heavy(const int radtype) { + return ((radtype == NU_HEAVY) || (radtype == ANTINU_HEAVY)); +} // for debugging void check_nu_type(const char *location) { diff --git a/core/scattering.c b/core/scattering.c index 84504c0..2e24772 100644 --- a/core/scattering.c +++ b/core/scattering.c @@ -841,7 +841,7 @@ double total_cross_lkup( } #else // Normal neutrino scattering { - if (type == NU_HEAVY && interaction == RSCATT_TYPE_E) { + if (nu_is_heavy(type) && interaction == RSCATT_TYPE_E) { return 0.0; // heavy cannot scatter off of electrons } double sigma = @@ -1032,7 +1032,7 @@ double total_cross_ions(double sigma_hc, double A, double Z) { // Burrows, Reddy, Thomson, arXiv:astro-ph/0404432 double nu_cross_delta(int type, int interaction) { // heavy cannot scatter off of electrons - if (type == NU_HEAVY && interaction == RSCATT_TYPE_E) + if (nu_is_heavy(type) && interaction == RSCATT_TYPE_E) return 0.0; if (interaction == RSCATT_TYPE_P) { double Cpv = 0.5 + 2.0 * S2THW; @@ -1054,7 +1054,7 @@ double nu_cross_delta(int type, int interaction) { double nu_cross_factor( double sigma, int type, int interaction, const struct of_microphysics *m) { // heavy cannot scatter off of electrons - if (type == NU_HEAVY && interaction == RSCATT_TYPE_E) + if (nu_is_heavy(type) && interaction == RSCATT_TYPE_E) return 0.0; if (interaction == RSCATT_TYPE_P) { return 0.25 * NUSIGMA0 * sigma * diff --git a/core/step.c b/core/step.c index dd3a059..983c16a 100644 --- a/core/step.c +++ b/core/step.c @@ -27,7 +27,13 @@ void step() { double ndt; #if RADIATION double dt_cool; -#endif +#if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS && \ + RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS + // not used for timestep control. Used to turn oscillations on or + // off. + double dt_osc = get_dt_oscillations(); +#endif // oscillations +#endif // radiation dtsave = dt; // Need both P_n and P_n+1 to calculate current @@ -65,7 +71,20 @@ void step() { interact(Ph, extra, t, dt); // check_nu_type("after interact"); // DEBUG bound_superphotons(Ph, t, dt); -// check_nu_type("after bound"); // DEBUG + // check_nu_type("after bound"); // DEBUG +#if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS && \ + RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS + int oscillations_active = (dt_osc <= dt); + if (mpi_io_proc()) { + printf("\t[Oscillations] Active? %d dt_osc = %.14e\n", oscillations_active, + dt_osc); + } + if (oscillations_active) { // TOOD(JMM): Some safety factor? + accumulate_local_angles(); + oscillate(local_moments, Gnu); + } + // check_nu_type("after oscillate"); // DEBUG +#endif // OSCILLATIONS #endif // Corrector step diff --git a/prob/oscillations/build.py b/prob/oscillations/build.py new file mode 100644 index 0000000..085beb5 --- /dev/null +++ b/prob/oscillations/build.py @@ -0,0 +1,166 @@ +################################################################################ +# # +# UNIT TEST FOR MULTISCATT # +# # +################################################################################ + +import sys +sys.dont_write_bytecode = True +sys.path.append('../../script/') +sys.path.append('../../script/analysis') +import bhlight as bhl +import make_tabulated_gamma as tab +from units import cgs +PROB = 'oscillations' + +MPI = '-mpi' in sys.argv +if '-idim' in sys.argv: + IDIM = int(sys.argv[sys.argv.index('-idim') +1]) +else: + IDIM = 1 + +NEUTRINOS = True +print("MPI = {}".format(MPI)) +print("IDIM = {}".format(IDIM)) + +NCPU = 2 if MPI else 1 +NTOT = 12 if MPI else 3 +TF = 0.1 +DTout = TF +T0 = TF/1e6 +Nsph_tot = 1e7 +Nph_tot = 1e55 +E_MEV = 25 + +FRACS = [0.2, 0.4, 0.2, 0.2] + +# Use a fake table for this test +GAMMA = 1.4 +TFINAL = 600 +YE = 0.5 + +CL = cgs['CL'] +RHOMIN, RHOMAX, NRHO = 1e-4, 1e20, 20 +UMIN, UMAX, NU = 1e-8, 1e8, 20 +YEMIN, YEMAX, NYE = 0.0, 0.6, 10 +CRASH_ON_SOUND_SPEED = False + +rhol = 1.0 +RHO_UNIT = 2.8e14 +L_UNIT = rhol/RHO_UNIT +M_UNIT = RHO_UNIT*(L_UNIT**3) + +print("rhol = {:.5}".format(rhol)) +print("RHO_UNIT = {:.5}".format(RHO_UNIT)) +print("L_UNIT = {:.5}".format(L_UNIT)) +print("M_UNIT = {:.5}".format(M_UNIT)) + +tablepath = tab.make_filename(GAMMA) +units = tab.UnitSystem(M_UNIT, L_unit = L_UNIT) +tab.make_table_u(RHOMIN, RHOMAX, NRHO, + UMIN, UMAX, NU, + YEMIN, YEMAX, NYE, + units, GAMMA, tablepath, + CRASH_ON_SOUND_SPEED) + + +RHO0 = 1. +UU0 = 1e-6*RHO0 +print("rho_cgs = {}\nu_cgs = {}\n".format(RHO0,UU0)) + + ### COMPILE TIME PARAMETERS ### + +# SPATIAL RESOLUTION AND MPI DECOMPOSITION +bhl.config.set_cparm('N1TOT', NTOT) +bhl.config.set_cparm('N2TOT', NTOT) +bhl.config.set_cparm('N3TOT', NTOT) +bhl.config.set_cparm('N1CPU', NCPU) +bhl.config.set_cparm('N2CPU', NCPU) +bhl.config.set_cparm('N3CPU', NCPU) + +# OPENMP PARALLELIZATION +bhl.config.set_cparm('OPENMP', True) + +# COORDINATES +bhl.config.set_cparm('METRIC', 'MINKOWSKI') + +# FLUID +bhl.config.set_cparm('RECONSTRUCTION', 'WENO') +bhl.config.set_cparm('X1L_GAS_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X1R_GAS_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X2L_GAS_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X2R_GAS_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X3L_GAS_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X3R_GAS_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X1L_INFLOW', False) +bhl.config.set_cparm('X1R_INFLOW', False) +bhl.config.set_cparm('X2L_INFLOW', False) +bhl.config.set_cparm('X2R_INFLOW', False) +bhl.config.set_cparm('X3L_INFLOW', False) +bhl.config.set_cparm('X3R_INFLOW', False) + +# EOS +bhl.config.set_cparm("EOS", "EOS_TYPE_TABLE") +bhl.config.set_cparm('NVAR_PASSIVE', 2) + +# RADIATION +bhl.config.set_cparm('RADIATION', 'RADTYPE_NEUTRINOS') +bhl.config.set_cparm('ESTIMATE_THETAE', False) +bhl.config.set_cparm('EMISSION', False) +bhl.config.set_cparm('ABSORPTION', False) +bhl.config.set_cparm('SCATTERING', False) +bhl.config.set_cparm('MULTISCATT_TEST', False) +bhl.config.set_cparm('NU_BINS', 5) +bhl.config.set_cparm('GRAYABSORPTION', False) +bhl.config.set_cparm('BREMSSTRAHLUNG', False) +bhl.config.set_cparm('SYNCHROTRON', False) +bhl.config.set_cparm('NU_BINS_SPEC', 5) +bhl.config.set_cparm("NTH", 8) +bhl.config.set_cparm("NPHI", 8) +bhl.config.set_cparm('X1L_RAD_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X1R_RAD_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X2L_RAD_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X2R_RAD_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X3L_RAD_BOUND', 'BC_PERIODIC') +bhl.config.set_cparm('X3R_RAD_BOUND', 'BC_PERIODIC') + +bhl.config.set_cparm('LOCAL_ANGULAR_DISTRIBUTIONS', True) +bhl.config.set_cparm('LOCAL_ANGLES_NMU', 32) +bhl.config.set_cparm('LOCAL_ANGLES_NX1', NTOT) +bhl.config.set_cparm('LOCAL_ANGLES_NX2', NTOT) +bhl.config.set_cparm('RAD_NUM_TYPES', 4) +bhl.config.set_cparm('NEUTRINO_OSCILLATIONS', True) +bhl.config.set_cparm('FORCE_EQUIPARTITION', False) + + ### RUNTIME PARAMETERS ### + +bhl.config.set_rparm('tf', 'double', default = TF) +bhl.config.set_rparm('dt', 'double', default = T0) +bhl.config.set_rparm('tune_emiss', 'double', default = 100.) +bhl.config.set_rparm('tune_scatt', 'double', default = 20) +bhl.config.set_rparm('t0_tune_emiss', 'double', default = 10.*TF) +bhl.config.set_rparm('t0_tune_scatt', 'double', default = 10.*TF) +bhl.config.set_rparm('L_unit', 'double', default = L_UNIT) +bhl.config.set_rparm('M_unit', 'double', default = M_UNIT) +bhl.config.set_rparm('DTl', 'double', default = DTout) +bhl.config.set_rparm('DTd', 'double', default = DTout) +bhl.config.set_rparm('DTr', 'double', default = 1e6) + +bhl.config.set_rparm('Nsph_tot', 'double', default = Nsph_tot) +bhl.config.set_rparm('Nph_tot', 'double', default = Nph_tot) +bhl.config.set_rparm('frac_e', 'double', default = FRACS[0]) +bhl.config.set_rparm('frac_ebar', 'double', default = FRACS[1]) +bhl.config.set_rparm('frac_x', 'double', default = FRACS[2]) +bhl.config.set_rparm('frac_xbar', 'double', default = FRACS[3]) + +bhl.config.set_rparm('E_MEV', 'double', default = E_MEV) +bhl.config.set_rparm('rho0', 'double', default = RHO0) +bhl.config.set_rparm('uu0', 'double', default = UU0) +bhl.config.set_rparm('eospath', 'string', default = tablepath) +bhl.config.set_rparm('numin', 'double', default = 1.e8) +bhl.config.set_rparm('numax', 'double', default = 1.e30) +bhl.config.set_rparm('direction', 'int', default = IDIM) + + ### CONFIGURE AND COMPILE ### + +bhl.build(PROB) diff --git a/prob/oscillations/problem.c b/prob/oscillations/problem.c new file mode 100644 index 0000000..b0d8049 --- /dev/null +++ b/prob/oscillations/problem.c @@ -0,0 +1,201 @@ +/****************************************************************************** + * * + * PROBLEM.C * + * * + * INITIAL CONDITIONS FOR UNIT TEST FOR MULTISCATT * + * * + ******************************************************************************/ + +#include "decs.h" + +static double Nsph_tot; +static double Nph_tot; + +static double fracs[RAD_NUM_TYPES]; +static double tot_frac; + +static double E_MEV; +static double rho0; +static double uu0; +static int direction; + +// Problem-specific variables to set at runtime +void set_problem_params() { + set_param("Nsph_tot", &Nsph_tot); + set_param("Nph_tot", &Nph_tot); + + // fraction of each species + set_param("frac_e", &fracs[NU_ELECTRON]); + set_param("frac_ebar", &fracs[ANTINU_ELECTRON]); + set_param("frac_x", &fracs[NU_HEAVY]); + set_param("frac_xbar", &fracs[ANTINU_HEAVY]); + + + // neutrino energy + set_param("E_MEV", &E_MEV); + + // background gas + set_param("rho0", &rho0); + set_param("uu0", &uu0); + + set_param("direction", &direction); +} + +double distmu(int type, double x) { + if (type == NU_ELECTRON) { + return 1.0 - (3. / 4.) * x * x; + } else if (type == ANTINU_ELECTRON) { + return (1. / 4.) + (3. / 4.) * x * x; + } else { + return 1.0; + } +} + +// Initialize dynamical variables +void init_prob() { + if ((RADIATION != RADTYPE_NEUTRINOS) || !LOCAL_ANGULAR_DISTRIBUTIONS || + !NEUTRINO_OSCILLATIONS || (RAD_NUM_TYPES != 4)) { + if (mpi_io_proc()) { + printf("Test needs the following physics active!\n" + "\tRADIATION == RADTYPE_NEUTRINOS: %d\n" + "\tLOCAL_ANGULAR_DISTRIBUTIONS: %d\n" + "\tNEUTRINO_OSCILLATIONS: %d\n" + "\tRAD_NUM_TYPES == 4: %d\n", + RADIATION, LOCAL_ANGULAR_DISTRIBUTIONS, NEUTRINO_OSCILLATIONS, + RAD_NUM_TYPES); + } + exit(1); + } + + tot_frac = 0; + TYPELOOP tot_frac += fracs[itp]; + if (fabs(tot_frac - 1) > 1e-5) { + if (mpi_io_proc()) { + printf("Fractions don't add up to 1!\n" + "%.14e %.14e %.14e %.14e\n", + fracs[0], fracs[1], fracs[2], fracs[3]); + } + exit(1); + } + + if ((direction < 1) || (direction > 2)) { + if (mpi_io_proc()) { + printf("direction may only be 1 or 2!\n"); + } + exit(1); + } + + double Nsph_per_cell = Nsph_tot / (N1TOT * N2TOT * N3TOT); + int Nsph_per_thread = (int)ceil(Nsph_per_cell/nthreads); + double nu_cgs = E_MEV * MEV / HPL; + double lnu = log10(nu_cgs); + double nu = pow(10., lnu); + + if (mpi_io_proc()) { + printf("E_MEV = %g\n" + "lnu = %g\n" + "nu = %g\n" + "rho0 = %g\n" + "uu0 = %g\n" + "Nphys = %g\n" + "Nnum = %g\n" + "Nsph_per_cell = %g\n" + "Nsph_per_cell_per_thread = %d\n" + "Fracs = [%g, %g, %g, %g]\n" + "direction = %d\n", + E_MEV, lnu, nu_cgs, rho0, uu0, + Nph_tot, Nsph_tot, Nsph_per_cell, Nsph_per_thread, + fracs[0], fracs[1], fracs[2], fracs[3], + direction); + } + +#pragma omp parallel for collapse(3) + ZLOOP { + P[i][j][k][RHO] = rho0; + P[i][j][k][UU] = uu0; + P[i][j][k][U1] = 0.0; + P[i][j][k][U2] = 0.0; + P[i][j][k][U3] = 0.0; + P[i][j][k][B1] = 0.0; + P[i][j][k][B2] = 0.0; + P[i][j][k][B3] = 0.0; +#if EOS == EOS_TYPE_TABLE + P[i][j][k][YE] = 0.5; + P[i][j][k][YE_EM] = 0.5; +#endif + } // ZLOOP + +#pragma omp parallel + { + double X[NDIM], K_tetrad[NDIM]; + struct of_photon *ph = photon_lists[omp_get_thread_num()]; + ZLOOP { + coord(i, j, k, CENT, X); + + for (int n = 0; n < Nsph_per_thread; n++) { + struct of_photon *phadd = safe_malloc(sizeof(struct of_photon)); + + // photon position centered on cell + phadd->X[2][0] = 0.; + for (int mu = 1; mu < NDIM; mu++) { + phadd->X[2][mu] = X[mu]; + } + + // sample particle type + double x = get_rand(); + double compare = 0; + TYPELOOP { + compare += fracs[itp]; + if (x <= compare) { + phadd->type = itp; + break; + } + } + + // compute mu, phi + double mu, f, fmax = 1.; + do { // rejection sample for mu + mu = 2 * get_rand() - 1; + x = get_rand() * fmax; + f = distmu(phadd->type, mu); + } while (x >= f); + // phi + double phi = 2 * M_PI * get_rand(); + // angles + double costh = mu; + double sinth = sqrt(1 - mu*mu); + + double E = HPL * nu / (ME * CL * CL); + K_tetrad[0] = -E; + if (direction == 1) { + K_tetrad[1] = E * costh; + K_tetrad[2] = E * sinth; + } else { // direction == 2 + K_tetrad[1] = E * sinth; + K_tetrad[2] = E * costh; + } + K_tetrad[3] = 0; + + // Can just use non-moving fluid frame + DLOOP1 { + phadd->Kcov[2][mu] = K_tetrad[mu]; + phadd->Kcon[2][mu] = K_tetrad[mu]; + } + // Minkowski background so Kcov/Kcon only differ by this + phadd->Kcon[2][0] *= -1.; + + phadd->t0 = 0.; + phadd->origin[0] = nstep; + phadd->origin[1] = i; + phadd->origin[2] = j; + phadd->origin[3] = k; + phadd->w = Nph_tot / Nsph_tot; + phadd->nscatt = 0; + + phadd->next = ph; + ph = phadd; + } + } + photon_lists[omp_get_thread_num()] = ph; + } // omp parallel +} diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py index 4493b44..ad2ac0e 100644 --- a/prob/torus_cbc/build.py +++ b/prob/torus_cbc/build.py @@ -115,7 +115,8 @@ if FORTRAN: OPACPATH = "opacity.SFHo.nohoro.juo.brem1.bin" else: - OPACPATH = "opacity.SFHo.nohoro.juo.brem1.h5" + #OPACPATH = "opacity.SFHo.nohoro.juo.brem1.h5" + OPACPATH = "NuLib_rho70_temp62_ye50_ng61_ns4_version1.0_20241120_bhlight.h5" OPACPARAM = "opacbin.LS220.evan.param" OPACPATH = "../../data/"+OPACPATH OPACPARAM = "../../data/"+OPACPARAM @@ -344,8 +345,8 @@ print("Make sure to move the generated table to your run directory!") # Radiation units -NUMIN_MEV = 1.0 -NUMAX_MEV = 100 +NUMIN_MEV = 0.1 +NUMAX_MEV = 500 NUMIN = NUMIN_MEV*cgs['MEV']/cgs['HPL'] NUMAX = NUMAX_MEV*cgs['MEV']/cgs['HPL'] if RADIATION: @@ -406,7 +407,7 @@ bhl.config.set_cparm('KILL_ALL_PACKETS', True) bhl.config.set_cparm('BURROWS_OPACITIES', FORTRAN) bhl.config.set_cparm('HDF5_OPACITIES', HDF) -bhl.config.set_cparm('NU_BINS', 200) +bhl.config.set_cparm('NU_BINS', 61) bhl.config.set_cparm('ESTIMATE_THETAE', False) bhl.config.set_cparm('GRAYABSORPTION', False) bhl.config.set_cparm('BREMSSTRAHLUNG', False) @@ -420,6 +421,12 @@ bhl.config.set_cparm('DIAGNOSTICS_USE_RADTYPES', True) # bhl.config.set_cparm('RECORD_DT_MIN', True) +bhl.config.set_cparm('LOCAL_ANGULAR_DISTRIBUTIONS', True) +bhl.config.set_cparm('LOCAL_ANGLES_NMU', 64) +bhl.config.set_cparm('LOCAL_ANGLES_NX1', 64) +bhl.config.set_cparm('LOCAL_ANGLES_NX2', 64) +bhl.config.set_cparm('RAD_NUM_TYPES', 4) + # Special. Don't turn this on if you don't need to if DIAGNOSTIC: bhl.config.set_cparm("EXIT_ON_INIT", True) diff --git a/script/config.py b/script/config.py index 07729aa..5c97b47 100644 --- a/script/config.py +++ b/script/config.py @@ -42,6 +42,13 @@ def print_config(key, var): def set_cparm(name, value): CPARMS[name] = value +def set_cparm_if_active(name, value): + if util.parm_is_active(CPARMS, name): + print_config(name + ' ', CPARMS[name]) + else: + set_cparm(name, value) + return + # SET RUNTIME PARAMETER. DO NOT OVERWRITE DEFAULT VALUES, AS THIS IS CALLED BY # PROBLEM FILE BEFORE CORE ROUTINE def set_rparm(name, datatype, default=None): @@ -330,8 +337,26 @@ def build(PROBLEM, PATHS): print_config("HDF5_OPACITIES", CPARMS["HDF5_OPACITIES"]) else: set_cparm("HDF5_OPACITIES", 0) + if util.parm_is_active(CPARMS, "RAD_NUM_TYPES"): + print_config("RAD_NUM_TYPES", CPARMS["RAD_NUM_TYPES"]) + else: + if CPARMS['RADIATION'] == 1: + set_cparm("RAD_NUM_TYPES", 1) + else: + set_cparm("RAD_NUM_TYPES", 3) + print_config("RAD_NUM_TYPES", CPARMS["RAD_NUM_TYPES"]) + if util.parm_is_active(CPARMS, "NEUTRINO_OSCILLATIONS"): + print_config("NEUTRINO_OSCILLATIONS ", CPARMS["NEUTRINO_OSCILLATIONS"]) + if util.parm_is_active(CPARMS, 'FORCE_EQUIPARTITION'): + print_config("FORCE_EQUIPARTITION ", CPARMS["FORCE_EQUIPARTITION"]) + else: + set_cparm('FORCE_EQUIPARTITION', 0) + else: + set_cparm('NEUTRINO_OSCILLATIONS', 0) + set_cparm('FORCE_EQUIPARTITION', 0) else: set_cparm("RADIATION", 0) + set_cparm("RAD_NUM_TYPES", 0) if util.parm_is_active(CPARMS, 'ELECTRONS'): print_config("ELECTRONS", CPARMS['ELECTRONS']) else: @@ -390,7 +415,12 @@ def build(PROBLEM, PATHS): if util.parm_is_active(CPARMS, 'RADIATION'): if util.parm_is_active(CPARMS, 'LOCAL_ANGULAR_DISTRIBUTIONS'): - print_config('LOCAL_ANGULAR_DISTRIBUTIONS', CPARMS['LOCAL_ANGULAR_DISTRIBUTIONS']) + print_config('LOCAL_ANGULAR_DISTRIBUTIONS ', + CPARMS['LOCAL_ANGULAR_DISTRIBUTIONS']) + # TODO(JMM): What should defaults be? + set_cparm_if_active('LOCAL_ANGLES_NMU', 64) + set_cparm_if_active('LOCAL_ANGLES_NX1', 64) + set_cparm_if_active('LOCAL_ANGLES_NX2', 64) else: set_cparm('LOCAL_ANGULAR_DISTRIBUTIONS', 0) diff --git a/script/machine/darwin.py b/script/machine/darwin.py index 9321f3a..f99b4f1 100644 --- a/script/machine/darwin.py +++ b/script/machine/darwin.py @@ -31,17 +31,17 @@ def matches_host(): host = os.uname()[1] - return 'darwin-fe' in host + return 'darwin-fe' in host# or 'cn' in host def get_options(): host = {} - local_root = os.path.join(os.environ['HOME'],'local','gcc') + local_root = os.path.join(os.environ['HOME'],'local','skylake-gold') host['NAME'] = os.uname()[1] - host['COMPILER'] = os.path.join(local_root,'hdf5-parallel','bin','h5pcc') + host['COMPILER'] = os.path.join(local_root,'bin','h5pcc') host['COMPILER_FLAGS'] = flags_base + ' ' + '-O3 -march=native' host['DEBUG_FLAGS'] = flags_base + ' ' + '-g -O0' - host['GSL_DIR'] = os.path.join(local_root,'gsl') + host['GSL_DIR'] = os.path.join(local_root,'lib') return host diff --git a/test/oscillations.py b/test/oscillations.py new file mode 100644 index 0000000..62aee9e --- /dev/null +++ b/test/oscillations.py @@ -0,0 +1,175 @@ +################################################################################ +# # +# UNIT TEST FOR Zaizen/Nakamura Oscillations # +# # +################################################################################ + +from __future__ import print_function, division +import os +import sys; sys.dont_write_bytecode = True +sys.path.insert(0, '../script/') +sys.path.insert(0, '../script/analysis') +from subprocess import call +import glob, util, h5py +import numpy as np +import hdf5_to_dict as io +from bhlight import bcall + +TMP_DIR = 'TMP' +util.safe_remove(TMP_DIR) +PROBLEM = 'oscillations' +AUTO = '-auto' in sys.argv +MPI = '-mpi' in sys.argv +gam = 1.4 + +# devnull +try: + from subprocess import DEVNULL # py3k +except ImportError: + import os + DEVNULL = open(os.devnull, 'wb') + +os.chdir('../prob/' + PROBLEM) + +# COMPILE CODE +args = [sys.executable, 'build.py', '-dir', TMP_DIR]#, '-idim', IDIM] +if MPI: + args += ['-mpi'] +print(args) +call(args) +call(['mv', 'sc_eos_gamma_{}.h5'.format(str(gam).replace('.','p')), + TMP_DIR]) +os.chdir('../../test') +call(['mv', '../prob/' + PROBLEM + '/' + TMP_DIR, './']) + +# Since this test is designed to run on a single machine (no batch scripts) +# set openmpi to only use a few threads. Let MPI handle the rest. +if MPI: + import multiprocessing + num_mpi = 8 + num_cpus = multiprocessing.cpu_count() + os.environ['OMP_NUM_THREADS'] = str(int(np.max([2,num_cpus/num_mpi]))) + +# RUN EXECUTABLE +os.chdir(TMP_DIR) +args = ['./bhlight', '-p', 'param_template.dat'] +if MPI: + bcall(args,int(num_mpi),stdin=DEVNULL) +else: + bcall(args) +os.chdir('../') + +# READ SIMULATION OUTPUT +dfiles = np.sort(glob.glob(os.path.join(TMP_DIR,'')+'/dumps/dump*.h5')) +geom = h5py.File(os.path.join(TMP_DIR,'')+'/dumps/grid.h5','r') +dumps = [h5py.File(f, 'r') for f in dfiles] + +b = 0 +mu = geom['local_angles_mu'][:] +angles = [d['local_angles'][b].sum(axis=(0,1)) for d in dumps] +moments = dumps[0]['local_moments'][b].sum(axis=-1).sum(axis=-1) +ntot = [a.sum() for a in angles] +A = moments[0] +B = moments[1] + +Gnu = [d['Gnu'][b].sum(axis=(0,1)) for d in dumps] +Gdiff = (Gnu[1].sum() - Gnu[0].sum()) +Greldiff = 2.0*Gdiff/(np.abs(Gnu[1].sum()) + np.abs(Gnu[0].sum()) + 1e-20) + +Nnu_tot = 1e55 +dmu = mu[1] - mu[0] +fracs = [0.2, 0.4, 0.2, 0.2] +e0 = dmu*Nnu_tot*fracs[0]*(2./3.)*(1 - (3./4.)*mu*mu) +ebar0 = dmu*Nnu_tot*fracs[1]*((1./4.) + (3./4.)*mu*mu) +x0 = dmu*fracs[2]*Nnu_tot*np.ones_like(mu)/2 +xbar0 = dmu*fracs[3]*Nnu_tot*np.ones_like(mu)/2 + +Amask = Gnu[0] < 0 +Bmask = Gnu[0] > 0 +e1 = e0.copy() +e1[Amask] = (1 - (2./3.)*(B/A))*e0[Amask] + (1./3.)*(B/A)*x0[Amask] +e1[Bmask] = (1./3.)*(e0[Bmask] + x0[Bmask]) +x1 = x0.copy() +x1[Amask] = (2./3.)*(B/A)*e0[Amask] + (1 - (1./3.)*(B/A))*x0[Amask] +x1[Bmask] = (2./3.)*(e0[Bmask] + x0[Bmask]) + +if AUTO: + data = {} + data['SOL'] = [0] + data['CODE'] = Gdiff + data['THRESHOLD'] = 0.005 + import pickle + pickle.dump(data, open('data.p', 'wb')) + # clean up + util.safe_remove(TMP_DIR) + sys.exit() + +print("Ntot = ",ntot) +print("Ntot per species (initial) = ", angles[0].sum(axis=-1)) +print("Ntot per species (final) = ", angles[1].sum(axis=-1)) +print("Frac per species (initial) = ", angles[0].sum(axis=-1)/ntot[0]) +print("Frac per species (final) = ", angles[1].sum(axis=-1)/ntot[1]) +print("frac e+x, frac ebar+xbar = ", + (angles[1].sum(axis=-1)/ntot[1])[0]+(angles[1].sum(axis=-1)/ntot[1])[2], + (angles[1].sum(axis=-1)/ntot[1])[1]+(angles[1].sum(axis=-1)/ntot[1])[3]) +print("Gnu start, Gnu end = ",Gnu[0].sum(), Gnu[1].sum()) +print("Relative diff in Gnu = {}".format(Greldiff)) + +# make figure +import matplotlib as mpl; mpl.use('Agg') +from matplotlib import pyplot as plt +mpl.rcParams.update({'font.size':18}) + +fig,axarr = plt.subplots(1,3,sharex=True,figsize=(12,4)) +axarr[0].plot(mu, Gnu[0]/1e53, label='initial') +axarr[0].plot(mu, Gnu[1]/1e53, label='final') +axarr[0].hlines(0,-1.1,1.1,color='k',linestyle='--') +axarr[0].legend() +axarr[0].set_ylabel(r'particle number$/10^{53}$') +axarr[0].set_xlim(-1.01, 1.01) + +cmap = plt.get_cmap("tab10") +for i,ax in enumerate(axarr[1:]): +# ax.plot(mu, angles[i][0]/1e53, +# label=r'$\nu_e$') +# ax.plot(mu, angles[i][1]/1e53, +# #linestyle='--', +# label=r'$\bar{\nu}_e$') +# ax.plot(mu, angles[i][2]/1e53, +# label=r'$\nu_X$') +# ax.plot(mu, angles[i][3]/1e53, +# #linestyle='--', +# label=r'$\bar{\nu}_x$') + ax.set_ylim(0,1.0) +axarr[1].plot(mu, angles[0][0]/1e53, + label=r'$\nu_e$') +axarr[2].plot(mu, angles[1][0]/1e53) + +axarr[1].plot(mu, angles[0][2]/1e53, + label=r'$\nu_x$') +axarr[2].plot(mu, angles[1][2]/1e53) + +axarr[2].plot(mu, e1/1e53, #color=cmap(1), + linestyle=':', label=r'$\nu_e$ analytic') +axarr[2].plot(mu, x1/1e53, #color=cmap(1), + linestyle=':', label=r'$\nu_x$ analytic') +axarr[1].legend(loc='lower center', + ncol=1) +axarr[2].legend(loc='lower center', + ncol=1) + + +for ax in axarr: + ax.set_xlabel(r'$\cos(\theta)$') + +axarr[0].set_title(r'$G_\nu$') +axarr[1].set_title(r'$f(t = 0)$') +axarr[2].set_title(r'$f(t = t_{final})$') + +plt.tight_layout() + +plt.savefig('oscillations_1zone.png', bbox_inches='tight') +plt.savefig('oscillations_1zone.pdf', bbox_inches='tight') + +# clean up +util.safe_remove(TMP_DIR)