From f5cd590612137e47de788eec28ef3576ec9921f7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Mar 2026 09:57:26 -0500 Subject: [PATCH 1/2] modify FFI treatment to use dt / tau --- core/oscillations.c | 10 +++++++++- core/step.c | 11 ++++------- 2 files changed, 13 insertions(+), 8 deletions(-) diff --git a/core/oscillations.c b/core/oscillations.c index d0a9e40..dde8f06 100644 --- a/core/oscillations.c +++ b/core/oscillations.c @@ -187,11 +187,19 @@ 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()]; + double X[NDIM], Kcov[NDIM], Kcon[NDIM]; 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]); + // find local FFI time scale nph already computed because of + // set_Rmunu which must be called each step + int i, j, k; + get_X_K_interp(ph, t, P, X, Kcov, Kcon); + Xtoijk(X, &i, &j, &k); + double tau = 1. / (T_unit * NUFERM * nph[i][j][k] + SMALL); + int b_osc = local_b_osc[ix1][ix2]; int imu = icosth[b_osc]; double A = local_moments[b_osc][MOMENTS_A][ix1][ix2]; @@ -220,7 +228,7 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) { 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 diff --git a/core/step.c b/core/step.c index 983c16a..3e32843 100644 --- a/core/step.c +++ b/core/step.c @@ -74,15 +74,12 @@ void step() { // 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); + printf("\t[Oscillations] tau_ffi = %.14e, dt/tau = %.14e\n", + dt_osc, dt / (dt_osc + SMALL)); } + accumulate_local_angles(); + oscillate(local_moments, Gnu); // check_nu_type("after oscillate"); // DEBUG #endif // OSCILLATIONS #endif From 0b75c9e840f182cecbd25e1014ceb4b221e458e0 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Mar 2026 22:57:27 -0500 Subject: [PATCH 2/2] update with Gails formula --- core/decs.h | 3 ++- core/oscillations.c | 34 +++++++++++++++++++++++++--------- core/step.c | 2 +- 3 files changed, 28 insertions(+), 11 deletions(-) diff --git a/core/decs.h b/core/decs.h index ef37661..5863724 100644 --- a/core/decs.h +++ b/core/decs.h @@ -1041,7 +1041,8 @@ void accumulate_local_angles(); void compute_local_gnu(grid_local_angles_type local_angles, 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(grid_local_angles_type f, + grid_local_moment_type local_moments, grid_Gnu_type gnu); #endif // RAD_NUM_TYPES >= 4 #endif // LOCAL_ANGULAR_DISTRIBUTIONS #endif // RADIATION diff --git a/core/oscillations.c b/core/oscillations.c index dde8f06..1595de9 100644 --- a/core/oscillations.c +++ b/core/oscillations.c @@ -182,7 +182,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(grid_local_angles_type f, + grid_local_moment_type local_moments, grid_Gnu_type gnu) { timer_start(TIMER_OSCILLATIONS); #pragma omp parallel { @@ -193,13 +194,6 @@ void oscillate(grid_local_moment_type local_moments, grid_Gnu_type gnu) { int ix1, ix2, icosth[LOCAL_NUM_BASES]; get_local_angle_bins(ph, &ix1, &ix2, &icosth[0], &icosth[1]); - // find local FFI time scale nph already computed because of - // set_Rmunu which must be called each step - int i, j, k; - get_X_K_interp(ph, t, P, X, Kcov, Kcon); - Xtoijk(X, &i, &j, &k); - double tau = 1. / (T_unit * NUFERM * nph[i][j][k] + SMALL); - int b_osc = local_b_osc[ix1][ix2]; int imu = icosth[b_osc]; double A = local_moments[b_osc][MOMENTS_A][ix1][ix2]; @@ -208,7 +202,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; diff --git a/core/step.c b/core/step.c index 3e32843..5ed6d12 100644 --- a/core/step.c +++ b/core/step.c @@ -79,7 +79,7 @@ void step() { dt_osc, dt / (dt_osc + SMALL)); } accumulate_local_angles(); - oscillate(local_moments, Gnu); + oscillate(local_angles, local_moments, Gnu); // check_nu_type("after oscillate"); // DEBUG #endif // OSCILLATIONS #endif