hyschive · hfhsieh · Feb 19, 2023 · Feb 19, 2023 · Feb 19, 2023 · Feb 19, 2023
diff --git a/configs/ncts.config b/configs/ncts.config
@@ -1,9 +1,11 @@
 # NCTS (default: openmpi-intel)
-CUDA_PATH       /cluster/cuda-11.6
+CUDA_PATH       /cluster/cuda-12.4
 FFTW2_PATH      /cluster/intel-2022.1/fftw-2.1.5_openmpi4
 MPI_PATH        /cluster/intel-2022.1/openmpi-4.1.2
 HDF5_PATH       /cluster/intel-2022.1/hdf5-1.10.8
-GSL_PATH        /usr
+GRACKLE_PATH
+GSL_PATH
+LIBYT_PATH
 
 # compilers
 CXX     icpc
@@ -26,6 +28,7 @@ NVCCFLAG_COM -O3
 NVCCFLAG_FLU -Xptxas -dlcm=ca -prec-div=false -ftz=true
 NVCCFLAG_POT -Xptxas -dlcm=ca
 
-#gpu
-GPU_COMPUTE_CAPABILITY 860    # RTX A6000
-#GPU_COMPUTE_CAPABILITY 750    # Quadro RTX 6000
+# gpu
+GPU_COMPUTE_CAPABILITY 750    # Quadro RTX 6000 (gpuserver01 and 02)
+#GPU_COMPUTE_CAPABILITY 860    # RTX A6000 (gpuserver03 and 04)
+#GPU_COMPUTE_CAPABILITY 890    # L40S (gpuserver05 and 06)
diff --git a/example/test_problem/Hydro/CCSN/Input__Flag_Lohner.Leakage b/example/test_problem/Hydro/CCSN/Input__Flag_Lohner.Leakage
@@ -0,0 +1,10 @@
+# Level    Threshold_Refine  Threshold_Derefine              Filter              Soften          MinDensity
+      0                0.80                0.80                0.01                0.00                0.00
+      1                0.80                0.80                0.01                0.00                0.00
+      2                0.80                0.80                0.01                0.00                0.00
+      3                0.80                0.80                0.01                0.00                0.00
+      4                0.80                0.80                0.01                0.00                0.00
+      5                0.80                0.80                0.01                0.00                0.00
+      6                0.80                0.80                0.01                0.00                0.00
+      7                0.80                0.80                0.01                0.00                0.00
+      8                0.80                0.80                0.01                0.00                0.00
diff --git a/example/test_problem/Hydro/CCSN/Input__Flag_User.Leakage b/example/test_problem/Hydro/CCSN/Input__Flag_User.Leakage
@@ -0,0 +1,10 @@
+# Level                Density
+      0                3.0e-9
+      1                1.0e-8
+      2                3.0e-8
+      3                1.0e-7
+      4                3.0e-7
+      5                3.0e-7
+      6                3.0e-7
+      7                3.0e-7
+      8                3.0e-7
diff --git a/example/test_problem/Hydro/CCSN/Input__Parameter.CoreCollapse b/example/test_problem/Hydro/CCSN/Input__Parameter.CoreCollapse
@@ -24,8 +24,8 @@ END_STEP                     -1           # end step (<0=auto -> must be set by
 
 
 # test problems
-TESTPROB_ID                  24           # test problem ID [0]
-                                          #   24: HYDRO CCSN (+GRAVITY & GREP & EOS_NUCLEAR) [+MHD]
+TESTPROB_ID                  50           # test problem ID [0]
+                                          #   50: HYDRO CCSN (+GRAVITY & GREP & EOS_NUCLEAR) [+MHD]
 
 
 # code units (in cgs)

diff --git a/example/test_problem/Hydro/CCSN/Input__Parameter.Leakage b/example/test_problem/Hydro/CCSN/Input__Parameter.Leakage
diff --git a/example/test_problem/Hydro/CCSN/Input__Parameter.Lightbulb b/example/test_problem/Hydro/CCSN/Input__Parameter.Lightbulb
@@ -24,8 +24,8 @@ END_STEP                     -1           # end step (<0=auto -> must be set by
 
 
 # test problems
-TESTPROB_ID                  24           # test problem ID [0]
-                                          #   24: HYDRO CCSN (+GRAVITY & GREP & EOS_NUCLEAR) [+MHD]
+TESTPROB_ID                  50           # test problem ID [0]
+                                          #   50: HYDRO CCSN (+GRAVITY & GREP & EOS_NUCLEAR) [+MHD]
 
 
 # code units (in cgs)

diff --git a/example/test_problem/Hydro/CCSN/Input__Parameter.MigrationTest b/example/test_problem/Hydro/CCSN/Input__Parameter.MigrationTest
@@ -24,8 +24,8 @@ END_STEP                     -1           # end step (<0=auto -> must be set by
 
 
 # test problems
-TESTPROB_ID                  24           # test problem ID [0]
-                                          #   24: HYDRO CCSN (+GRAVITY & GREP & EOS_GAMMA) [+MHD]
+TESTPROB_ID                  50           # test problem ID [0]
+                                          #   50: HYDRO CCSN (+GRAVITY & GREP & EOS_GAMMA) [+MHD]
 
 
 # code units (in cgs)

diff --git a/example/test_problem/Hydro/CCSN/Input__TestProb.CoreCollapse b/example/test_problem/Hydro/CCSN/Input__TestProb.CoreCollapse
@@ -28,7 +28,7 @@ CCSN_CC_MaxRefine_Dens1   1.0e11                   # central density threshold t
 CCSN_CC_MaxRefine_Dens2   1.0e12                   # central density threshold that reduces the maximum refinement level 2  [1.0e12]
 CCSN_CC_CentralDensFac    1.0e13                   # factor that reduces the dt constrained by the central density (in cgs) during the core collapse  [1.0e13]
 CCSN_CC_Red_DT            1.0e-5                   # reduced time step (in s) when the central density exceeds CCSN_CC_CentralDensFac before bounce  [1e-5]
-CCSN_LB_TimeFac           0.10                     # factor that scales the dt constrained by lightbulb scheme  [0.1]
+CCSN_NuHeat_TimeFac       0.10                     # factor that scales the dt constrained by the lightbulb/leakage scheme  [0.1]
 
 CCSN_CC_Rot               2                        # mode for rotational profile (0:off, 1:analytical, 2:table)                         [2]
                                                    # --> analytical formula: Omega(r)=Omega_0*[R_0^2/(r^2+R_0^2)], where r is the spherical radius

diff --git a/example/test_problem/Hydro/CCSN/Input__TestProb.Leakage b/example/test_problem/Hydro/CCSN/Input__TestProb.Leakage
@@ -0,0 +1,51 @@
+# problem-specific runtime parameters
+CCSN_Prob        2                                 # target CCSN problem
+                                                   #    0 : GREP migration test
+                                                   #    1 : Post bounce test
+                                                   #    2 : Core collapse test
+
+CCSN_Prof_File   IC/s20_W2007                      # filename of input profile
+
+CCSN_Mag         1                                 # target magnetic field profile                     [1]
+                                                   # ( 0: Ap = B0* varpi^2 * (1 - rho / rho_max)^np * (P / P_max)
+                                                   #   1: Ap = 0.5 * B0 * ( R0^3 / (r^3 + R0^3) ) * r * sin(theta) )
+CCSN_Mag_B0      1.0e14                            # magnetic field strength (in gauss)                [1.e14]
+CCSN_Mag_np      0.0                               # dependence of magnetic field on density           [0.0]
+CCSN_Mag_R0      1.0e8                             # characteristic radius of magnetic field (in cm)   [1.0e8]
+
+CCSN_GW_OUTPUT   1                                 # output GW signal (0=off, 1=on)                    [0]
+CCSN_GW_DT       0.1                               # output interval of GW signal                      [1.0]
+
+CCSN_Eint_Mode   1                                 # Mode of obtaining internal energy in SetGridIC()  [2]
+                                                   # ( 1=Temp Mode: Eint(dens, temp, [Ye])
+                                                   #   2=Pres Mode: Eint(dens, pres, [Ye]) )
+
+CCSN_CC_MaxRefine_Flag1   1                        # enable limiting maximum refinement level 1  [0]
+CCSN_CC_MaxRefine_Flag2   1                        # enable limiting maximum refinement level 2  [0]
+CCSN_CC_MaxRefine_LV1     6                        # reduced maximum refinement level 1 to this value  [MAX_LEVEL-2]
+CCSN_CC_MaxRefine_LV2     7                        # reduced maximum refinement level 2 to this value  [MAX_LEVEL-1]
+CCSN_CC_MaxRefine_Dens1   1.0e11                   # central density threshold that reduces the maximum refinement level 1  [1.0e11]
+CCSN_CC_MaxRefine_Dens2   1.0e12                   # central density threshold that reduces the maximum refinement level 2  [1.0e12]
+CCSN_CC_CentralDensFac    1.0e13                   # factor that reduces the dt constrained by the central density (in cgs) during the core collapse  [1.0e13]
+CCSN_CC_Red_DT            1.0e-5                   # reduced time step (in s) when the central density exceeds CCSN_CC_CentralDensFac before bounce  [1e-5]
+CCSN_NuHeat_TimeFac       0.01                     # factor that scales the dt constrained by the lightbulb/leakage scheme  [0.1]
+
+CCSN_CC_Rot               0                        # mode for rotational profile (0:off, 1:analytical, 2:table)                         [2]
+                                                   # --> analytical formula: Omega(r)=Omega_0*[R_0^2/(r^2+R_0^2)], where r is the spherical radius
+CCSN_CC_Rot_R0            1.0e8                    # characteristic radius R_0 (in cm) in the analytical rotational profile             [2.0e8]
+CCSN_CC_Rot_Omega0        0.0                      # central angular frequency Omega_0 (in rad/s) in the analytical rotational profile  [0.5]
+CCSN_CC_Rot_Fac          -1.0                      # multiplication factor for the tabular rotational profile (<=0 -> off)              [-1.0]
+
+CCSN_REF_RBase           -1.0;                     # reference distance for determining a maximum refinement level based on distance from the box center (in cm) [1.25e7]
+
+CCSN_Is_PostBounce        0                        # boolean that indicates whether core bounce has occurred  [0]
+
+CCSN_DT_YE                1                        # dt criterion on Ye (1=Ye-Ye_min/Ye_max, 2=Ye, 3=none) [1]
+
+CCSN_MaxRefine_Rad        3.0e6                    # radius in cm within which to refine to the maximum allowed level [3.0e6]
+CCSN_AngRes_Min          -1.0                      # minimum angular resolution in degrees (<=0=off) [-1.0]
+CCSN_AngRes_Max          -1.0                      # maximum angular resolution in degrees (<=0=off) [-1.0]
+
+CCSN_Shock_ThresFac_Pres  0.5                      # pressure threshold factor for detecting postbounce shock  [0.5]
+CCSN_Shock_ThresFac_Vel   0.1                      # velocity threshold facotr for detecting postbounce shock  [0.1]
+CCSN_Shock_Weight         2                        # weighting method of the averaged shock radius (1:volume, 2:1/volume) [2]
diff --git a/example/test_problem/Hydro/CCSN/Input__TestProb.Lightbulb b/example/test_problem/Hydro/CCSN/Input__TestProb.Lightbulb
@@ -20,7 +20,7 @@ CCSN_Eint_Mode   1                                 # Mode of obtaining internal
                                                    # ( 1=Temp Mode: Eint(dens, temp, [Ye])
                                                    #   2=Pres Mode: Eint(dens, pres, [Ye]) )
 
-CCSN_LB_TimeFac           0.10                     # factor that scales the dt constrained by lightbulb scheme  [0.1]
+CCSN_NuHeat_TimeFac       0.10                     # factor that scales the dt constrained by the lightbulb/leakage scheme  [0.1]
 
 CCSN_Is_PostBounce        0                        # boolean that indicates whether core bounce has occurred    [0]
 

diff --git a/example/test_problem/Hydro/CCSN/README b/example/test_problem/Hydro/CCSN/README
@@ -30,11 +30,17 @@ Default setup:
 4. For Lightbulb:
    --> OPT__DT_USER = 1
    --> OPT__RESET_FLUID = 1
+   --> NEUTRINO_SCHEME = LIGHTBULB
 
 5. For Core collapse:
    --> MAX_LEVEL = 9
    --> OPT__RESET_FLUID = 1
 
+6. For Leakage:
+   --> OPT__DT_USER = 1
+   --> OPT__RESET_FLUID = 1
+   --> NEUTRINO_SCHEME = LEAKAGE
+
 Note:
 ========================================
 1. A suite of test problems for core-collapse supernova simulations

diff --git a/example/test_problem/Hydro/CCSN/clean.sh b/example/test_problem/Hydro/CCSN/clean.sh
@@ -3,4 +3,5 @@ rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Du
       Diag* Box* BaseXYslice* BaseYZslice* BaseXZslice* BaseXline* BaseYline* BaseZline* BaseDiag* \
       PowerSpec_* Particle_* nohup.out Record__Performance Record__TimingMPI_* \
       Record__ParticleCount Record__User Patch_* Record__NCorrUnphy FailedPatchGroup* *.pyc Record__LoadBalance \
-      GRACKLE_INFO Record__DivB Record__CentralQuant Record__QuadMom_2nd GREP_Lv*
+      GRACKLE_INFO Record__DivB Record__Center \
+      Record__CentralQuant Record__QuadMom_2nd Leakage_Lum_* GREP_Lv*
diff --git a/example/test_problem/Hydro/CCSN/generate_make.sh b/example/test_problem/Hydro/CCSN/generate_make.sh
@@ -2,6 +2,6 @@
 
 PYTHON=python3
 
-${PYTHON} configure.py --mpi=true --hdf5=true --fftw=FFTW3 --gpu=true --debug=false \
+${PYTHON} configure.py --mpi=true --hdf5=true --fftw=FFTW3 --gpu=true --gpu_regcount_flu=255 --debug=false \
                        --model=HYDRO --flu_scheme=MHM_RP --slope=PPM --mhd=false --gravity=true --unsplit_gravity=false \
                        --eos=NUCLEAR --nuc_table=TEMP --nuc_solver=ORIG --neutrino=LIGHTBULB --grep=true "$@"
diff --git a/example/test_problem/Template/Input__Parameter b/example/test_problem/Template/Input__Parameter
@@ -48,7 +48,7 @@ TESTPROB_ID                   0           # test problem ID [0]
                                           #   20: HYDRO MHD Cosmic Ray Soundwave
                                           #   21: HYDRO MHD Cosmic Ray Shocktube
                                           #   23: HYDRO MHD Cosmic Ray Diffusion
-                                          #   24: HYDRO CCSN (+GRAVITY & GREP & EOS_GAMMA/EOS_NUCLEAR) [+MHD]
+                                          #   50: HYDRO CCSN (+GRAVITY & GREP & EOS_GAMMA/EOS_NUCLEAR) [+MHD]
                                           #  100: HYDRO CDM cosmological simulation (+GRAVITY & COMOVING & PARTICLE)
                                           #  101: HYDRO Zeldovich pancake collapse (+GRAVITY & COMOVING & PARTICLE)
                                           # 1000: ELBDM external potential (+GRAVITY)
@@ -208,6 +208,18 @@ SRC_DELEP_YEC                 0.0601082693    # deleptonization parameter Ye_c [
 SRC_LIGHTBULB                 0           # lightbulb scheme for neutrino heating and cooling [0] ##HYDRO ONLY##
 SRC_LIGHTBULB_LNUE            4.0e52      # electron neutrino luminosity in erg/s [1e52] ##HYDRO ONLY##
 SRC_LIGHTBULB_TNUE            4.0         # electron neutrino temperature in MeV [4.0] ##HYDRO ONLY##
+SRC_LEAKAGE                   0           # multiflavour neutrino leakage scheme [0]
+SRC_LEAKAGE_NRADIUS           768         # number of bin in the radial    direction [768]
+SRC_LEAKAGE_NTHETA            1           # number of ray in the polar     direction [1]
+SRC_LEAKAGE_NPHI              1           # number of ray in the azimuthal direction [1]
+                                          # (must be 1 or a multiple of 2)
+SRC_LEAKAGE_BINSIZE_RADIUS    1.0e5       # bin size of linear bins in radius in cm [1.0e5]
+SRC_LEAKAGE_RADIUSMAX         3.0e8       # maximum radius in cm [3.0e8]
+SRC_LEAKAGE_RADIUSMIN_LOG     3.0e7       # minimum radius of logarithmic bins in cm [3.0e7]
+SRC_LEAKAGE_NUHEAT            1           # include a parameterized neutrino heating scheme [1]
+SRC_LEAKAGE_NUHEAT_FAC        1.0         # rescale the parameterized heating rate by SRC_LEAKAGE_NUHEAT_FAC [1.0]
+SRC_LEAKAGE_OPT_TEMP          2           # temperature computation scheme [2]
+                                          # (1=individual cell, 2=binned data)
 SRC_USER                      0           # user-defined source terms -> edit "Src_User.cpp" [0]
 SRC_GPU_NPGROUP              -1           # number of patch groups sent into the CPU/GPU source-term solver (<=0=auto) [-1]
 
@@ -261,6 +273,8 @@ MOLECULAR_WEIGHT              0.6         # mean molecular weight [0.6]
 MU_NORM                      -1.0         # normalization of MOLECULAR_WEIGHT (<0=m_H, 0=amu, >0=input manually) [-1.0]
 ISO_TEMP                      1.0e4       # isothermal temperature in kelvin ##EOS_ISOTHERMAL ONLY##
 NUC_TABLE                     LS220_eps.h5   # nuclear EoS table ##EOS_NUCLEAR ONLY##
+NUC_INT_SCHEME_MAIN           1           # interpolation scheme for nuclear EoS table (1=linear, 2=cubic) [1] ##EOS_NUCLEAR ONLY##
+NUC_INT_SCHEME_AUX            1           # interpolation scheme for auxiliary   table (1=linear, 2=cubic) [1] ##EOS_NUCLEAR ONLY##
 MINMOD_COEFF                  1.5         # coefficient of the generalized MinMod limiter (1.0~2.0) [1.5]
 MINMOD_MAX_ITER               0           # maximum number of iterations to reduce MINMOD_COEFF when data reconstruction fails (0=off) [0]
 OPT__LR_LIMITER              -1           # slope limiter of data reconstruction in the MHM/MHM_RP/CTU schemes:

diff --git a/include/CUDA_ConstMemory.h b/include/CUDA_ConstMemory.h
@@ -36,6 +36,8 @@ SET_GLOBAL( __constant__ double c_Src_Dlep_AuxArray_Flt     [SRC_NAUX_DLEP     ]
 SET_GLOBAL( __constant__ int    c_Src_Dlep_AuxArray_Int     [SRC_NAUX_DLEP     ] );
 SET_GLOBAL( __constant__ double c_Src_Lightbulb_AuxArray_Flt[SRC_NAUX_LIGHTBULB] );
 SET_GLOBAL( __constant__ int    c_Src_Lightbulb_AuxArray_Int[SRC_NAUX_LIGHTBULB] );
+SET_GLOBAL( __constant__ double c_Src_Leakage_AuxArray_Flt  [SRC_NAUX_LEAKAGE  ] );
+SET_GLOBAL( __constant__ int    c_Src_Leakage_AuxArray_Int  [SRC_NAUX_LEAKAGE  ] );
 #endif
 SET_GLOBAL( __constant__ double c_Src_User_AuxArray_Flt     [SRC_NAUX_USER     ] );
 SET_GLOBAL( __constant__ int    c_Src_User_AuxArray_Int     [SRC_NAUX_USER     ] );

diff --git a/include/Field.h b/include/Field.h
@@ -25,6 +25,9 @@ SET_GLOBAL( FieldIdx_t Idx_MomZ,          Idx_Undefined );
 SET_GLOBAL( FieldIdx_t Idx_Engy,          Idx_Undefined );
 #if ( EOS == EOS_NUCLEAR )
 SET_GLOBAL( FieldIdx_t Idx_Ye,            Idx_Undefined );
+#if ( NEUTRINO_SCHEME == LEAKAGE )
+SET_GLOBAL( FieldIdx_t Idx_dYedt_Nu,      Idx_Undefined );
+#endif
 SET_GLOBAL( FieldIdx_t Idx_dEdt_Nu,       Idx_Undefined );
 #if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP )
 SET_GLOBAL( FieldIdx_t Idx_Temp_IG,       Idx_Undefined );

diff --git a/include/Global.h b/include/Global.h
@@ -356,6 +356,8 @@ extern double     Src_Dlep_AuxArray_Flt     [SRC_NAUX_DLEP     ];
 extern int        Src_Dlep_AuxArray_Int     [SRC_NAUX_DLEP     ];
 extern double     Src_Lightbulb_AuxArray_Flt[SRC_NAUX_LIGHTBULB];
 extern int        Src_Lightbulb_AuxArray_Int[SRC_NAUX_LIGHTBULB];
+extern double     Src_Leakage_AuxArray_Flt  [SRC_NAUX_LEAKAGE  ];
+extern int        Src_Leakage_AuxArray_Int  [SRC_NAUX_LEAKAGE  ];
 #endif
 extern double     Src_User_AuxArray_Flt     [SRC_NAUX_USER     ];
 extern int        Src_User_AuxArray_Int     [SRC_NAUX_USER     ];
@@ -491,6 +493,15 @@ extern real       (*h_Mag_Array_S_In [2])[NCOMP_MAG ][ SRC_NXT_P1*SQR(SRC_NXT) ]
 #endif
 extern double     (*h_Corner_Array_S[2])[3];
 
+#if ( MODEL == HYDRO )
+extern real        *h_SrcLeakage_Radius;
+extern real        *h_SrcLeakage_tau;
+extern real        *h_SrcLeakage_chi;
+extern real        *h_SrcLeakage_HeatFlux;
+extern real        *h_SrcLeakage_HeatERms;
+extern real        *h_SrcLeakage_HeatEAve;
+#endif
+
 
 
 // 4/5. GPU (device) global memory arrays and timers