Skip to content

modify FFI treatment to use dt / tau #50

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
Yurlungur merged 2 commits into from
Apr 2, 2026
Merged

modify FFI treatment to use dt / tau #50

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
f5cd590
modify FFI treatment to use dt / tau
jonahm-LANL Mar 4, 2026
0b75c9e
update with Gails formula
jonahm-LANL Mar 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion core/decs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
30 changes: 27 additions & 3 deletions core/oscillations.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -182,11 +182,13 @@ 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
{
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];
Expand All @@ -200,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;
Expand All @@ -220,7 +244,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
Expand Down
11 changes: 4 additions & 7 deletions core/step.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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_angles, local_moments, Gnu);
// check_nu_type("after oscillate"); // DEBUG
#endif // OSCILLATIONS
#endif
Expand Down