Skip to content
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.

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

local angular distributions #36

Merged
merged 55 commits into from
Dec 11, 2024
Merged
Show file tree
Hide file tree
Changes from 51 commits
Commits
Show all changes
55 commits
Select commit Hold shift + click to select a range
6c5351f
add ability to generate histograms
jonahm-LANL Jul 10, 2024
1e4ad19
atomic!
jonahm-LANL Jul 10, 2024
bafaacc
oscillations in situ machinery seems to run
jonahm-LANL Nov 17, 2024
021287c
loop over photons broken with tracers. Now fixed.
jonahm-LANL Nov 18, 2024
14d9b3c
add local angular distributions to torus_cbc
jonahm-LANL Nov 18, 2024
5538e6e
darwin
jonahm-LANL Nov 18, 2024
e951146
fundamental flaw... coords only specified for Minkowski. now fixed.
jonahm-LANL Nov 18, 2024
437fbe3
Merge branch 'jmm/local-angular-distributions' of github.com:lanl/nub…
jonahm-LANL Nov 18, 2024
deaee4c
testing
jonahm-LANL Nov 19, 2024
2902ca7
whee
jonahm-LANL Nov 19, 2024
558dd53
infinity...
jonahm-LANL Nov 19, 2024
7445dfb
why is G = 0
jonahm-LANL Nov 19, 2024
9a8f349
pass in struct pointer???
jonahm-LANL Nov 19, 2024
722b6aa
fabs
jonahm-LANL Nov 19, 2024
8f97909
progress
jonahm-LANL Nov 19, 2024
993a836
try normalizing kspace
jonahm-LANL Nov 19, 2024
2b22cde
cleanup...
jonahm-LANL Nov 19, 2024
426e855
one more try
jonahm-LANL Nov 19, 2024
31d8c80
i > mu
jonahm-LANL Nov 19, 2024
534e4d5
should be gcov not gcon
jonahm-LANL Nov 19, 2024
8269258
cleanup
jonahm-LANL Nov 19, 2024
cdfff36
Support 4 species transport
jonahm-LANL Nov 21, 2024
4ec9374
fornax 1zone test modified for 4 species NuLib table
jonahm-LANL Nov 21, 2024
1bf5e6e
output local angles grid positions and mu
jonahm-LANL Nov 21, 2024
4da3b69
Define Gnu, moment types but don't set them yet. Refactor binning cod…
jonahm-LANL Nov 21, 2024
e9592fd
compute Gnu and moments
jonahm-LANL Nov 22, 2024
52b3dfe
Better include guards
jonahm-LANL Nov 22, 2024
90a3010
Cleanup...
jonahm-LANL Nov 22, 2024
cebaaed
dt and oscillations
jonahm-LANL Nov 23, 2024
ada6bdc
OSCILLATIONS
jonahm-LANL Nov 23, 2024
37c4ae9
oscillations doesn't depend on time. Make sure dt is in code units.
jonahm-LANL Nov 23, 2024
7b6f7e5
thread it through step.c
jonahm-LANL Nov 23, 2024
547ea4b
better formatting
jonahm-LANL Nov 23, 2024
fa23e77
multiscatt getting crazier...
jonahm-LANL Nov 23, 2024
5f824d9
starting oscillations setup. not done yet.
jonahm-LANL Nov 23, 2024
fa0112e
oops didn't mean to commit that
jonahm-LANL Nov 23, 2024
039644e
NuLib opacities
jonahm-LANL Nov 23, 2024
4bcc147
debugging, performance, cleanup
jonahm-LANL Nov 23, 2024
420f5f8
oscillations problem test works
jonahm-LANL Nov 23, 2024
903c927
simplify/cleanup
jonahm-LANL Nov 23, 2024
b55f1f8
automated test
jonahm-LANL Nov 23, 2024
be16433
revert intermediate fornax and multiscatt tests
jonahm-LANL Nov 23, 2024
f4631f1
can do it badly
jonahm-LANL Nov 25, 2024
c09f8a1
fix
jonahm-LANL Nov 25, 2024
cbb3f38
NOW the tests look right
jonahm-LANL Nov 28, 2024
7da4cb8
update to use Payels formula for std devation, with correct units
jonahm-LANL Dec 8, 2024
7756675
fix formula
jonahm-LANL Dec 8, 2024
1e2eeb2
add timers
jonahm-LANL Dec 8, 2024
0decec3
local moments
jonahm-LANL Dec 8, 2024
7f6ca37
add local angles positions to geom reader
jonahm-LANL Dec 8, 2024
7a6e870
update
jonahm-LANL Dec 9, 2024
e6c3104
kalunds comments
jonahm-LANL Dec 11, 2024
2e2bb2a
revert energy range for torus_cbc
jonahm-LANL Dec 11, 2024
3c84807
oscillations flag
jonahm-LANL Dec 11, 2024
84cc9c0
oscillations flag for opacity too
jonahm-LANL Dec 11, 2024
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
4 changes: 4 additions & 0 deletions README.md
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 4 additions & 0 deletions core/constants.h
Original file line number Diff line number Diff line change
Expand Up @@ -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)
kelslund marked this conversation as resolved.
Show resolved Hide resolved
#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
Expand Down
12 changes: 11 additions & 1 deletion core/coord.c
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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));
Expand Down
56 changes: 50 additions & 6 deletions core/decs.h
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -294,9 +292,10 @@
#define TIMER_MAKE (7)
#define TIMER_PUSH (8)
#define TIMER_INTERACT (9)
#define TIMER_MICRO (10)
#define TIMER_ALL (11)
#define NUM_TIMERS (12)
#define TIMER_OSCILLATIONS (10)
#define TIMER_MICRO (11)
#define TIMER_ALL (12)
#define NUM_TIMERS (13)

// Units
#define NEED_UNITS (RADIATION || EOS == EOS_TYPE_TABLE)
Expand Down Expand Up @@ -368,6 +367,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)
Comment on lines +372 to +374
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

might want to move these up...

#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]
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
[LOCAL_ANGLES_NX2];
typedef int grid_local_basis_idx_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_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
Expand Down Expand Up @@ -671,6 +694,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.)
Expand Down Expand Up @@ -981,6 +1010,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_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);
#endif // RAD_NUM_TYPES >= 4
#endif // LOCAL_ANGULAR_DISTRIBUTIONS
#endif // RADIATION

// passive.c
Expand Down Expand Up @@ -1065,6 +1108,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);
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

needed for generalizing to 4 species transport

#endif // NEUTRINOS
#endif // RADIATION

Expand Down
10 changes: 10 additions & 0 deletions core/defs.h
Original file line number Diff line number Diff line change
Expand Up @@ -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_Ns, local_wsqr;
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
Expand Down
9 changes: 8 additions & 1 deletion core/diag.c
Original file line number Diff line number Diff line change
Expand Up @@ -298,6 +298,9 @@ void diag(int call_code) {
get_time_per_step(TIMER_ALL));
#if ELECTRONS
fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_ELECTRON));
#endif
#if NEUTRINO_OSCILLATIONS || LOCAL_ANGULAR_DISTRIBUTIONS
Copy link
Collaborator

Choose a reason for hiding this comment

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

What is the difference between these two conditions?

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

One can compute G, A, B, etc without actually running the oscillations code. So NEUTRINO_OSCILLATIONS requires LOCAL_ANGULAR_DISTRIBUTIONS but not vice versa. Would be sufficient to include only the latter, but I was being careful.

fprintf(ener_file, "%15.8g ", get_time_per_step(TIMER_OSCILLATIONS));
#endif
fprintf(ener_file, "\n");
fflush(ener_file);
Expand Down Expand Up @@ -583,8 +586,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
Yurlungur marked this conversation as resolved.
Show resolved Hide resolved
fprintf(stdout, " ANTIX %.2f %%\n",
rad_type_counts[ANTINU_HEAVY]);
#endif
fprintf(stdout, "*********************************\n\n");
}
timer_stop(TIMER_DIAG);
Expand Down
Loading
Loading