Skip to content

[WIP] CFI #51

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

Open
Yurlungur wants to merge 17 commits into
base: master
Choose a base branch
from
Open

[WIP] CFI #51

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
17 commits
Select commit Hold shift + click to select a range
f5cd590
modify FFI treatment to use dt / tau
jonahm-LANL Mar 4, 2026
d0b0eb5
clarify names for ffi vs CFI
jonahm-LANL Mar 4, 2026
93163f5
rename more things to ffi to distinguish from cfi
jonahm-LANL Mar 4, 2026
06a9a55
add 1zone ffi test to automated suite
jonahm-LANL Mar 4, 2026
ec66cdc
in progress
jonahm-LANL Mar 5, 2026
56db136
argh
jonahm-LANL Mar 5, 2026
0b75c9e
update with Gails formula
jonahm-LANL Mar 5, 2026
9ad6159
merge jmm/ffi-timescale
jonahm-LANL Mar 5, 2026
fa0f575
oops bad merge fixed now
jonahm-LANL Mar 5, 2026
cddc621
cleanup
jonahm-LANL Mar 5, 2026
5e6e2a2
cleanup
jonahm-LANL Mar 5, 2026
17c5607
rename oscillations -> ffi
jonahm-LANL Mar 5, 2026
cf3bd2c
CFI correct name.
jonahm-LANL Mar 5, 2026
f2470ad
add some todos for this project
jonahm-LANL Mar 5, 2026
cb0f94f
start building the averaging machinery
jonahm-LANL Mar 5, 2026
1e27120
compute average gammas and kappas for growth rate calc
jonahm-LANL Mar 5, 2026
eed23c2
Bnu inv now supports Fermi-Dirac
jonahm-LANL Mar 6, 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
5 changes: 5 additions & 0 deletions core/constants.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
79 changes: 52 additions & 27 deletions core/decs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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);
Comment thread
kelslund marked this conversation as resolved.
#endif // NEUTRINOS
#endif // RADIATION

Expand Down Expand Up @@ -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);
Expand Down
18 changes: 16 additions & 2 deletions core/defs.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand Down
16 changes: 8 additions & 8 deletions core/diag.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
}

Expand Down Expand Up @@ -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");
Expand Down Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions core/emissivity.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.;

Expand All @@ -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.;

Expand All @@ -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;
Expand All @@ -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
{
Expand Down
10 changes: 5 additions & 5 deletions core/io.c
View file
Open in desktop

@kelslund kelslund Apr 1, 2026

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should there be something similar int neutrino_oscillations_ffi = NEUTRINO_OSCILLATIONS_FFI for CFI? Same for diag.c

@Yurlungur Yurlungur Apr 1, 2026

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes. This MR actually isn't done. :) Can you review the one below too? That's the one I was hoping for a review of.
#50

Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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);
Expand All @@ -1169,7 +1169,7 @@ void dump() {
TYPE_DBL);
#undef RANK
}
#endif // NEUTRINO_OSCILLATIONS
#endif // NEUTRINO_OSCILLATIONS_FFI
#endif // LOCAL_ANGULAR_DISTRIBUTIONS

#if RZ_HISTOGRAMS
Expand All @@ -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
Expand Down
4 changes: 4 additions & 0 deletions core/make_superphotons.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.;
Expand Down Expand Up @@ -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;

Expand Down Expand Up @@ -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
Expand All @@ -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) {
Expand Down
Loading