diff --git a/core/constants.h b/core/constants.h index d4fe9b0..5d98507 100644 --- a/core/constants.h +++ b/core/constants.h @@ -32,6 +32,11 @@ #define ROOT2 (1.4142135623730951) #define NUFERM ((ROOT2*GFERM) / HBAR) +// Cos(2 theta21) where theta_m is the vacuum mixing angle theta21 +#define C2TH21 (0.386) +// square of mass gap between m1 and m2 for neutrino masses +#define DM21SQR_EV (7.5e-5) // eV + // Unit conversions #define EV (1.60217653e-12) // Electron-volt #define MEV (1.0e6 * EV) // Mega-Electron-Volt diff --git a/core/decs.h b/core/decs.h index ef37661..d149ab9 100644 --- a/core/decs.h +++ b/core/decs.h @@ -100,12 +100,6 @@ // in scattering stability criterion #define SCATT_BIAS_SAFETY (0.5) // (0.9/RAD_SCATT_TYPES) -//#define NUMIN (1.e10) -//#define NUMAX (1.e25) -//#define NU_BINS_SPEC (200) -//#define NTH (8) -//#define NPHI (8) - // Whether to move polar axis slightly off of coordinate singularity #define COORDSINGFIX 1 #define SINGSMALL (1.E-20) @@ -171,17 +165,20 @@ #define ANTINU_HEAVY (3) #if MULTISCATT_TEST #define RAD_SCATT_TYPES (3) -#else +#else // NOT MULTISCATT TEST #define RAD_SCATT_TYPES (4) // TODO: Should be 5, including electrons #define RSCATT_TYPE_P (0) #define RSCATT_TYPE_N (1) #define RSCATT_TYPE_A (2) #define RSCATT_TYPE_ALPHA (3) #define RSCATT_TYPE_E (4) // TOOD: implement me -#endif +#endif // NOT MULTISCATT TEST #define NRADCOMP (2) #define RADG_YE (4) #define RADG_YE_EM (5) +#define CFI_INACTIVE (0) // For collisional flavor instability +#define CFI_TYPE_GAPPED (1) // Omega_- mode +#define CFI_TYPE_GAPLESS (2) // Omega_+ mode #elif RADIATION == RADTYPE_LIGHT #define RAD_SCATT_TYPES (1) #define NRADCOMP (0) @@ -342,9 +339,9 @@ extern grid_int_type Nsph; extern grid_double_type nph; #if RZ_HISTOGRAMS extern rz_hist_type rz_r_orig_hist, rz_z_orig_hist; -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI extern rz_hist_type osc_rz_r_orig_hist, osc_rz_z_orig_hist; -#endif // NEUTRINO_OSCILLATIONS +#endif // NEUTRINO_OSCILLATIONS_FFI #endif // RZ_HISTOGRAMS extern struct of_photon **photon_lists; @@ -396,13 +393,31 @@ typedef double grid_local_moment_type[LOCAL_NUM_BASES][LOCAL_NUM_MOMENTS] [LOCAL_ANGLES_NX1][LOCAL_ANGLES_NX2]; typedef int grid_local_basis_idx_type[LOCAL_ANGLES_NX1][LOCAL_ANGLES_NX2]; typedef double grid_local_count_type[LOCAL_ANGLES_NX1][LOCAL_ANGLES_NX2]; -extern grid_Gnu_type Gnu, local_Ns, local_wsqr; -extern grid_local_moment_type local_moments; +extern grid_Gnu_type Gnu, local_Ns, local_wsqr; +extern grid_local_moment_type local_moments; extern grid_local_basis_idx_type local_b_osc; -extern grid_local_count_type local_osc_count; +extern grid_local_count_type local_osc_count; #endif // #if RAD_NUM_TYPES >= 4 #endif // LOCAL_ANGULAR_DISTRIBUTIONS +#if NEUTRINO_OSCILLATIONS_CFI +typedef double grid_symm_double_type[N1 + 2 * NG][N2 + 2 * NG]; +typedef int grid_symm_int_type[N1 + 2 * NG][N2 + 2 * NG]; +typedef double grid_symm_radtype_type[N1 + 2 * NG][N2 + 2 * NG][RAD_NUM_TYPES]; +typedef double grid_CFI_Gamma_type[N1 + 2 * NG][N2 + 2 * NG][2]; + +// number of neutrinos per flavor +extern grid_symm_radtype_type nph_flavor; +// individual distribution function-weighted opacity averages +extern grid_symm_radtype_type kappa_avg; +extern grid_CFI_Gamma_type cfi_Gamma; + +// Which CFI mode is active, if any +extern grid_symm_int_type cfi_active_mode; +// time scale for asymptotic state +extern grid_symm_double_type cfi_tau_asymp; +#endif // NEUTRINO_OSCILLATIONS_CFI + #endif // RADIATION // Default initialization is 0, which in this case is @@ -582,11 +597,11 @@ struct of_photon { // radiation type. For neutrinos, flavor. Always active. // Not always important. // TODO: make sure to always set type when it's needed. - int type; - int nscatt; - int origin[NDIM]; - double t0; - int is_tracked; + int type; + int nscatt; + int origin[NDIM]; + double t0; + int is_tracked; // Only relevant for neutrino oscillations int osc_count; struct of_photon *next; @@ -711,8 +726,8 @@ 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) \ +#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 @@ -1033,26 +1048,35 @@ double alpha_nu_hdf(double nu, int type, const struct of_microphysics *m); // oscillations.c #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS -double get_dt_oscillations(); +// FFI +double get_dt_ffi(); void get_local_angle_bins( - struct of_photon *ph, int *pi, int *pj, int *pmu1, int *pmu2); + 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_Ns, grid_Gnu_type local_wsqr, grid_Gnu_type gnu); + grid_Gnu_type local_Ns, grid_Gnu_type local_wsqr, 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); +void oscillate_ffi(grid_local_angles_type f, + grid_local_moment_type local_moments, grid_Gnu_type gnu); #endif // RAD_NUM_TYPES >= 4 #endif // LOCAL_ANGULAR_DISTRIBUTIONS + +// CFI +#if NEUTRINO_OSCILLATIONS_CFI +void compute_cfi_symmetrized_avgs(grid_symm_radtype_type nph_flavor, + grid_symm_radtype_type kappa_avg, + grid_CFI_Gamma_type cfi_gamma); +void compute_cfi_active_mode(grid_int_type cfi_active_mode); +#endif // NEUTRINO_OSCILLATIONS_CFI + #endif // RADIATION // passive.c -//#if NVAR_PASSIVE > 0 void fixup_passive( int i, int j, int k, double pv[NVAR], double pv_prefloor[NVAR]); void init_passives(); void name_passives(); -//#endif // phys.c void primtoflux(double *pr, struct of_state *q, int dir, int magnetic, @@ -1129,6 +1153,7 @@ 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); +int is_antiparticle(const struct of_photon *ph); #endif // NEUTRINOS #endif // RADIATION @@ -1235,7 +1260,7 @@ double find_median(double *array, int size); double interp_1d(double x, const double xmin, const double xmax, const int imin, const int imax, const double *restrict tab_x, const double *restrict tab_y); int find_index(double value, const double *array, int size); -void * safe_malloc(size_t size); +void *safe_malloc(size_t size); void safe_system(const char *command); void safe_fscanf(FILE *stream, const char *format, ...); int is_practically_nan(double v); diff --git a/core/defs.h b/core/defs.h index 6dd006b..eef677a 100644 --- a/core/defs.h +++ b/core/defs.h @@ -29,9 +29,9 @@ grid_int_type Nsph; grid_double_type nph; #if RZ_HISTOGRAMS rz_hist_type rz_r_orig_hist, rz_z_orig_hist; -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI rz_hist_type osc_rz_r_orig_hist, osc_rz_z_orig_hist; -#endif // NEUTRINO_OSCILLATIONS +#endif // NEUTRINO_OSCILLATIONS_FFI #endif // RZ_HISTOGRAMS struct of_photon **photon_lists; @@ -67,6 +67,20 @@ grid_local_basis_idx_type local_b_osc; grid_local_count_type local_osc_count; #endif // RAD_NUM_TYPES #endif // LOCAL_ANGULAR_DISTRIBUTIONS + +#if NEUTRINO_OSCILLATIONS_CFI +// number of neutrinos per flavor +grid_symm_radtype_type nph_flavor; +// individual distribution function-weighted opacity averages +grid_symm_radtype_type kappa_avg; +grid_CFI_Gamma_type cfi_Gamma; + +// Which CFI mode is active, if any +grid_symm_int_type cfi_active_mode; +// time scale for asymptotic state +grid_symm_double_type cfi_tau_asymp; +#endif // NEUTRINO_OSCILLATIONS_CFI + #endif // RADIATION #if ELECTRONS diff --git a/core/diag.c b/core/diag.c index 10c4da6..416ef2e 100644 --- a/core/diag.c +++ b/core/diag.c @@ -38,14 +38,14 @@ void reset_dump_variables() { #if RZ_HISTOGRAMS memset(rz_r_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double)); memset(rz_z_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double)); -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI memset(osc_rz_r_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double)); memset(osc_rz_z_orig_hist, 0, RZ_HISTOGRAMS_N * sizeof(double)); -#endif // NEUTRINO_OSCILLATIONS +#endif // NEUTRINO_OSCILLATIONS_FFI #endif // RZ_HISTOGRAMS -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI memset(local_osc_count, 0, LOCAL_ANGLES_NX1*LOCAL_ANGLES_NX2*sizeof(double)); -#endif // NEUTRINO_OSCILLATIONS +#endif // NEUTRINO_OSCILLATIONS_FFI #endif // RADIATION } @@ -310,7 +310,7 @@ void diag(int call_code) { #if ELECTRONS fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_ELECTRON)); #endif -#if NEUTRINO_OSCILLATIONS || LOCAL_ANGULAR_DISTRIBUTIONS +#if NEUTRINO_OSCILLATIONS_FFI || LOCAL_ANGULAR_DISTRIBUTIONS fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_OSCILLATIONS)); #endif fprintf(ener_file, "\n"); @@ -596,21 +596,21 @@ void generate_rz_histograms() { #pragma omp atomic rz_z_orig_hist[iz] += ph->w; -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI if (ph->osc_count) { #pragma omp atomic osc_rz_r_orig_hist[ir] += ph->w; #pragma omp atomic osc_rz_z_orig_hist[iz] += ph->w; } -#endif // NEUTRINO_OSCILLATIONS +#endif // NEUTRINO_OSCILLATIONS_FFI ph = ph->next; } } mpi_dbl_allreduce_array((double *)rz_r_orig_hist, RZ_HISTOGRAMS_N); mpi_dbl_allreduce_array((double *)rz_z_orig_hist, RZ_HISTOGRAMS_N); -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI mpi_dbl_allreduce_array((double *)osc_rz_r_orig_hist, RZ_HISTOGRAMS_N); mpi_dbl_allreduce_array((double *)osc_rz_z_orig_hist, RZ_HISTOGRAMS_N); #endif diff --git a/core/emissivity.c b/core/emissivity.c index 0ce6e26..8c1166e 100644 --- a/core/emissivity.c +++ b/core/emissivity.c @@ -46,6 +46,7 @@ void init_emissivity() { #endif } +// todo(JMM): Pass in a rescaling option double jnu(double nu, int type, const struct of_microphysics *m, double theta) { double jnu = 0.; @@ -71,6 +72,7 @@ double jnu(double nu, int type, const struct of_microphysics *m, double theta) { return jnu; } +// todo(JMM): Pass in a rescaling option double Jnu(double nu, int type, const struct of_microphysics *m) { double Jnu = 0.; @@ -96,6 +98,7 @@ double Jnu(double nu, int type, const struct of_microphysics *m) { return Jnu; } +// todo(JMM): Pass in a rescaling option double integrandJ(double x, void *params) { struct of_J_params * p = (struct of_J_params *)params; struct of_microphysics *m = p->microphysics; @@ -110,6 +113,7 @@ double integrandJ(double x, void *params) { return Jsamp; } +// todo(JMM): Pass in a rescaling option double get_J(struct of_microphysics *m) { #if RADIATION == RADTYPE_LIGHT { diff --git a/core/io.c b/core/io.c index b6aa310..9e03cad 100644 --- a/core/io.c +++ b/core/io.c @@ -878,8 +878,8 @@ void dump() { WRITE_HDR(numin, TYPE_DBL); WRITE_HDR(numax, TYPE_DBL); - int neutrino_oscillations = NEUTRINO_OSCILLATIONS; - WRITE_HDR(neutrino_oscillations, TYPE_INT); + int neutrino_oscillations_ffi = NEUTRINO_OSCILLATIONS_FFI; + WRITE_HDR(neutrino_oscillations_ffi, TYPE_INT); int force_equipartition = FORCE_EQUIPARTITION; WRITE_HDR(force_equipartition, TYPE_INT); int local_angular_distributions = LOCAL_ANGULAR_DISTRIBUTIONS; @@ -1145,7 +1145,7 @@ void dump() { #undef RANK } -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI { mpi_dbl_allreduce_array( (double *)local_osc_count, LOCAL_ANGLES_NX1 * LOCAL_ANGLES_NX2); @@ -1169,7 +1169,7 @@ void dump() { TYPE_DBL); #undef RANK } -#endif // NEUTRINO_OSCILLATIONS +#endif // NEUTRINO_OSCILLATIONS_FFI #endif // LOCAL_ANGULAR_DISTRIBUTIONS #if RZ_HISTOGRAMS @@ -1186,7 +1186,7 @@ void dump() { } WRITE_ARRAY(rz_r_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); WRITE_ARRAY(rz_z_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); -#if NEUTRINO_OSCILLATIONS +#if NEUTRINO_OSCILLATIONS_FFI WRITE_ARRAY(osc_rz_r_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); WRITE_ARRAY(osc_rz_z_orig_hist, RANK, fdims, fstart, fcount, mdims, mstart, TYPE_DBL); #endif diff --git a/core/make_superphotons.c b/core/make_superphotons.c index 4f77cd7..8eaa734 100644 --- a/core/make_superphotons.c +++ b/core/make_superphotons.c @@ -153,6 +153,7 @@ void sample_photon(int i, int j, int k, double t, double dt, int type, double weight = get_wgt(nu, dtau); // Sample emissivity in solid angle + // todo(JMM): Pass in a rescaling option double jmax = jnu(nu, type, m, 0.5 * M_PI); do { cth[0] = 2. * get_rand() - 1.; @@ -288,6 +289,7 @@ void get_dndlnu(int i, int j, int k, double dt, double dndlnu[NU_BINS + 1], double dndlnu_max = -1.e100; for (int n = 0; n <= NU_BINS; n++) { + // todo(JMM): Pass in a rescaling option double Jsamp = Jnu(nusamp[n], type, m); Jsamp *= dx[1] * dx[2] * dx[3] * pow(L_unit, 3.) * ggeom[i][j][CENT].g; @@ -330,6 +332,7 @@ void set_weight(grid_prim_type Prad, grid_eosvar_type extra) { get_fluid_zone(i, j, k, Prad, extra, &m, Ucon, Ucov, Bcon, Bcov); TYPELOOP { for (int n = 0; n <= NU_BINS; n++) { + // todo(JMM): Pass in a rescaling option Jtot += Jnu(nusamp[n], itp, &m) * zoneVol * ggeom[i][j][CENT].g; } } // TYPELOOP @@ -353,6 +356,7 @@ double f(double x, void *params) { int type = p->type; double nu = exp(x); + // todo(JMM): Pass in a rescaling option double Jsamp = Jnu(nu, type, m) * nu; double wgt = get_wgt(nu, get_dtau(nu, type, dt, m)); if (isinf(Jsamp) || wgt < SMALL) { diff --git a/core/oscillations.c b/core/oscillations.c index d0a9e40..5b1c23b 100644 --- a/core/oscillations.c +++ b/core/oscillations.c @@ -10,7 +10,7 @@ #if RADIATION == RADTYPE_NEUTRINOS #if LOCAL_ANGULAR_DISTRIBUTIONS -double get_dt_oscillations() { +double get_dt_ffi() { timer_start(TIMER_OSCILLATIONS); set_Rmunu(); // So we have Nsph and nph double nph_max = 0; @@ -51,7 +51,7 @@ void accumulate_local_angles() { #pragma omp atomic local_Ns[b][ix1][ix2][icosth[b]] += 1.; #pragma omp atomic - local_wsqr[b][ix1][ix2][icosth[b]] += (ph->w)*(ph->w); + local_wsqr[b][ix1][ix2][icosth[b]] += (ph->w) * (ph->w); } } ph = ph->next; @@ -106,7 +106,7 @@ void get_local_angle_bins( 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( @@ -151,7 +151,8 @@ void compute_local_gnu(grid_local_angles_type f, grid_Gnu_type local_Ns, // 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) { +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 { @@ -182,7 +183,8 @@ 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) { +void oscillate_ffi(grid_local_angles_type f, + grid_local_moment_type local_moments, grid_Gnu_type gnu) { timer_start(TIMER_OSCILLATIONS); #pragma omp parallel { @@ -190,6 +192,7 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) { while (ph != NULL) { if (ph->type != TYPE_TRACER) { int ix1, ix2, icosth[LOCAL_NUM_BASES]; + double X[NDIM], Kcov[NDIM], Kcon[NDIM]; get_local_angle_bins(ph, &ix1, &ix2, &icosth[0], &icosth[1]); int b_osc = local_b_osc[ix1][ix2]; @@ -200,7 +203,29 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) { // gnu == 0 when we activated stddev trigger. Don't oscillate. if (((G != 0) || FORCE_EQUIPARTITION) && (A != 0) && (B != 0)) { - // if ((A != 0) && (B != 0)) { + // finite volumes cell + int i, j, k; + get_X_K_interp(ph, t, P, X, Kcov, Kcon); + Xtoijk(X, &i, &j, &k); + + // compute crossing depth + double n_tot_supercell = 0; + TYPELOOP { + LOCALMULOOP { + // total neutrinos of all flavors + n_tot_supercell += fabs(f[0][ix1][ix2][itp][imu]); + } + } + // crossing depth in units of number/cm^3 + // sqrt(AB)/n_tot_supercell is normalized unitless crossing + // depth but units should be in number/cm^3. + double crossing_depth = sqrt(fabs(A * B))*nph[i][j][k] / (n_tot_supercell + SMALL); + + // find local FFI time scale nph already computed because of + // Rmunu This crossing depth formula is slightly more + // sophisticated than just 1 / G_F n + double tau = 1. / (T_unit * NUFERM * crossing_depth + SMALL); + // If A == B then which region we treat as shallow is // unimportant. Psurvive = 1/3 for both regions. int A_is_shallow = A < B; @@ -213,14 +238,14 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) { int in_shallow = (A_is_shallow && g_in_A) || (B_is_shallow && g_in_B); - double peq = nu_is_heavy(ph->type) ? (2./3.) : (1./3.); + 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) * shallow / (deep + SMALL)); #endif // FORCE_EQUIPARTITION - double p_osc = 1. - p_survival; + double p_osc = MY_MIN(1.0, dt / (tau + SMALL)) * (1. - p_survival); if (get_rand() < p_osc) { // JMM: // Type order is NUE, NUEBAR, NUX, NUXBAR @@ -241,6 +266,91 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) { } #endif // RAD_NUM_TYPES >= 4 - #endif // LOCAL_ANGULAR_DISTRIBUTIONS + +#if NEUTRINO_OSCILLATIONS_CFI +void compute_cfi_symmetrized_avgs(grid_symm_radtype_type nph_flavor, + grid_symm_radtype_type kappa_avg, + grid_CFI_Gamma_type cfi_gamma) { + timer_start(TIMER_OSCILLATIONS); + + static const int SYMM_SIZE = (N1 + 2*NG)*(N2 + 2*NG); + + memset(nph_flavor, 0, SYMM_SIZE * RAD_NUM_TYPES * sizeof(double)); + memset(kappa_avg, 0, SYMM_SIZE * RAD_NUM_TYPES * sizeof(double)); + memset(cfi_Gamma, 0, 2 * SYMM_SIZE * sizeof(double)); + +#pragma omp parallel + { + struct of_photon *ph = photon_lists[omp_get_thread_num()]; + while (ph != NULL) { + if (ph->type != TYPE_TRACER) { + double X[NDIM], Kcov[NDIM], Kcon[NDIM]; + double kappa[RAD_NUM_TYPES]; + int i, j, k; + get_X_K_interp(ph, t, P, X, Kcov, Kcon); + Xtoijk(X, &i, &j, &k); + + #pragma omp atomic + nph_flavor[i][j][ph->type] += ph->w; + + double theta = get_bk_angle( + X, Kcon, Ucov_grd[i][j][k], + Bcov_grd[i][j][k], m_grd[i][j][k].B); + double nu = get_fluid_nu(X, Kcov, Ucon_grd[i][j][k]); + TYPELOOP { + kappa[itp] = alpha_inv_abs(nu, itp, + &(m_grd[i][j][k]), theta); + kappa[itp] /= (nu + SMALL); + } + + #pragma omp atomic + kappa_avg[i][j][ph->type] += kappa[ph->type] * ph->w; + + // electrons and antis positive, heavies and antis negative + int gamma_sign = (ph->type < NU_HEAVY) ? 1 : -1; + + // Check for the need to divide neutrino count for heavies by 2 + double heavy_dupfac = ((RAD_NUM_TYPES < 4) && (ph->type >= NU_HEAVY))? 0.5 : 1.0; + + // check for antiparticles + int isanti = is_antiparticle(ph); + + int kappa_e_idx = NU_ELECTRON + isanti; // kappa_e for es, kappa_ebar for ebars + // kappa_x for xs, kappa_xbar for xbars. But account for x = xbar + int kappa_x_idx = NU_HEAVY + (RAD_NUM_TYPES == 4)*isanti; + + // es contribute + to this sum, xs minus. Need to divide by 2 + // if heavies and their antis are bundled together. + double weight = gamma_sign * heavy_dupfac * ph->w; + + // gamma_0 = (1/2) (kappa_e + kappa_x) + // gamma_1 = (1/2) (kappa_ebar + kappa_xbar) + // + // Method A + // avg gamma_0 = \int E^2 dE (kappa_e + kappa_x) (f_e - f_x) / \int E^2 / \int E^2 dE (f_e - f_x) + // avg gamma_1 = \int E^2 dE (kappa_ebar + kappa_xbar) (f_ebar - f_xbar) / \int E^2 / \int E^2 dE (f_ebar - f_xbar) + #pragma omp atomic + cfi_Gamma[i][j][isanti] += 0.5*(kappa[kappa_e_idx] + kappa[kappa_x_idx])*weight; + } + ph = ph->next; + } + } + + mpi_dbl_allreduce_array((double *)nph_flavor, SYMM_SIZE * RAD_NUM_TYPES); + mpi_dbl_allreduce_array((double *)kappa_avg, SYMM_SIZE * RAD_NUM_TYPES); + mpi_dbl_allreduce_array((double *)cfi_Gamma, 2 * SYMM_SIZE); + + #pragma omp parallel for collapse(2) + ILOOP { + JLOOP { + cfi_Gamma[i][j][0] /= (nph_flavor[i][j][NU_ELECTRON] - nph_flavor[i][j][NU_HEAVY]); + cfi_Gamma[i][j][1] /= (nph_flavor[i][j][ANTINU_ELECTRON] - nph_flavor[i][j][ANTINU_HEAVY]); + } + } + + timer_stop(TIMER_OSCILLATIONS); +} +#endif // NEUTRINO_OSCILLATIONS_CFI + #endif // RADIATION diff --git a/core/rad_utils.c b/core/rad_utils.c index eddcfe5..9bbf6da 100644 --- a/core/rad_utils.c +++ b/core/rad_utils.c @@ -749,7 +749,7 @@ double get_min_dt_cool(grid_prim_type P, grid_eosvar_type extra) { } #if RADIATION == RADTYPE_NEUTRINOS -void record_lepton_flux(const struct of_photon *ph) { +void record_lepton_flux(const struct of_photon *ph) { #pragma omp atomic lepton_lost_local += (ph->w) * get_lepton_sign(ph); } @@ -765,6 +765,14 @@ int nu_is_heavy(const int radtype) { return ((radtype == NU_HEAVY) || (radtype == ANTINU_HEAVY)); } +int is_antiparticle(const struct of_photon *ph) { +#if RAD_NUM_TYPES == 4 + return (ph->type == ANTINU_ELECTRON) || (ph->type == ANTINU_HEAVY); +#else + return (ph->type == ANTINU_ELECTRON) || (ph->type == NU_HEAVY); +#endif +} + // for debugging void check_nu_type(const char *location) { #pragma omp parallel diff --git a/core/radiation.c b/core/radiation.c index aa56f22..d324236 100644 --- a/core/radiation.c +++ b/core/radiation.c @@ -9,24 +9,35 @@ #include "decs.h" #if RADIATION -double Bnu_inv(double nu, - const struct of_microphysics - *m) { // unused for neutrinos. Use tables by Burrows et al. -#if RADIATION == RADTYPE_LIGHT +double Bnu_inv(double nu, const struct of_microphysics *m) { double x; +#if RADIATION == RADTYPE_LIGHT x = HPL * nu / (ME * CL * CL * m->Thetae); +#elif RADIATION == RADTYPE_NEUTRINOS + x = HPL * nu / (KBOL * m->T); +#endif +#if RADIATION == RADTYPE_LIGHT + // Planck (i.e., Bose-Einstein) distribution if (x < 1.e-3) // Taylor-expand small arguments for numerical accuracy return ( (2. * HPL / (CL * CL)) / (x / 24. * (24. + x * (12. + x * (4. + x))))); else return ((2. * HPL / (CL * CL)) / (exp(x) - 1.)); -#endif -#if RADIATION == RADTYPE_NEUTRINOS - return 0.; // STUB +#elif RADIATION == RADTYPE_NEUTRINOS + // Fermi-Dirac distribution + double dist; + if (x < 1.e-3) { // Taylor-expand small arguments for numerical accuracy + double x2 = x*x; + dist = 0.5*(1 - 0.5*x*(1 - (x2/12)*(1 - 0.1*x2))); + } else { + dist = 1. / (exp(x) + 1.); + } + return (2. * HPL / (CL * CL)) * dist; #endif } +// todo(JMM): Pass in a rescaling option double jnu_inv( double nu, int type, const struct of_microphysics *m, double theta) { double j; @@ -47,6 +58,7 @@ double alpha_inv_scatt( } // Invariant absorption opacity +// TODO(JMM): Need to pass in a rescaling option double alpha_inv_abs( double nu, int type, const struct of_microphysics *m, double theta) { #if RADIATION == RADTYPE_NEUTRINOS diff --git a/core/step.c b/core/step.c index 983c16a..25cf4a8 100644 --- a/core/step.c +++ b/core/step.c @@ -28,10 +28,10 @@ void step() { #if RADIATION double dt_cool; #if RADIATION == RADTYPE_NEUTRINOS && LOCAL_ANGULAR_DISTRIBUTIONS && \ - RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS + RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS_FFI // not used for timestep control. Used to turn oscillations on or // off. - double dt_osc = get_dt_oscillations(); + double dt_ffi = get_dt_ffi(); #endif // oscillations #endif // radiation dtsave = dt; @@ -73,16 +73,13 @@ void step() { bound_superphotons(Ph, t, dt); // 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); + RAD_NUM_TYPES >= 4 && NEUTRINO_OSCILLATIONS_FFI 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); + printf("\t[Oscillations] tau_ffi = %.14e, dt/tau_ffi = %.14e\n", + dt_ffi, dt / (dt_ffi + SMALL)); } + accumulate_local_angles(); + oscillate_ffi(local_angles, local_moments, Gnu); // check_nu_type("after oscillate"); // DEBUG #endif // OSCILLATIONS #endif diff --git a/core/timing.c b/core/timing.c index 6a5b3c3..00296a7 100644 --- a/core/timing.c +++ b/core/timing.c @@ -91,7 +91,7 @@ void report_performance() { fprintf(stdout, " INTERACT: %8.4g s (%.4g %%)\n", times[TIMER_INTERACT] / dnstep, 100. * times[TIMER_INTERACT] / times[TIMER_ALL]); -#if NEUTRINO_OSCILLATIONS || LOCAL_ANGULAR_DISTRIBUTIONS +#if NEUTRINO_OSCILLATIONS_FFI || LOCAL_ANGULAR_DISTRIBUTIONS fprintf(stdout, " OSCILL: %8.4g s (%.4g %%)\n", times[TIMER_OSCILLATIONS] / dnstep, 100. * times[TIMER_OSCILLATIONS] / times[TIMER_ALL]); diff --git a/prob/oscillations/build.py b/prob/ffi/build.py similarity index 99% rename from prob/oscillations/build.py rename to prob/ffi/build.py index ac6367b..a4c23e2 100644 --- a/prob/oscillations/build.py +++ b/prob/ffi/build.py @@ -129,7 +129,7 @@ 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('NEUTRINO_OSCILLATIONS_FFI', True) bhl.config.set_cparm('FORCE_EQUIPARTITION', False) bhl.config.set_cparm("RZ_HISTOGRAMS", True) diff --git a/prob/oscillations/problem.c b/prob/ffi/problem.c similarity index 97% rename from prob/oscillations/problem.c rename to prob/ffi/problem.c index b0d8049..198ded1 100644 --- a/prob/oscillations/problem.c +++ b/prob/ffi/problem.c @@ -54,14 +54,14 @@ double distmu(int type, double x) { // Initialize dynamical variables void init_prob() { if ((RADIATION != RADTYPE_NEUTRINOS) || !LOCAL_ANGULAR_DISTRIBUTIONS || - !NEUTRINO_OSCILLATIONS || (RAD_NUM_TYPES != 4)) { + !NEUTRINO_OSCILLATIONS_FFI || (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" + "\tNEUTRINO_OSCILLATIONS_FFI: %d\n" "\tRAD_NUM_TYPES == 4: %d\n", - RADIATION, LOCAL_ANGULAR_DISTRIBUTIONS, NEUTRINO_OSCILLATIONS, + RADIATION, LOCAL_ANGULAR_DISTRIBUTIONS, NEUTRINO_OSCILLATIONS_FFI, RAD_NUM_TYPES); } exit(1); diff --git a/prob/torus_cbc/build.py b/prob/torus_cbc/build.py index b39e2d5..1ac3561 100644 --- a/prob/torus_cbc/build.py +++ b/prob/torus_cbc/build.py @@ -40,7 +40,7 @@ TRACERTEST = '-tracertest' in sys.argv RESTARTTEST = '-restarttest' in sys.argv HDF = '-hdf' in sys.argv -OSCILLATIONS = "-oscillations" in sys.argv +FFI = "-ffi" in sys.argv N1N2N3CPU_FROM_CLI = '-n1n2n3cpu' in sys.argv N1N2N3TOT_FROM_CLI = '-n1n2n3tot' in sys.argv @@ -116,7 +116,7 @@ if FORTRAN: OPACPATH = "opacity.SFHo.nohoro.juo.brem1.bin" else: - if OSCILLATIONS: + if FFI: OPACPATH = "NuLib_rho70_temp62_ye50_ng61_ns4_version1.0_20241120_bhlight.h5" else: OPACPATH = "opacity.SFHo.nohoro.juo.brem1.h5" @@ -426,13 +426,13 @@ bhl.config.set_cparm('RZ_HISTOGRAMS', True) # bhl.config.set_cparm('RECORD_DT_MIN', True) -if OSCILLATIONS: +if FFI: bhl.config.set_cparm('LOCAL_ANGULAR_DISTRIBUTIONS', True) bhl.config.set_cparm('LOCAL_ANGLES_NMU', 32) bhl.config.set_cparm('LOCAL_ANGLES_NX1', 32) bhl.config.set_cparm('LOCAL_ANGLES_NX2', 32) bhl.config.set_cparm('RAD_NUM_TYPES', 4) - bhl.config.set_cparm('NEUTRINO_OSCILLATIONS', True) + bhl.config.set_cparm('NEUTRINO_OSCILLATIONS_FFI', True) # Special. Don't turn this on if you don't need to if DIAGNOSTIC: diff --git a/script/config.py b/script/config.py index 892091d..91ea80f 100644 --- a/script/config.py +++ b/script/config.py @@ -345,15 +345,19 @@ def build(PROBLEM, PATHS): 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, "NEUTRINO_OSCILLATIONS_FFI"): + print_config("NEUTRINO_OSCILLATIONS_FFI ", CPARMS["NEUTRINO_OSCILLATIONS_FFI"]) 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('NEUTRINO_OSCILLATIONS_FFI', 0) set_cparm('FORCE_EQUIPARTITION', 0) + if util.parm_is_active(CPARMS, "NEUTRINO_OSCILLATIONS_CFI"): + print_config("NEUTRINO_OSCILLATIONS_CFI ", CPARMS["NEUTRINO_OSCILLATIONS_CFI"]) + else: + set_cparm('NEUTRINO_OSCILLATIONS_CFI', 0) if util.parm_is_active(CPARMS, "RZ_HISTOGRAMS"): print_config("RZ_HISTOGRAMS", CPARMS["RZ_HISTOGRAMS"]) if not util.parm_is_active(CPARMS, "RZ_HISTOGRAMS_N"): diff --git a/test/oscillations.py b/test/ffi.py similarity index 97% rename from test/oscillations.py rename to test/ffi.py index 2da64b3..2df2766 100644 --- a/test/oscillations.py +++ b/test/ffi.py @@ -17,7 +17,7 @@ TMP_DIR = 'TMP' util.safe_remove(TMP_DIR) -PROBLEM = 'oscillations' +PROBLEM = 'ffi' AUTO = '-auto' in sys.argv MPI = '-mpi' in sys.argv gam = 1.4 @@ -175,8 +175,8 @@ plt.tight_layout() -plt.savefig('oscillations_1zone.png', bbox_inches='tight') -plt.savefig('oscillations_1zone.pdf', bbox_inches='tight') +plt.savefig('ffi_1zone.png', bbox_inches='tight') +plt.savefig('ffi_1zone.pdf', bbox_inches='tight') # clean up util.safe_remove(TMP_DIR) diff --git a/test/test_all.py b/test/test_all.py index 5784054..aa51e2a 100644 --- a/test/test_all.py +++ b/test/test_all.py @@ -45,7 +45,7 @@ M_NU_TESTS = ['yedecay.py', 'yedecay.py -antinu', 'yedecay.py -mpi', # 'yedecay.py -mpi -antinu', 'multiscatt.py', - 'tracers1d.py'] + 'tracers1d.py', 'ffi.py'] M_LIGHT_TESTS = ['binning.py', 'brem.py', 'thermalization.py', 'thermalization_mpi.py', 'comptonization.py']