diff --git a/configs/ncts.config b/configs/ncts.config index e3f156876..958d3c90d 100644 --- 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 new file mode 100644 index 000000000..b77bea38e --- /dev/null +++ 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 new file mode 100644 index 000000000..ef047328f --- /dev/null +++ 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 index 819bff7fd..d60fad79b 100644 --- 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 new file mode 100644 index 000000000..5223e06fa --- /dev/null +++ b/example/test_problem/Hydro/CCSN/Input__Parameter.Leakage @@ -0,0 +1,351 @@ + + +# ================================================================================================================= +# NOTE: +# 1. Comment symbol: # +# 2. [*]: defaults +# 3. Parameters set to "auto" (usually by setting to a negative value) do not have deterministic default values +# and will be set according to the adopted compilation options and/or other runtime parameters +# 4. To add new parameters, please edit "Init/Init_Load_Parameter.cpp" +# 5. All dimensional variables should be set consistently with the code units (set by UNIT_L/M/T/V/D) unless +# otherwise specified (e.g., SF_CREATE_STAR_MIN_GAS_DENS & SF_CREATE_STAR_MIN_STAR_MASS) +# 6. For boolean options: 0/1 -> off/on +# ================================================================================================================= + + +# simulation scale +BOX_SIZE 20000.0 # box size along the longest side (in Mpc/h if COMOVING is adopted) +NX0_TOT_X 160 # number of base-level cells along x +NX0_TOT_Y 160 # number of base-level cells along y +NX0_TOT_Z 160 # number of base-level cells along z +OMP_NTHREAD 32 # number of OpenMP threads (<=0=auto) [-1] ##OPENMP ONLY## +END_T 2000.0 # end physical time (<0=auto -> must be set by test problems or restart) [-1.0] +END_STEP -1 # end step (<0=auto -> must be set by test problems or restart) [-1] + + +# test problems +TESTPROB_ID 50 # test problem ID [0] + # 50: HYDRO CCSN (+GRAVITY & GREP & EOS_GAMMA/EOS_NUCLEAR) [+MHD] + + +# code units (in cgs) +OPT__UNIT 1 # specify code units -> must set exactly 3 basic units below [0] ##USELESS FOR COMOVING## +UNIT_L 1.0e5 # length unit (<=0 -> set to UNIT_V*UNIT_T or (UNIT_M/UNIT_D)^(1/3)) [-1.0] +UNIT_M 1.989e33 # mass unit (<=0 -> set to UNIT_D*UNIT_L^3) [-1.0] +UNIT_T 1.0e-3 # time unit (<=0 -> set to UNIT_L/UNIT_V) [-1.0] +UNIT_V -1.0 # velocity unit (<=0 -> set to UNIT_L/UNIT_T) [-1.0] +UNIT_D -1.0 # mass density unit (<=0 -> set to UNIT_M/UNIT_L^3) [-1.0] + + +# boundary conditions +OPT__BC_FLU_XM 2 # fluid boundary condition at the -x face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_XP 2 # fluid boundary condition at the +x face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_YM 2 # fluid boundary condition at the -y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_YP 2 # fluid boundary condition at the +y face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_ZM 2 # fluid boundary condition at the -z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_FLU_ZP 2 # fluid boundary condition at the +z face: (1=periodic, 2=outflow, 3=reflecting, 4=user, 5=diode) ##2/3/5 for HYDRO ONLY## +OPT__BC_POT 2 # gravity boundary condition: (1=periodic, 2=isolated) +GFUNC_COEFF0 -1.0 # Green's function coefficient at the origin for the isolated BC (<0=auto) [-1.0] + + +# time-step +DT__MAX 0.1 # dt criterion: maximum allowed dt (<0=off) [-1.0] +DT__FLUID -1.0 # dt criterion: fluid solver CFL factor (<0=auto) [-1.0] +DT__FLUID_INIT 1.0e-5 # dt criterion: DT__FLUID at the first step (<0=auto) [-1.0] +DT__SPEED_OF_LIGHT 0 # dt criterion: speed of light [0] ##SRHD ONLY## +DT__GRAVITY -1.0 # dt criterion: gravity solver safety factor (<0=auto) [-1.0] +DT__PHASE 0.0 # dt criterion: phase rotation safety factor (0=off) [0.0] ##ELBDM ONLY## +DT__PARVEL 0.5 # dt criterion: particle velocity safety factor [0.5] +DT__PARVEL_MAX -1.0 # dt criterion: maximum allowed dt from particle velocity (<0=off) [-1.0] +DT__PARACC 0.5 # dt criterion: particle acceleration safety factor (0=off) [0.5] ##STORE_PAR_ACC ONLY## +DT__MAX_DELTA_A 0.01 # dt criterion: maximum variation of the cosmic scale factor [0.01] +DT__SYNC_PARENT_LV 0.1 # dt criterion: allow dt to adjust by (1.0+DT__SYNC_PARENT) in order to synchronize + # with the parent level (for OPT__DT_LEVEL==3 only) [0.1] +DT__SYNC_CHILDREN_LV 0.1 # dt criterion: allow dt to adjust by (1.0-DT__SYNC_CHILDREN) in order to synchronize + # with the children level (for OPT__DT_LEVEL==3 only; 0=off) [0.1] +OPT__DT_USER 1 # dt criterion: user-defined -> edit "Mis_GetTimeStep_UserCriteria.cpp" [0] +OPT__DT_LEVEL 3 # dt at different AMR levels (1=shared, 2=differ by two, 3=flexible) [3] +OPT__RECORD_DT 1 # record info of the dt determination [1] +AUTO_REDUCE_DT 1 # reduce dt automatically when the program fails (for OPT__DT_LEVEL==3 only) [1] +AUTO_REDUCE_DT_FACTOR 0.8 # reduce dt by a factor of AUTO_REDUCE_DT_FACTOR when the program fails [1.0] +AUTO_REDUCE_DT_FACTOR_MIN 0.005 # minimum allowed AUTO_REDUCE_DT_FACTOR after consecutive failures [0.1] +AUTO_REDUCE_MINMOD_FACTOR 0.8 # reduce MINMOD_COEFF by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] ##HYDRO ONLY## +AUTO_REDUCE_MINMOD_MIN 1.0e-2 # minimum allowed MINMOD_COEFF after consecutive failures [1.0e-2] ##HYDRO ONLY## +AUTO_REDUCE_INT_MONO_FACTOR 0.8 # reduce INT_MONO_COEFF(_B) by this factor together with AUTO_REDUCE_DT (1.0=off) [0.8] +AUTO_REDUCE_INT_MONO_MIN 1.0e-2 # minimum allowed INT_MONO_COEFF(_B) after consecutive failures [1.0e-2] + + +# grid refinement (examples of Input__Flag_XXX tables are put at "example/input/") +REGRID_COUNT 4 # refine every REGRID_COUNT sub-step [4] +REFINE_NLEVEL 1 # number of new AMR levels to be created at once during refinement [1] +FLAG_BUFFER_SIZE 8 # number of buffer cells for the flag operation (0~PATCH_SIZE; <0=auto -> PATCH_SIZE) [-1] +FLAG_BUFFER_SIZE_MAXM1_LV 4 # FLAG_BUFFER_SIZE at the level MAX_LEVEL-1 (<0=auto -> REGRID_COUNT) [-1] +FLAG_BUFFER_SIZE_MAXM2_LV -1 # FLAG_BUFFER_SIZE at the level MAX_LEVEL-2 (<0=auto) [-1] +MAX_LEVEL 8 # maximum refinement level (0~NLEVEL-1) [NLEVEL-1] +OPT__FLAG_RHO 0 # flag: density (Input__Flag_Rho) [0] +OPT__FLAG_RHO_GRADIENT 0 # flag: density gradient (Input__Flag_RhoGradient) [0] +OPT__FLAG_PRES_GRADIENT 0 # flag: pressure gradient (Input__Flag_PresGradient) [0] ##HYDRO ONLY## +OPT__FLAG_VORTICITY 0 # flag: vorticity (Input__Flag_Vorticity) [0] ##HYDRO ONLY## +OPT__FLAG_JEANS 0 # flag: Jeans length (Input__Flag_Jeans) [0] ##HYDRO ONLY## +OPT__FLAG_CURRENT 0 # flag: current density in MHD (Input__Flag_Current) [0] ##MHD ONLY## +OPT__FLAG_ENGY_DENSITY 0 # flag: energy density (Input_Flag_EngyDensity) [0] ##ELBDM ONLY## +OPT__FLAG_LOHNER_DENS 1 # flag: Lohner for mass density (Input__Flag_Lohner) [0] ##BOTH HYDRO AND ELBDM## +OPT__FLAG_LOHNER_ENGY 0 # flag: Lohner for energy density (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_PRES 1 # flag: Lohner for pressure (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_TEMP 0 # flag: Lohner for temperature (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_ENTR 1 # flag: Lohner for entropy (Input__Flag_Lohner) [0] ##HYDRO ONLY## +OPT__FLAG_LOHNER_CRAY 0 # flag: Lohner for cosmic-ray energy (Input__Flag_Lohner) [0] ##COSMIC_RAY ONLY## +OPT__FLAG_LOHNER_FORM 2 # form of Lohner: (1=FLASH-1, 2=FLASH-2, 3=form-invariant-1, 4=form-invariant-2) [2] +OPT__FLAG_USER 1 # flag: user-defined (Input__Flag_User) -> edit "Flag_User.cpp" [0] +OPT__FLAG_USER_NUM 1 # number of threshold values in user-defined table (Input__Flag_User) [1] +OPT__FLAG_REGION 1 # flag: specify the regions **allowed** to be refined -> edit "Flag_Region.cpp" [0] +OPT__FLAG_NPAR_PATCH 0 # flag: # of particles per patch (Input__Flag_NParPatch): (0=off, 1=itself, 2=itself+siblings) [0] +OPT__FLAG_NPAR_CELL 0 # flag: # of particles per cell (Input__Flag_NParCell) [0] +OPT__FLAG_PAR_MASS_CELL 0 # flag: total particle mass per cell (Input__Flag_ParMassCell) [0] +OPT__NO_FLAG_NEAR_BOUNDARY 0 # flag: disallow refinement near the boundaries [0] +OPT__PATCH_COUNT 1 # record the # of patches at each level: (0=off, 1=every step, 2=every sub-step) [1] +OPT__PARTICLE_COUNT 1 # record the # of particles at each level: (0=off, 1=every step, 2=every sub-step) [1] +OPT__REUSE_MEMORY 2 # reuse patch memory to reduce memory fragmentation: (0=off, 1=on, 2=aggressive) [2] +OPT__MEMORY_POOL 0 # preallocate patches for OPT__REUSE_MEMORY=1/2 (Input__MemoryPool) [0] + + +# load balance (LOAD_BALANCE only) +LB_INPUT__WLI_MAX 0.1 # weighted-load-imbalance (WLI) threshold for redistributing all patches [0.1] +LB_INPUT__PAR_WEIGHT 0.0 # load-balance weighting of one particle over one cell [0.0] +OPT__RECORD_LOAD_BALANCE 1 # record the load-balance info [1] +OPT__MINIMIZE_MPI_BARRIER 0 # minimize MPI barriers to improve load balance, especially with particles [0] + # (STORE_POT_GHOST, PAR_IMPROVE_ACC=1, OPT__TIMING_BARRIER=0 only; recommend AUTO_REDUCE_DT=0) + + +# source terms +SRC_DELEPTONIZATION 1 # deleptonization (for simulations of stellar core collapse) [0] ##HYDRO ONLY## +SRC_DELEP_ENU 10.0 # deleptonization parameter E_nu in MeV [10.0] ##HYDRO ONLY## +SRC_DELEP_RHO1 4.9316595742e7 # deleptonization parameter rho_1 in g/cm^3 [3.0e7 ] ##HYDRO ONLY## +SRC_DELEP_RHO2 1.5882185396e13 # deleptonization parameter rho_2 in g/cm^3 [2.0e13] ##HYDRO ONLY## +SRC_DELEP_YE1 0.5 # deleptonization parameter Ye_1 [0.5 ] ##HYDRO ONLY## +SRC_DELEP_YE2 0.265508866 # deleptonization parameter Ye_2 [0.278] ##HYDRO ONLY## +SRC_DELEP_YEC 0.0601082693 # deleptonization parameter Ye_c [0.035] ##HYDRO ONLY## +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] + + +# fluid solver in HYDRO (MODEL==HYDRO only) +GAMMA 1.666666667 # ratio of specific heats (i.e., adiabatic index) [5.0/3.0] ##EOS_GAMMA ONLY## +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: + # (-1=auto, 0=none, 1=vanLeer, 2=generalized MinMod, 3=vanAlbada, 4=vanLeer+generalized MinMod, 6=central) [-1] +OPT__1ST_FLUX_CORR -1 # correct unphysical results (defined by MIN_DENS/PRES) by the 1st-order fluxes: + # (<0=auto, 0=off, 1=3D, 2=3D+1D) [-1] ##MHM/MHM_RP/CTU ONLY## +OPT__1ST_FLUX_CORR_SCHEME -1 # Riemann solver for OPT__1ST_FLUX_CORR (-1=auto, 0=none, 1=Roe, 2=HLLC, 3=HLLE, 4=HLLD) [-1] +DUAL_ENERGY_SWITCH 2.0e-2 # apply dual-energy if E_int/E_kin < DUAL_ENERGY_SWITCH [2.0e-2] ##DUAL_ENERGY ONLY## +OPT__SAME_INTERFACE_B 0 # ensure B field consistency on the shared interfaces between sibling patches (for debugging) [0] ##MHD ONLY## + + +# fluid solvers in all models +FLU_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU fluid solver (<=0=auto) [-1] +GPU_NSTREAM -1 # number of CUDA streams for the asynchronous memory copy in GPU (<=0=auto) [-1] +OPT__FIXUP_FLUX 1 # correct coarse grids by the fine-grid boundary fluxes [1] ##HYDRO and ELBDM ONLY## +OPT__FIXUP_ELECTRIC 1 # correct coarse grids by the fine-grid boundary electric field [1] ##MHD ONLY## +OPT__FIXUP_RESTRICT 1 # correct coarse grids by averaging the fine-grid data [1] +OPT__CORR_AFTER_ALL_SYNC -1 # apply various corrections after all levels are synchronized (see "Flu_CorrAfterAllSync"): + # (-1=auto, 0=off, 1=every step, 2=before dump) [-1] +OPT__NORMALIZE_PASSIVE 1 # ensure "sum(passive_scalar_density) == gas_density" [1] +OPT__INT_FRAC_PASSIVE_LR 1 # convert specified passive scalars to mass fraction during data reconstruction [1] +OPT__OVERLAP_MPI 0 # overlap MPI communication with CPU/GPU computations [0] ##NOT SUPPORTED YET## +OPT__RESET_FLUID 1 # reset fluid variables after each update -> edit "Flu_ResetByUser.cpp" [0] +OPT__RESET_FLUID_INIT -1 # reset fluid variables during initialization (<0=auto -> OPT__RESET_FLUID, 0=off, 1=on) [-1] +OPT__FREEZE_FLUID 0 # do not evolve fluid at all [0] +OPT__CHECK_PRES_AFTER_FLU -1 # check unphysical pressure at the end of the fluid solver (<0=auto) [-1] +OPT__LAST_RESORT_FLOOR 1 # apply floor values as the last resort when the fluid solver fails [1] ##HYDRO and MHD ONLY## +MIN_DENS 0.0 # minimum mass density (must >= 0.0) [0.0] ##HYDRO, MHD, and ELBDM ONLY## +MIN_PRES 0.0 # minimum pressure (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_EINT 0.0 # minimum internal energy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_TEMP 0.0 # minimum temperature in K (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +MIN_ENTR 0.0 # minimum entropy (must >= 0.0) [0.0] ##HYDRO and MHD ONLY## +JEANS_MIN_PRES 0 # minimum pressure estimated from the Jeans length [0] ##HYDRO/MHD and GRAVITY ONLY## +JEANS_MIN_PRES_LEVEL -1 # for JEANS_MIN_PRES; ensure Jeans length is resolved by JEANS_MIN_PRES_NCELL*dh[JEANS_MIN_PRES_LEVEL] + # (<0=auto -> MAX_LEVEL) [-1] +JEANS_MIN_PRES_NCELL 4 # for JEANS_MIN_PRES; see JEANS_MIN_PRES_LEVEL [4] + + +# gravity solvers in all models +NEWTON_G 1.0 # gravitational constant (will be overwritten if OPT__UNIT or COMOVING is on) +SOR_OMEGA -1.0 # over-relaxation parameter in SOR: (<0=auto) [-1.0] +SOR_MAX_ITER -1 # maximum number of iterations in SOR: (<0=auto) [-1] +SOR_MIN_ITER -1 # minimum number of iterations in SOR: (<0=auto) [-1] +MG_MAX_ITER -1 # maximum number of iterations in multigrid: (<0=auto) [-1] +MG_NPRE_SMOOTH -1 # number of pre-smoothing steps in multigrid: (<0=auto) [-1] +MG_NPOST_SMOOTH -1 # number of post-smoothing steps in multigrid: (<0=auto) [-1] +MG_TOLERATED_ERROR -1.0 # maximum tolerated error in multigrid (<0=auto) [-1.0] +POT_GPU_NPGROUP -1 # number of patch groups sent into the CPU/GPU Poisson solver (<=0=auto) [-1] +OPT__GRA_P5_GRADIENT 0 # 5-points gradient in the Gravity solver (must have GRA/USG_GHOST_SIZE_G>=2) [0] +OPT__GRAVITY_EXTRA_MASS 0 # add extra mass source when computing gravity [0] +OPT__SELF_GRAVITY 1 # add self-gravity [1] +OPT__EXT_ACC 0 # add external acceleration (0=off, 1=function, 2=table) [0] ##HYDRO ONLY## + # --> 2 (table) is not supported yet +OPT__EXT_POT 3 # add external potential (0=off, 1=function, 2=table, 3=GREP) [0] + # --> for 3 (GREP), one must also enable GREP in the Makefile + # --> for 2/3 (table/GREP), edit the corresponding parameters below too + + +# GREP +GREP_CENTER_METHOD 3 # center of the GREP profile [3] + # (1=box center, 2=density maximum, 3=potential minimum, 4=CoM) +GREP_MAXITER 1000 # number of iterations for constructing mass_TOV and Gamma_TOV [1000] +GREP_LOGBIN 1 # log/linear bins [true] + # (true=logarithmic, false=linear) +GREP_LOGBINRATIO 1.01 # ratio of adjacent log bins [1.25] +GREP_MAXRADIUS -1.0 # maximum radius in the radial profile [-1.0] + # (<=0=sqrt(3.0)*BOX_SIZE) +GREP_MINBINSIZE -1.0 # minimum bin size [-1.0] + # (<=0=amr->dh[MAX_LEVEL]) +GREP_OPT_FIXUP 1 # correct the profiles after grid allocation/deallocation [1] +GREP_OPT_PRES 2 # pressure computation scheme [2] + # (1=individual cell, 2=binned data) + + +# initialization +OPT__INIT 1 # initialization option: (1=FUNCTION, 2=RESTART, 3=FILE->"UM_IC") +OPT__INIT_BFIELD_BYVECPOT 2 # initialize the magnetic field from vector potential + # (0=off, 1=external disk file named "B_IC", see tool/inits/gen_vec_pot.py for example, 2=function) [0] ##MHD ONLY## +RESTART_LOAD_NRANK 1 # number of parallel I/O (i.e., number of MPI ranks) for restart [1] +OPT__RESTART_RESET 0 # reset some simulation status parameters (e.g., current step and time) during restart [0] +OPT__UM_IC_LEVEL 0 # starting AMR level in UM_IC [0] +OPT__UM_IC_NLEVEL 1 # number of AMR levels UM_IC [1] --> edit "Input__UM_IC_RefineRegion" if >1 +OPT__UM_IC_NVAR -1 # number of variables in UM_IC: (1~NCOMP_TOTAL; <=0=auto) [HYDRO=5+passive/ELBDM=2] +OPT__UM_IC_FORMAT 1 # data format of UM_IC: (1=vzyx, 2=zyxv; row-major and v=field) [1] +OPT__UM_IC_DOWNGRADE 1 # downgrade UM_IC from level OPT__UM_IC_LEVEL to 0 [1] +OPT__UM_IC_REFINE 1 # refine UM_IC from level OPT__UM_IC_LEVEL to MAX_LEVEL [1] +OPT__UM_IC_LOAD_NRANK 1 # number of parallel I/O (i.e., number of MPI ranks) for loading UM_IC [1] +OPT__INIT_RESTRICT 1 # restrict all data during the initialization [1] +OPT__INIT_GRID_WITH_OMP 1 # enable OpenMP when assigning the initial condition of each grid patch [1] +OPT__GPUID_SELECT -1 # GPU ID selection mode: (-3=Laohu, -2=CUDA, -1=MPI rank, >=0=input) [-1] +INIT_SUBSAMPLING_NCELL 4 # perform sub-sampling during initialization: (0=off, >0=# of sub-sampling cells) [0] +OPT__FFTW_STARTUP -1 # initialise fftw plans: (-1=auto, 0=ESTIMATE, 1=MEASURE, 2=PATIENT (only FFTW3)) [-1] + + +# interpolation schemes: (-1=auto, 1=MinMod-3D, 2=MinMod-1D, 3=vanLeer, 4=CQuad, 5=Quad, 6=CQuar, 7=Quar) +OPT__INT_TIME 1 # perform "temporal" interpolation for OPT__DT_LEVEL == 2/3 [1] +OPT__INT_PRIM 1 # switch to primitive variables when the interpolation on conserved variables fails [1] ##HYDRO ONLY## +OPT__INT_PHASE 1 # interpolation on phase (does not support MinMod-1D) [1] ##ELBDM ONLY## +OPT__FLU_INT_SCHEME -1 # ghost-zone fluid variables for the fluid solver [-1] +OPT__REF_FLU_INT_SCHEME -1 # newly allocated fluid variables during grid refinement [-1] +OPT__MAG_INT_SCHEME 4 # ghost-zone magnetic field for the MHD solver (2,3,4,6 only) [4] +OPT__REF_MAG_INT_SCHEME 4 # newly allocated magnetic field during grid refinement (2,3,4,6 only) [4] +OPT__POT_INT_SCHEME 5 # ghost-zone potential for the Poisson solver (only supports 4 & 5) [4] +OPT__RHO_INT_SCHEME 4 # ghost-zone mass density for the Poisson solver [4] +OPT__GRA_INT_SCHEME 5 # ghost-zone potential for the gravity solver (for UNSPLIT_GRAVITY as well) [4] +OPT__REF_POT_INT_SCHEME 5 # newly allocated potential during grid refinement [4] +INT_MONO_COEFF 2.0 # coefficient for ensuring the interpolation monotonicity (1.0~4.0) [2.0] +INT_MONO_COEFF_B 2.0 # coefficient for ensuring the interpolation monotonicity of B field (1.0~4.0) [2.0] ##MHD ONLY## +MONO_MAX_ITER 15 # maximum number of iterations to reduce INT_MONO_COEFF when interpolation fails (0=off) [10] +INT_OPP_SIGN_0TH_ORDER 1 # switch to 0th-order interpolation if adjacent values change signs [HYDRO:1; ELBDM:0] + + +# data dump +OPT__OUTPUT_TOTAL 1 # output the simulation snapshot: (0=off, 1=HDF5, 2=C-binary) [1] +OPT__OUTPUT_PART 0 # output a single line or slice: (0=off, 1=xy, 2=yz, 3=xz, 4=x, 5=y, 6=z, 7=diag) [0] +OPT__OUTPUT_TEXT_FORMAT_FLT %24.16e # string format of output text files [%24.16e] +OPT__OUTPUT_USER 0 # output the user-specified data -> edit "Output_User.cpp" [0] +OPT__OUTPUT_PAR_MODE 0 # output the particle data: (0=off, 1=text-file, 2=C-binary) [0] ##PARTICLE ONLY## +OPT__OUTPUT_BASEPS 0 # output the base-level power spectrum [0] +OPT__OUTPUT_BASE 0 # only output the base-level data [0] ##OPT__OUTPUT_PART ONLY## +OPT__OUTPUT_POT 1 # output gravitational potential [1] ##OPT__OUTPUT_TOTAL ONLY## +OPT__OUTPUT_PAR_DENS 1 # output the particle or total mass density on grids: + # (0=off, 1=particle mass density, 2=total mass density) [1] ##OPT__OUTPUT_TOTAL ONLY## +OPT__OUTPUT_CC_MAG 1 # output **cell-centered** magnetic field (necessary for yt analysis) [1] ##MHD ONLY## +OPT__OUTPUT_PRES 1 # output gas pressure [0] ##HYDRO ONLY## +OPT__OUTPUT_TEMP 1 # output gas temperature [0 (HD) or 1 (SRHD)] ##HYDRO ONLY## +OPT__OUTPUT_ENTR 1 # output gas entropy [0] ##HYDRO ONLY## +OPT__OUTPUT_CS 0 # output sound speed [0] ##HYDRO ONLY## +OPT__OUTPUT_DIVVEL 0 # output divergence(velocity) [0] ##HYDRO ONLY## +OPT__OUTPUT_MACH 0 # output mach number [0] ##HYDRO ONLY## +OPT__OUTPUT_DIVMAG 0 # output |divergence(B)*dh/|B|| [0] ##MHD ONLY## +OPT__OUTPUT_LORENTZ 0 # output Lorentz factor [0] ##SRHD ONLY## +OPT__OUTPUT_3VELOCITY 0 # output 3-velocities [0] ##SRHD ONLY## +OPT__OUTPUT_ENTHALPY 1 # output reduced enthalpy [1] ##SRHD ONLY## +OPT__OUTPUT_USER_FIELD 0 # output user-defined derived fields [0] -> edit "Flu_DerivedField_User.cpp" +OPT__OUTPUT_MODE 2 # (1=const step, 2=const dt, 3=dump table) -> edit "Input__DumpTable" for 3 +OPT__OUTPUT_RESTART 0 # output data immediately after restart [0] +OUTPUT_STEP 5 # output data every OUTPUT_STEP step ##OPT__OUTPUT_MODE==1 ONLY## +OUTPUT_DT 10.0 # output data every OUTPUT_DT time interval ##OPT__OUTPUT_MODE==2 ONLY## +OUTPUT_WALLTIME -1.0 # output data every OUTPUT_WALLTIME walltime (<=0.0=off) [-1.0] +OUTPUT_WALLTIME_UNIT 0 # unit of OUTPUT_WALLTIME (0=second, 1=minute, 2=hour, 3=day) [0] +OUTPUT_PART_X -1.0 # x coordinate for OPT__OUTPUT_PART [-1.0] +OUTPUT_PART_Y -1.0 # y coordinate for OPT__OUTPUT_PART [-1.0] +OUTPUT_PART_Z -1.0 # z coordinate for OPT__OUTPUT_PART [-1.0] +INIT_DUMPID -1 # set the first dump ID (<0=auto) [-1] + + +# yt inline analysis (SUPPORT_LIBYT only) +YT_SCRIPT yt_inline # yt inline analysis script (do not include the ".py" file extension) +YT_VERBOSE 1 # verbose level of yt (0=off, 1=info, 2=warning, 3=debug) [1] +YT_FIG_BASENAME Fig # figure basename [Fig] +YT_JUPYTER_USE_CONNECTION_FILE 0 # use user-provided connection file when using libyt Jupyter UI [0] + + +# miscellaneous +OPT__VERBOSE 0 # output the simulation progress in detail [0] +OPT__TIMING_BARRIER -1 # synchronize before timing -> more accurate, but may slow down the run (<0=auto) [-1] +OPT__TIMING_BALANCE 0 # record the max/min elapsed time in various code sections for checking load balance [0] +OPT__TIMING_MPI 0 # record the MPI bandwidth achieved in various code sections [0] ##LOAD_BALANCE ONLY## +OPT__RECORD_NOTE 1 # take notes for the general simulation info [1] +OPT__RECORD_UNPHY 1 # record the number of cells with unphysical results being corrected [1] +OPT__RECORD_MEMORY 1 # record the memory consumption [1] +OPT__RECORD_PERFORMANCE 1 # record the code performance [1] +OPT__MANUAL_CONTROL 1 # support manually dump data, stop run, or pause run during the runtime + # (by generating the file DUMP_GAMER_DUMP, STOP_GAMER_STOP, PAUSE_GAMER_PAUSE, respectively) [1] +OPT__RECORD_CENTER 1 # record the position of maximum density, minimum potential, and center of mass [0] +COM_CEN_X -1.0 # x coordinate as an initial guess for determining center of mass (if one of COM_CEN_X/Y/Z < 0 -> peak density position x) [-1.0] +COM_CEN_Y -1.0 # y coordinate as an initial guess for determining center of mass (if one of COM_CEN_X/Y/Z < 0 -> peak density position y) [-1.0] +COM_CEN_Z -1.0 # z coordinate as an initial guess for determining center of mass (if one of COM_CEN_X/Y/Z < 0 -> peak density position z) [-1.0] +COM_MAX_R -1.0 # maximum radius for determining center of mass (<0=auto -> __FLT_MAX__) [-1.0] +COM_MIN_RHO 0.0 # minimum density for determining center of mass (must >= 0.0) [0.0] +COM_TOLERR_R -1.0 # maximum tolerated error of deviation in radius during the iterations of determining the center of mass (<0=auto -> amr->dh[MAX_LEVEL]) [-1.0] +COM_MAX_ITER 10 # maximum number of iterations for determining the center of mass (must >= 1) [10] +OPT__RECORD_USER 1 # record the user-specified info -> edit "Aux_RecordUser.cpp" [0] +OPT__OPTIMIZE_AGGRESSIVE 0 # apply aggressive optimizations (experimental) [0] +OPT__SORT_PATCH_BY_LBIDX 1 # sort patches to improve bitwise reproducibility [SERIAL:0, LOAD_BALACNE:1] + + +# checks +OPT__CK_REFINE 0 # check the grid refinement [0] +OPT__CK_PROPER_NESTING 0 # check the proper-nesting condition [0] +OPT__CK_CONSERVATION 1 # check the conservation law [0] +ANGMOM_ORIGIN_X -1.0 # x coordinate of the origin for angular momentum (<0=auto -> BoxCenter) [-1.0] +ANGMOM_ORIGIN_Y -1.0 # y coordinate of the origin for angular momentum (<0=auto -> BoxCenter) [-1.0] +ANGMOM_ORIGIN_Z -1.0 # z coordinate of the origin for angular momentum (<0=auto -> BoxCenter) [-1.0] +OPT__CK_NORMALIZE_PASSIVE 0 # check the normalization of passive scalars [0] ##OPT__NORMALIZE_PASSIVE ONLY## +OPT__CK_RESTRICT 0 # check the data restriction [0] +OPT__CK_FINITE 0 # check if all variables are finite [0] +OPT__CK_PATCH_ALLOCATE 0 # check if all patches are properly allocated [0] +OPT__CK_FLUX_ALLOCATE 0 # check if all flux arrays are properly allocated [0] ##HYDRO and ELBDM ONLY## +OPT__CK_NEGATIVE 0 # check the negative values: (0=off, 1=density, 2=pressure and entropy, 3=both) [0] ##HYDRO ONLY## +OPT__CK_MEMFREE 1.0 # check the free memory in GB (0=off, >0=threshold) [1.0] +OPT__CK_PARTICLE 0 # check the particle allocation [0] +OPT__CK_INTERFACE_B 0 # check the consistency of patch interface B field [0] ##MHD ONLY## +OPT__CK_DIVERGENCE_B 0 # check the divergence-free constraint on B field (0=off, 1=on, 2=on+verbose) [0] ##MHD ONLY## +OPT__CK_INPUT_FLUID 0 # check the input data of the fluid solver [0] + diff --git a/example/test_problem/Hydro/CCSN/Input__Parameter.Lightbulb b/example/test_problem/Hydro/CCSN/Input__Parameter.Lightbulb index 3e3b408fd..99b335db5 100644 --- 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 index a135d6de5..3db672a4b 100644 --- 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 index 7d28474c9..319aa1a09 100644 --- 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 new file mode 100644 index 000000000..188ddadcf --- /dev/null +++ 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 index 9c10e1d7c..c864626af 100644 --- 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 index d0de98e6e..64e8f9715 100644 --- 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 index 52bf1581b..770422af2 100644 --- 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 index 6993c8951..f07cf9b05 100644 --- 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 index b8c427918..a7b0bc3bb 100644 --- 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 index 07f14fc5a..dca49a130 100644 --- 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 index e34e2c15b..f4effe4e3 100644 --- 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 index a5200424d..8e318a9c3 100644 --- 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 diff --git a/include/Macro.h b/include/Macro.h index d114f8d3e..2b3d40bdd 100644 --- a/include/Macro.h +++ b/include/Macro.h @@ -122,6 +122,7 @@ #define LIGHTBULB 1 #define IDSA 2 #define M1 3 +#define LEAKAGE 4 // ELBDM schemes @@ -217,13 +218,22 @@ # define NCOMP_PASSIVE_BUILTIN1 0 # endif -// electron fraction (Ye), neutrino heating/cooling rate, and temperature initial guess (TEMP_IG) +// electron fraction (YE), neutrino heating/cooling rate (DEDT_NU), Ye change rate (DYEDT_NU, leakage scheme only), +// and temperature initial guess (TEMP_IG) # if ( EOS == EOS_NUCLEAR ) +# if ( NEUTRINO_SCHEME == LEAKAGE ) +# if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) +# define NCOMP_PASSIVE_BUILTIN2 4 +# else +# define NCOMP_PASSIVE_BUILTIN2 3 +# endif +# else # if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) # define NCOMP_PASSIVE_BUILTIN2 3 # else # define NCOMP_PASSIVE_BUILTIN2 2 # endif +# endif // # if ( NEUTRINO_SCHEME == LEAKAGE ) ... else ... # else # define NCOMP_PASSIVE_BUILTIN2 0 # endif @@ -341,12 +351,17 @@ # if ( EOS == EOS_NUCLEAR ) # define YE ( PASSIVE_NEXT_IDX2 ) +# if ( NEUTRINO_SCHEME == LEAKAGE ) +# define DYEDT_NU ( YE - 1 ) +# define DEDT_NU ( YE - 2 ) +# else # define DEDT_NU ( YE - 1 ) +# endif # if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) -# define TEMP_IG ( YE - 2 ) -# define PASSIVE_NEXT_IDX3 ( YE - 3 ) +# define TEMP_IG ( DEDT_NU - 1 ) +# define PASSIVE_NEXT_IDX3 ( DEDT_NU - 2 ) # else -# define PASSIVE_NEXT_IDX3 ( YE - 2 ) +# define PASSIVE_NEXT_IDX3 ( DEDT_NU - 1 ) # endif # else # define PASSIVE_NEXT_IDX3 ( PASSIVE_NEXT_IDX2 ) @@ -389,16 +404,21 @@ # endif # if ( EOS == EOS_NUCLEAR ) -# define FLUX_YE ( FLUX_NEXT_IDX2 ) -# define FLUX_DEDT_NU ( FLUX_YE - 1 ) +# define FLUX_YE ( FLUX_NEXT_IDX2 ) +# if ( NEUTRINO_SCHEME == LEAKAGE ) +# define FLUX_DYEDT_NU ( FLUX_YE - 1 ) +# define FLUX_DEDT_NU ( FLUX_YE - 2 ) +# else +# define FLUX_DEDT_NU ( FLUX_YE - 1 ) +# endif # if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) -# define FLUX_TEMP_IG ( FLUX_YE - 2 ) -# define FLUX_NEXT_IDX3 ( FLUX_YE - 3 ) +# define FLUX_TEMP_IG ( FLUX_DEDT_NU - 2 ) +# define FLUX_NEXT_IDX3 ( FLUX_DEDT_NU - 3 ) # else -# define FLUX_NEXT_IDX3 ( FLUX_YE - 2 ) +# define FLUX_NEXT_IDX3 ( FLUX_DEDT_NU - 2 ) # endif # else -# define FLUX_NEXT_IDX3 ( FLUX_NEXT_IDX2 ) +# define FLUX_NEXT_IDX3 ( FLUX_NEXT_IDX2 ) # endif #endif // #if ( NCOMP_PASSIVE > 0 ) @@ -424,10 +444,13 @@ # endif # if ( EOS == EOS_NUCLEAR ) -# define _YE ( 1L << YE ) -# define _DEDT_NU ( 1L << DEDT_NU ) +# define _YE ( 1L << YE ) +# if ( NEUTRINO_SCHEME == LEAKAGE ) +# define _DYEDT_NU ( 1L << DYEDT_NU ) +# endif +# define _DEDT_NU ( 1L << DEDT_NU ) # if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) -# define _TEMP_IG ( 1L << TEMP_IG ) +# define _TEMP_IG ( 1L << TEMP_IG ) # endif # endif @@ -461,10 +484,13 @@ # endif # if ( EOS == EOS_NUCLEAR ) -# define _FLUX_YE ( 1L << FLUX_YE ) -# define _FLUX_DEDT_NU ( 1L << FLUX_DEDT_NU ) +# define _FLUX_YE ( 1L << FLUX_YE ) +# if ( NEUTRINO_SCHEME == LEAKAGE ) +# define _FLUX_DYEDT_NU ( 1L << FLUX_DYEDT_NU ) +# endif +# define _FLUX_DEDT_NU ( 1L << FLUX_DEDT_NU ) # if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) -# define _FLUX_TEMP_IG ( 1L << FLUX_TEMP_IG ) +# define _FLUX_TEMP_IG ( 1L << FLUX_TEMP_IG ) # endif # endif @@ -973,9 +999,11 @@ #if ( MODEL == HYDRO ) # define SRC_NAUX_DLEP 8 // SrcTerms.Dlep_AuxArray_Flt/Int[] # define SRC_NAUX_LIGHTBULB 2 // SrcTerms.Lightbulb_AuxArray_Flt/Int[] +# define SRC_NAUX_LEAKAGE 11 // SrcTerms.Leakage_AuxArray_Flt/Int[] #else # define SRC_NAUX_DLEP 0 # define SRC_NAUX_LIGHTBULB 0 +# define SRC_NAUX_LEAKAGE 0 #endif # define SRC_NAUX_USER 10 // SrcTerms.User_AuxArray_Flt/Int[] diff --git a/include/Prototype.h b/include/Prototype.h index 18d99beb7..47d18806f 100644 --- a/include/Prototype.h +++ b/include/Prototype.h @@ -46,6 +46,10 @@ void Aux_FindExtrema( Extrema_t *Extrema, const ExtremaMode_t Mode, const int Mi const PatchType_t PatchType ); void Aux_FindWeightedAverageCenter( double WeightedAverageCenter[], const double Center_ref[], const double MaxR, const double MinWD, const long WeightingDensityField, const double TolErrR, const int MaxIter, double *FinaldR, int *FinalNIter ); +void Aux_ComputeRay( Profile_t *Ray[], const double Center[], const double Edge[], + const int NRadius_Linear, const int NRadius, const int NTheta, const int NPhi, + const double BinSize_Linear, const double MaxRad_Linear, const double MaxRad, + const long TVarBitIdx[], const int NProf, const double PrepTimeIn ); #ifndef SERIAL void Aux_Record_BoundaryPatch( const int lv, int *NList, int **IDList, int **PosList ); #endif diff --git a/include/SrcTerms.h b/include/SrcTerms.h index fa4d49cee..79dddaf0a 100644 --- a/include/SrcTerms.h +++ b/include/SrcTerms.h @@ -27,6 +27,7 @@ typedef void (*SrcFunc_t)( real fluid[], const real B[], // Data Member : Any : True if at least one of the source terms is activated // Deleptonization : SRC_DELEPTONIZATION // Lightbulb : SRC_LIGHTBULB +// Leakage : SRC_LEAKAGE // User : SRC_USER // BoxCenter : Simulation box center // Unit_* : Code units @@ -35,8 +36,6 @@ typedef void (*SrcFunc_t)( real fluid[], const real B[], // *_AuxArrayDevPtr_* : Auxiliary array pointers // --> For GPU, these pointers store the addresses of constant memory arrays, // which should NOT be used by host -// Lightbulb_Lnue : Electron neutrino luminosity in erg/s -// Lightbulb_Tnue : Electron neutrino temperature in MeV // // Method : None --> It seems that CUDA does not support functions in a struct //------------------------------------------------------------------------------------------------------- @@ -46,6 +45,7 @@ struct SrcTerms_t bool Any; bool Deleptonization; bool Lightbulb; + bool Leakage; bool User; double BoxCenter[3]; @@ -87,7 +87,31 @@ struct SrcTerms_t int *Lightbulb_AuxArrayDevPtr_Int; double Lightbulb_Lnue; double Lightbulb_Tnue; + +// leakage + SrcFunc_t Leakage_FuncPtr; + SrcFunc_t Leakage_CPUPtr; +# ifdef GPU + SrcFunc_t Leakage_GPUPtr; # endif + double *Leakage_AuxArrayDevPtr_Flt; + int *Leakage_AuxArrayDevPtr_Int; + int Leakage_NRadius; + int Leakage_NTheta; + int Leakage_NPhi; + double Leakage_BinSize_Radius; + double Leakage_RadiusMax; + double Leakage_RadiusMin_Log; + bool Leakage_NuHeat; + double Leakage_NuHeat_Fac; + int Leakage_Opt_Temp; + real *Leakage_Radius_DevPtr; + real *Leakage_tau_DevPtr; + real *Leakage_chi_DevPtr; + real *Leakage_HeatFlux_DevPtr; + real *Leakage_HeatERms_DevPtr; + real *Leakage_HeatEAve_DevPtr; +# endif // if ( MODEL == HYDRO ) // user-specified source term SrcFunc_t User_FuncPtr; diff --git a/include/Typedef.h b/include/Typedef.h index f71a461b6..e55d56835 100644 --- a/include/Typedef.h +++ b/include/Typedef.h @@ -93,7 +93,7 @@ const TestProbID_t TESTPROB_HYDRO_CR_SOUNDWAVE = 20, TESTPROB_HYDRO_CR_SHOCKTUBE = 21, TESTPROB_HYDRO_CR_DIFFUSION = 23, - TESTPROB_HYDRO_CCSN = 24, + TESTPROB_HYDRO_CCSN = 50, TESTPROB_HYDRO_BARRED_POT = 51, TESTPROB_HYDRO_JET_ICM_WALL = 52, TESTPROB_HYDRO_CDM_LSS = 100, @@ -593,6 +593,18 @@ const LoadParaMode_t LOAD_HDF5_OUTPUT = 2; +// leakage source term +typedef int Leakage_TempScheme_t; +const Leakage_TempScheme_t + LEAK_TEMP_INDCELL = 1, + LEAK_TEMP_BINDATA = 2; + +typedef int Leakage_Mode_t; +const Leakage_Mode_t + LEAK_MODE_EVOLVE = 1, + LEAK_MODE_RECORD = 2; + + // function pointers typedef real (*EoS_GUESS_t) ( const real Con[], real* const Constant, const double AuxArray_Flt[], const int AuxArray_Int[], const real *const Table[EOS_NTABLE_MAX] ); diff --git a/src/Auxiliary/Aux_Check_Parameter.cpp b/src/Auxiliary/Aux_Check_Parameter.cpp index 8836dd426..bc5f48f2a 100644 --- a/src/Auxiliary/Aux_Check_Parameter.cpp +++ b/src/Auxiliary/Aux_Check_Parameter.cpp @@ -925,8 +925,8 @@ void Aux_Check_Parameter() # endif # endif // #ifdef BAROTROPIC_EOS ... else ... -# if ( defined NEUTRINO_SCHEME && NEUTRINO_SCHEME != LIGHTBULB && NEUTRINO_SCHEME != IDSA && NEUTRINO_SCHEME != M1 ) -# error : ERROR : unsupported neutrino updating scheme (LIGHTBULB/IDSA/M1) !! +# if ( defined NEUTRINO_SCHEME && NEUTRINO_SCHEME != LIGHTBULB && NEUTRINO_SCHEME != IDSA && NEUTRINO_SCHEME != M1 && NEUTRINO_SCHEME != LEAKAGE ) +# error : ERROR : unsupported neutrino updating scheme (LIGHTBULB/IDSA/M1/LEAKAGE) !! # endif if ( OPT__1ST_FLUX_CORR != FIRST_FLUX_CORR_NONE ) @@ -1866,6 +1866,38 @@ void Aux_Check_Parameter() if ( SrcTerms.Lightbulb ) Aux_Error( ERROR_INFO, "SRC_LIGHTBULB is only supported in HYDRO and must work with EOS_NUCLEAR !!\n" ); + + if ( SrcTerms.Leakage ) + Aux_Error( ERROR_INFO, "SRC_LEAKAGE is only supported in HYDRO and must work with EOS_NUCLEAR !!\n" ); +# endif + + if ( SrcTerms.Lightbulb && SrcTerms.Leakage ) + Aux_Error( ERROR_INFO, "SRC_LIGHTBULB and SRC_LEAKAGE cannot be enabled simultaneously !!\n" ); + +# if ( NEUTRINO_SCHEME != LIGHTBULB ) + if ( SrcTerms.Lightbulb ) + Aux_Error( ERROR_INFO, "NEUTRINO_SCHEME != LIGHTBULB for the lightbulb scheme !!\n" ); +# endif + +# if ( NEUTRINO_SCHEME != LEAKAGE ) + if ( SrcTerms.Leakage ) + Aux_Error( ERROR_INFO, "NEUTRINO_SCHEME != LEAKAGE for the leakage scheme !!\n" ); +# else + if ( SrcTerms.Leakage_RadiusMax < SrcTerms.Leakage_BinSize_Radius ) + Aux_Error( ERROR_INFO, "%s (%14.7e) <= %s (%14.7e) !!\n", + "SRC_LEAKAGE_RADIUSMAX", SrcTerms.Leakage_RadiusMax, "SRC_LEAKAGE_BINSIZE_RADIUS", SrcTerms.Leakage_BinSize_Radius ); + + if ( SrcTerms.Leakage_RadiusMax < SrcTerms.Leakage_RadiusMin_Log ) + Aux_Error( ERROR_INFO, "%s (%14.7e) <= %s (%14.7e) !!\n", + "SRC_LEAKAGE_RADIUSMAX", SrcTerms.Leakage_RadiusMax, "SRC_LEAKAGE_RADIUSMIN_LOG", SrcTerms.Leakage_RadiusMin_Log ); + + if ( SrcTerms.Leakage_NPhi > 1 && SrcTerms.Leakage_NPhi % 2 ) + Aux_Error( ERROR_INFO, "(%d) must be 1 or a multiple of 2 !!\n", + "SRC_LEAKAGE_NPHI", SrcTerms.Leakage_NPhi ); + + if ( int( SrcTerms.Leakage_RadiusMin_Log / SrcTerms.Leakage_BinSize_Radius ) > SrcTerms.Leakage_NRadius ) + Aux_Error( ERROR_INFO, "%s (%14.7e) / %s (%14.7e) > %s (%d) !!\n", + "SRC_LEAKAGE_RADIUSMIN_LOG", SrcTerms.Leakage_RadiusMin_Log, "SRC_LEAKAGE_BINSIZE_RADIUS", SrcTerms.Leakage_BinSize_Radius, "SRC_LEAKAGE_NRADIUS", SrcTerms.Leakage_NRadius ); # endif // warning diff --git a/src/Auxiliary/Aux_ComputeRay.cpp b/src/Auxiliary/Aux_ComputeRay.cpp new file mode 100644 index 000000000..882798ce8 --- /dev/null +++ b/src/Auxiliary/Aux_ComputeRay.cpp @@ -0,0 +1,732 @@ +#include "GAMER.h" + +static const double TWOPI = 2.0 * M_PI; + + +static void Aux_GetMinMax_Vertex_Radius( const double x, const double y, const double z, const double HalfWidth, + double *RadMin, double *RadMax ); + +static void Aux_GetMinMax_Vertex_Angle( const double x, const double y, const double z, const double HalfWidth, + double *ThtMin, double *ThtMax, double *PhiMin, double *PhiMax ); + +extern void SetTempIntPara( const int lv, const int Sg0, const double PrepTime, const double Time0, const double Time1, + bool &IntTime, int &Sg, int &Sg_IntT, real &Weighting, real &Weighting_IntT ); + + + + +//------------------------------------------------------------------------------------------------------- +// Function : Aux_ComputeRay +// Description : Maps the target field(s) from cells to rays +// +// Note : 1. Results will be stored in the input "Prof" object +// --> with index "Idx_Radius + NRadius * ( Idx_Theta + NTheta * Idx_Phi )" +// 2. Support hybrid OpenMP/MPI parallelization +// --> All ranks will share the same ray data after invoking this function +// 3. Weighting of each cell: +// --> Cell mass : gas velocity, gravitational potential +// Cell volume: other fields +// 4. Support computing multiple fields +// --> The order of fields to be returned follows TVarBitIdx[] +// 5. This routine is thread-unsafe when the temporal interpolation set by PrepTime and OPT__INT_TIME +// are inconsistent with each other +// --> But it shouldn't be a big issue since this routine itself has been parallelized with OpenMP +// +// Parameter : Ray : Profile_t object array to store the results +// Center : Target center coordinates +// Edge : Edge of bins in the radial direction +// NRadius_Linear : Number of linear bins in the radial direction +// NRadius : Number of bins in the radial direction +// NTheta : Number of bins in the polar direction +// NPhi : Number of bins in the azimuthal direction +// BinSize_Linear : Bin size of linear bins in the radial direction +// MaxRad_Linear : Maximum radius of linear bins +// MaxRad : Maximum radius in each ray profile +// TVarBitIdx : Bitwise indices of target variables for computing the rays +// --> Supported indices (defined in Macro.h): +// HYDRO : _DENS, _MOMX, _MOMY, _MOMZ, _ENGY, _VELX, _VELY, _VELZ, _VELR, +// _PRES, _TEMP, _ENTR, _EINT +// [, _DUAL, _CRAY, _POTE, __MAGX_CC, _MAGY_CC, _MAGZ_CC, _MAGE_CC] +// ELBDM : _DENS, _REAL, _IMAG [, _POTE] +// --> All fields supported by Prepare_PatchData() are also supported here +// NProf : Number of target fields in Ray +// PrepTimeIn : Target physical time to prepare data +// --> If PrepTimeIn<0, turn off temporal interpolation and always use the most recent data +// +// Example : const double Center[3] = { amr->BoxCenter[0], amr->BoxCenter[1], amr->BoxCenter[2] }; +// const int NRadius_Linear = 256; +// const int NRadius = 256; +// const int NTheta = 1; +// const int NPhi = 1; +// const double BinSize_Linear = amr->dh[MAX_LEVEL]; +// const double MaxRad_Linear = 0.5*amr->BoxSize[0]; +// const double MaxRad = 0.5*amr->BoxSize[0]; +// const long TVarBitIdx[] = { _DENS, _TEMP }; +// const int NProf = 2; +// const double PrepTimeIn = -1.0; +// double Edge[NRadius+1]; +// +// Profile_t Ray_Dens, Ray_Temp; +// Profile_t *Ray[NProf] = { &Ray_Dens, &Ray_Temp }; +// +// for (int i=0; i all fields supported by Prepare_PatchData() should be supported here + long SupportedFields = ( _TOTAL | _DERIVED ); +# ifdef GRAVITY + SupportedFields |= _POTE; +# endif +# ifdef PARTICLE + SupportedFields |= _PAR_DENS; + SupportedFields |= _TOTAL_DENS; +# endif + + for (int p=0; pNBin = NBin; + Ray[p]->AllocateMemory(); + } + + +// allocate memory for the per-thread arrays +# ifdef OPENMP + const int NT = OMP_NTHREAD; +# else + const int NT = 1; +# endif + + double ***OMP_Data=NULL, ***OMP_Weight=NULL; + long ***OMP_NCell=NULL; + + Aux_AllocateArray3D( OMP_Data, NProf, NT, NBin ); + Aux_AllocateArray3D( OMP_Weight, NProf, NT, NBin ); + Aux_AllocateArray3D( OMP_NCell, NProf, NT, NBin ); + +// initialize profile arrays +// --> use memset() instead to improve performance +// OMP_Data and OMP_Weight are initialized during their first assignment + for (int p=0; pBoxSize[0], 0.5*amr->BoxSize[1], 0.5*amr->BoxSize[2] }; + const bool Periodic[3] = { OPT__BC_FLU[0] == BC_FLU_PERIODIC, + OPT__BC_FLU[2] == BC_FLU_PERIODIC, + OPT__BC_FLU[4] == BC_FLU_PERIODIC }; + + +// temporarily overwrite OPT__INT_TIME +// --> necessary because SetTempIntPara() called by Prepare_PatchData() relies on OPT__INT_TIME +// --> must restore it before exiting this routine +// --> note that modifying OPT__INT_TIME renders this routine thread-unsafe +//###REVISE: make temporal interpolation a function parameter in Prepare_PatchData() to solve this thread-safety issue + const bool IntTimeBackup = OPT__INT_TIME; + OPT__INT_TIME = ( PrepTimeIn >= 0.0 ) ? true : false; + + +// loop over all levels + for (int lv=0; lv<=MAX_LEVEL; lv++) + { + if ( NPatchTotal[lv] == 0 ) continue; + + const double dh = amr->dh[lv]; + const double dh_half = 0.5*amr->dh[lv]; + const double dv = CUBE( dh ); + + +// determine the temporal interpolation parameters +// --> mainly for computing cell mass for weighting; Prepare_PatchData() needs PrepTime + const int FluSg0 = amr->FluSg[lv]; + const double PrepTime = ( PrepTimeIn >= 0.0 ) ? PrepTimeIn : amr->FluSgTime[lv][FluSg0]; + + bool FluIntTime; + int FluSg, FluSg_IntT; + real FluWeighting, FluWeighting_IntT; + + SetTempIntPara( lv, FluSg0, PrepTime, amr->FluSgTime[lv][FluSg0], amr->FluSgTime[lv][1-FluSg0], + FluIntTime, FluSg, FluSg_IntT, FluWeighting, FluWeighting_IntT ); + + +// initialize the particle density array (rho_ext) and collect particles to the target level +# ifdef PARTICLE + const bool TimingSendPar_No = false; + const bool JustCountNPar_No = false; +# ifdef LOAD_BALANCE + const bool PredictPos = amr->Par->PredictPos; + const bool SibBufPatch = true; + const bool FaSibBufPatch = true; +# else + const bool PredictPos = false; + const bool SibBufPatch = NULL_BOOL; + const bool FaSibBufPatch = NULL_BOOL; +# endif + + if ( NeedPar ) + { +// these two routines should NOT be put inside an OpenMP parallel region + Par_CollectParticle2OneLevel( lv, _PAR_MASS|_PAR_POSX|_PAR_POSY|_PAR_POSZ, _PAR_TYPE, PredictPos, + PrepTime, SibBufPatch, FaSibBufPatch, JustCountNPar_No, TimingSendPar_No ); + + Prepare_PatchData_InitParticleDensityArray( lv, PrepTime ); + } // if ( NeedPar ) +# endif // #ifdef PARTICLE + + +// different OpenMP threads and MPI processes first compute profiles independently +// --> their data will be combined later +# pragma omp parallel + { +# ifdef OPENMP + const int TID = omp_get_thread_num(); +# else + const int TID = 0; +# endif + +// use the "static" schedule for reproducibility +# pragma omp for schedule( static ) + for (int PID0=0; PID0NPatchComma[lv][1]; PID0+=8) + { +// skip non-leaf patches + bool SkipPatch[8], SkipPatchGroup=true; + + for (int LocalID=0; LocalID<8; LocalID++) + { + const int PID = PID0 + LocalID; + SkipPatch[LocalID] = false; + + if ( amr->patch[0][lv][PID]->son != -1 && lv != MAX_LEVEL ) SkipPatch[LocalID] = true; + if ( ! SkipPatch[LocalID] ) SkipPatchGroup = false; + } // for (int LocalID=0; LocalID<8; LocalID++) + + if ( SkipPatchGroup ) continue; + + +// store the bin indices associated with each cell +// --> do it before looping over all target fields to avoid redundant calculations + for (int LocalID=0; LocalID<8; LocalID++) + { + if ( SkipPatch[LocalID] ) continue; + + const int PID = PID0 + LocalID; + const double x0 = amr->patch[0][lv][PID]->EdgeL[0] + dh_half - Center[0]; + const double y0 = amr->patch[0][lv][PID]->EdgeL[1] + dh_half - Center[1]; + const double z0 = amr->patch[0][lv][PID]->EdgeL[2] + dh_half - Center[2]; + + for (int k=0; k +HalfBox[2] ) { dz -= amr->BoxSize[2]; } + else if ( dz < -HalfBox[2] ) { dz += amr->BoxSize[2]; } + } + for (int j=0; j +HalfBox[1] ) { dy -= amr->BoxSize[1]; } + else if ( dy < -HalfBox[1] ) { dy += amr->BoxSize[1]; } + } + for (int i=0; i +HalfBox[0] ) { dx -= amr->BoxSize[0]; } + else if ( dx < -HalfBox[0] ) { dx += amr->BoxSize[0]; } + } +// find the minimum/maximum radius among the eight cell vertices + double MinRad_Cell, MaxRad_Cell; + + Aux_GetMinMax_Vertex_Radius( dx, dy, dz, dh_half, &MinRad_Cell, &MaxRad_Cell ); + + + if ( MinRad_Cell >= MaxRad ) + { + Patch_BinMin_Rad[TID][LocalID][k][j][i] = CellSkip; + Patch_BinMax_Rad[TID][LocalID][k][j][i] = CellSkip; + } + + else + { +// find the minimum/maximum azimuthal angle and polar angle among the eight cell vertices + double MinTht_Cell, MaxTht_Cell, MinPhi_Cell, MaxPhi_Cell; + + Aux_GetMinMax_Vertex_Angle( dx, dy, dz, dh_half, &MinTht_Cell, &MaxTht_Cell, &MinPhi_Cell, &MaxPhi_Cell ); + + +// find the index range of bins overlapping with this cell +// --> use Cell_Skip for handling the vertex lies on the bin edge +// case (a): the center coordinate is inside the cell +// --> contribute to all rays + if ( ( fabs(dx) < dh_half ) && + ( fabs(dy) < dh_half ) && + ( fabs(dz) < dh_half ) ) + { + Patch_BinMin_Rad[TID][LocalID][k][j][i] = 0; + Patch_BinMax_Rad[TID][LocalID][k][j][i] = ( MaxRad_Cell <= MaxRad_Linear ) + ? int( MaxRad_Cell * _dr ) + : Mis_BinarySearch_Real( Edge, NRadius_Linear, NRadius-1, MaxRad_Cell ); + Patch_BinMin_Tht[TID][LocalID][k][j][i] = 0; + Patch_BinMax_Tht[TID][LocalID][k][j][i] = NTheta - 1; + Patch_BinMin_Phi[TID][LocalID][k][j][i] = 0; + Patch_BinMax_Phi[TID][LocalID][k][j][i] = NPhi - 1; + } + + else + { + Patch_BinMin_Rad[TID][LocalID][k][j][i] = ( MinRad_Cell <= MaxRad_Linear ) + ? int( MinRad_Cell * _dr ) + : Mis_BinarySearch_Real( Edge, NRadius_Linear, NRadius-1, MinRad_Cell ); + Patch_BinMax_Rad[TID][LocalID][k][j][i] = ( MaxRad_Cell <= MaxRad_Linear ) + ? int( MaxRad_Cell * _dr ) + : Mis_BinarySearch_Real( Edge, NRadius_Linear, NRadius-1, MaxRad_Cell ); + +// case (b): the cell intersects the reference z-axis +// --> contributes to rays at any azimuthal angle + if ( ( fabs(dx) < dh_half ) && ( fabs(dy) < dh_half ) ) + { + Patch_BinMin_Tht[TID][LocalID][k][j][i] = ( dz > 0.0 ) ? 0 : int( MinTht_Cell * _dtheta ); + Patch_BinMax_Tht[TID][LocalID][k][j][i] = ( dz > 0.0 ) ? int( MaxTht_Cell * _dtheta ) : NTheta - 1; + Patch_BinMin_Phi[TID][LocalID][k][j][i] = 0; + Patch_BinMax_Phi[TID][LocalID][k][j][i] = NPhi - 1; + } + + else + { + Patch_BinMin_Tht[TID][LocalID][k][j][i] = int( MinTht_Cell * _dtheta ); + Patch_BinMax_Tht[TID][LocalID][k][j][i] = int( MaxTht_Cell * _dtheta ); + +// case (c): the cell intersects the positive xz-plane (azimuthal discontinuity) +// --> handle the index range of the azimuthal angle carefully + if ( ( dx > dh_half ) && ( fabs(dy) < dh_half ) ) + { + MinPhi_Cell = atan2( dy - dh_half, dx - dh_half ) + TWOPI; + MaxPhi_Cell = atan2( dy + dh_half, dx - dh_half ); + +// add NPhi to enable looping later + Patch_BinMin_Phi[TID][LocalID][k][j][i] = int( MinPhi_Cell * _dphi ); + Patch_BinMax_Phi[TID][LocalID][k][j][i] = ( NPhi > 1 ) + ? int( MaxPhi_Cell * _dphi ) + NPhi + : int( MaxPhi_Cell * _dphi ); + } + +// case (d): general case + else + { + Patch_BinMin_Phi[TID][LocalID][k][j][i] = int( MinPhi_Cell * _dphi ); + Patch_BinMax_Phi[TID][LocalID][k][j][i] = int( MaxPhi_Cell * _dphi ); + } + } + } + +// handle the maximum bin indices carefully to prevent redundant mapping + double BinMax_Tht = MaxTht_Cell * _dtheta; + double BinMax_Phi = MaxPhi_Cell * _dphi; + + if ( Mis_CompareRealValue( floor(BinMax_Tht), BinMax_Tht, NULL, false ) ) + Patch_BinMax_Tht[TID][LocalID][k][j][i] -= 1; + + if ( Mis_CompareRealValue( floor(BinMax_Phi), BinMax_Phi, NULL, false ) ) + Patch_BinMax_Phi[TID][LocalID][k][j][i] -= 1; + } // if ( MaxRad_Cell > MaxRad ) ... else ... + }}} // i,j,k + } // for (int LocalID=0; LocalID<8; LocalID++) + + +// compute one field at a time + for (int p=0; ppatch[ FluSg ][lv][PID]->fluid; + const real (*FluidPtr_IntT)[PS1][PS1][PS1] = ( FluIntTime ) ? amr->patch[ FluSg_IntT ][lv][PID]->fluid : NULL; + + const double x0 = amr->patch[0][lv][PID]->EdgeL[0] + dh_half - Center[0]; + const double y0 = amr->patch[0][lv][PID]->EdgeL[1] + dh_half - Center[1]; + const double z0 = amr->patch[0][lv][PID]->EdgeL[2] + dh_half - Center[2]; + + for (int k=0; k +HalfBox[2] ) { dz -= amr->BoxSize[2]; } + else if ( dz < -HalfBox[2] ) { dz += amr->BoxSize[2]; } + } + for (int j=0; j +HalfBox[1] ) { dy -= amr->BoxSize[1]; } + else if ( dy < -HalfBox[1] ) { dy += amr->BoxSize[1]; } + } + for (int i=0; i +HalfBox[0] ) { dx -= amr->BoxSize[0]; } + else if ( dx < -HalfBox[0] ) { dx += amr->BoxSize[0]; } + } + + if ( Patch_BinMin_Rad[TID][LocalID][k][j][i] == CellSkip ) continue; + + const double r = sqrt( SQR(dx) + SQR(dy) + SQR(dz) ); + const real _Dens = (real)1.0 / FluidPtr [DENS][k][j][i]; + const real _Dens_IntT = ( FluIntTime ) ? (real)1.0 / FluidPtr_IntT[DENS][k][j][i] : NULL_REAL; + + real VelR; + if ( r == 0.0 ) { + VelR = (real)0.0; // take care of the corner case where the profile center coincides with a cell center + } + + else { + VelR = ( FluIntTime ) + ? ( FluWeighting *( FluidPtr [MOMX][k][j][i]*dx + + FluidPtr [MOMY][k][j][i]*dy + + FluidPtr [MOMZ][k][j][i]*dz )*_Dens + + FluWeighting_IntT*( FluidPtr_IntT[MOMX][k][j][i]*dx + + FluidPtr_IntT[MOMY][k][j][i]*dy + + FluidPtr_IntT[MOMZ][k][j][i]*dz )*_Dens_IntT ) / r + : ( FluidPtr [MOMX][k][j][i]*dx + + FluidPtr [MOMY][k][j][i]*dy + + FluidPtr [MOMZ][k][j][i]*dz )*_Dens / r; + } + + Patch_Data[TID][LocalID][k][j][i] = VelR; + }}} // i,j,k + } // for (int LocalID=0; LocalID<8; LocalID++) + break; // _VELR +# endif // #ifdef _VELR + + default: + const int NGhost = 0; + const int NPG = 1; + const bool IntPhase_No = false; + const real MinDens_No = -1.0; + const real MinPres_No = -1.0; + const real MinTemp_No = -1.0; + const real MinEntr_No = -1.0; + const bool DE_Consistency_Yes = true; + + Prepare_PatchData( lv, PrepTime, &Patch_Data[TID][0][0][0][0], NULL, NGhost, NPG, &PID0, + TVarBitIdx[p], _NONE, INT_NONE, INT_NONE, UNIT_PATCH, NSIDE_00, IntPhase_No, + OPT__BC_FLU, BC_POT_NONE, MinDens_No, MinPres_No, MinTemp_No, MinEntr_No, + DE_Consistency_Yes ); + break; // default + } // switch ( TVarBitIdx[p] ) + + +// set the weight field +//###REVISE: allow users to choose the weight field + int WeightField=-1; + + switch ( TVarBitIdx[p] ) + { +# ifdef _VELX + case _VELX : WeightField = WeightByMass; break; +# endif +# ifdef _VELY + case _VELY : WeightField = WeightByMass; break; +# endif +# ifdef _VELZ + case _VELZ : WeightField = WeightByMass; break; +# endif +# ifdef _VELR + case _VELR : WeightField = WeightByMass; break; +# endif +# ifdef _POTE + case _POTE : WeightField = WeightByMass; break; +# endif + default : WeightField = WeightByVolume; break; + } // switch ( TVarBitIdx[p] ) + + +// compute the ray profile + for (int LocalID=0; LocalID<8; LocalID++) + { + if ( SkipPatch[LocalID] ) continue; + + const int PID = PID0 + LocalID; + const real (*DensPtr )[PS1][PS1] = amr->patch[ FluSg ][lv][PID]->fluid[DENS]; + const real (*DensPtr_IntT)[PS1][PS1] = ( FluIntTime ) ? amr->patch[ FluSg_IntT ][lv][PID]->fluid[DENS] : NULL; + + for (int k=0; kNPatchComma[lv][1]; PID0+=8) + } // OpenMP parallel region + + +// free particle resources +// --> these two routines should NOT be put inside an OpenMP parallel region +# ifdef PARTICLE + if ( NeedPar ) + { + Par_CollectParticle2OneLevel_FreeMemory( lv, SibBufPatch, FaSibBufPatch ); + + Prepare_PatchData_FreeParticleDensityArray( lv ); + } +# endif + } // for (int lv=MinLv; lv<=MaxLv; lv++) + + +# pragma omp parallel + { +// initialize Data and Weight of empty bins in the first OpenMP thread to zero for later summation +# pragma omp for schedule( static ) collapse( 2 ) + for (int p=0; pNCell[b] = OMP_NCell[p][0][b]; + Ray[p]->Data [b] = ( Ray[p]->NCell[b] > 0L ) ? OMP_Data[p][0][b] / OMP_Weight[p][0][b] : 0.0; + } + } // OpenMP parallel region + + +// free per-thread arrays + Aux_DeallocateArray3D( OMP_Data ); + Aux_DeallocateArray3D( OMP_Weight ); + Aux_DeallocateArray3D( OMP_NCell ); + + delete [] Patch_Data; + delete [] Patch_BinMin_Rad; + delete [] Patch_BinMax_Rad; + delete [] Patch_BinMin_Tht; + delete [] Patch_BinMax_Tht; + delete [] Patch_BinMin_Phi; + delete [] Patch_BinMax_Phi; + + +// restore the original temporal interpolation setup + OPT__INT_TIME = IntTimeBackup; + +} // FUNCTION : Aux_ComputeRay + + + +//------------------------------------------------------------------------------------------------------- +// Function : Aux_GetMinMax_Vertex_Radius +// Description : Find the minimum and maximum radius among the cell vertices +// +// Parameter : x/y/z : Physical coordinates of cell center +// HalfWidth : Half cell width +// RadMin : Variable to store the output minimum radius +// RadMax : Variable to store the output maximum radius +// +// Return : RadMin, RadMax +//------------------------------------------------------------------------------------------------------- +void Aux_GetMinMax_Vertex_Radius( const double x, const double y, const double z, const double HalfWidth, + double *RadMin, double *RadMax ) +{ + + *RadMin = HUGE_NUMBER; + *RadMax = -HUGE_NUMBER; + + for (int k=-1; k<=1; k+=2) { const double zz = z + k * HalfWidth; const double zz2 = SQR(zz); + for (int j=-1; j<=1; j+=2) { const double yy = y + j * HalfWidth; const double yy2 = SQR(yy); + for (int i=-1; i<=1; i+=2) { const double xx = x + i * HalfWidth; const double xx2 = SQR(xx); + + const double Rad = sqrt( xx2 + yy2 + zz2 ); + + *RadMin = fmin( *RadMin, Rad ); + *RadMax = fmax( *RadMax, Rad ); + + }}} + +} // FUNCTION : Aux_GetMinMax_Vertex_Radius + + + +//------------------------------------------------------------------------------------------------------- +// Function : Aux_GetMinMax_Vertex_Angle +// Description : Find the minimum and maximum azimuthal angle and polar angle among the cell vertices +// +// Parameter : x/y/z : Physical coordinates of cell center +// HalfWidth : Half cell width +// ThtMin : Variable to store the output minimum polar angle +// ThtMax : Variable to store the output maximum polar angle +// PhiMin : Variable to store the output minimum azimuthal angle +// PhiMax : Variable to store the output maximum azimuthal angle +// +// Return : ThtMin, ThtMax, PhiMin, PhiMax +//------------------------------------------------------------------------------------------------------- +void Aux_GetMinMax_Vertex_Angle( const double x, const double y, const double z, const double HalfWidth, + double *ThtMin, double *ThtMax, double *PhiMin, double *PhiMax ) +{ + + *ThtMin = HUGE_NUMBER; + *ThtMax = -HUGE_NUMBER; + *PhiMin = HUGE_NUMBER; + *PhiMax = -HUGE_NUMBER; + + for (int k=-1; k<=1; k+=2) { const double zz = z + k * HalfWidth; const double zz2 = SQR(zz); + for (int j=-1; j<=1; j+=2) { const double yy = y + j * HalfWidth; const double yy2 = SQR(yy); + for (int i=-1; i<=1; i+=2) { const double xx = x + i * HalfWidth; const double xx2 = SQR(xx); + + const double Rad = sqrt( xx2 + yy2 + zz2 ); + +// exclude the corner case where the vertex lies on the reference z-axis + if ( xx != 0.0 || yy != 0.0 ) + { + double Phi = ( yy >= 0.0 ) + ? atan2( yy, xx ) + : atan2( yy, xx ) + TWOPI; + + if ( ( Phi == 0.0 ) && ( y < 0.0 ) ) Phi += TWOPI; + + *PhiMin = fmin( *PhiMin, Phi ); + *PhiMax = fmax( *PhiMax, Phi ); + } + + +// exclude the corner case where the vertex coincides with the reference center + if ( Rad != 0.0 ) + { + const double Tht = acos( zz / Rad ); + + *ThtMin = fmin( *ThtMin, Tht ); + *ThtMax = fmax( *ThtMax, Tht ); + } + }}} // i,j,k + +} // FUNCTION : Aux_GetMinMax_Vertex_Angle diff --git a/src/Auxiliary/Aux_TakeNote.cpp b/src/Auxiliary/Aux_TakeNote.cpp index 99c637f3f..f9dc8eeab 100644 --- a/src/Auxiliary/Aux_TakeNote.cpp +++ b/src/Auxiliary/Aux_TakeNote.cpp @@ -239,6 +239,8 @@ void Aux_TakeNote() fprintf( Note, "NEUTRINO_SCHEME IDSA\n" ); # elif ( NEUTRINO_SCHEME == M1 ) fprintf( Note, "NEUTRINO_SCHEME M1\n" ); +# elif ( NEUTRINO_SCHEME == LEAKAGE ) + fprintf( Note, "NEUTRINO_SCHEME LEAKAGE\n" ); # else fprintf( Note, "NEUTRINO_SCHEME UNKNOWN\n" ); # endif @@ -806,6 +808,7 @@ void Aux_TakeNote() # if ( MODEL == HYDRO ) fprintf( Note, "#define SRC_NAUX_DLEP % d\n", SRC_NAUX_DLEP ); fprintf( Note, "#define SRC_NAUX_LIGHTBULB % d\n", SRC_NAUX_LIGHTBULB ); + fprintf( Note, "#define SRC_NAUX_LEAKAGE % d\n", SRC_NAUX_LEAKAGE ); # endif fprintf( Note, "#define SRC_NAUX_USER % d\n", SRC_NAUX_USER ); # ifdef GPU @@ -1147,21 +1150,35 @@ void Aux_TakeNote() // record the parameters of source terms fprintf( Note, "Parameters of Source Terms\n" ); fprintf( Note, "***********************************************************************************\n" ); - fprintf( Note, "SRC_ANY % d\n", SrcTerms.Any ); - fprintf( Note, "SRC_DELEPTONIZATION % d\n", SrcTerms.Deleptonization ); - fprintf( Note, "SRC_LIGHTBULB % d\n", SrcTerms.Lightbulb ); + fprintf( Note, "SRC_ANY % d\n", SrcTerms.Any ); + fprintf( Note, "SRC_DELEPTONIZATION % d\n", SrcTerms.Deleptonization ); if ( SrcTerms.Deleptonization ) { - fprintf( Note, "SRC_DELEP_ENU % 14.7e\n", SrcTerms.Dlep_Enu ); - fprintf( Note, "SRC_DELEP_RHO1 % 14.7e\n", SrcTerms.Dlep_Rho1 ); - fprintf( Note, "SRC_DELEP_RHO2 % 14.7e\n", SrcTerms.Dlep_Rho2 ); - fprintf( Note, "SRC_DELEP_YE1 % 14.7e\n", SrcTerms.Dlep_Ye1 ); - fprintf( Note, "SRC_DELEP_YE2 % 14.7e\n", SrcTerms.Dlep_Ye2 ); - fprintf( Note, "SRC_DELEP_YEC % 14.7e\n", SrcTerms.Dlep_Yec ); } - if ( SrcTerms.Lightbulb ) { - fprintf( Note, "SRC_LIGHTBULB_LNUE % 14.7e\n", SrcTerms.Lightbulb_Lnue ); - fprintf( Note, "SRC_LIGHTBULB_TNUE % 14.7e\n", SrcTerms.Lightbulb_Tnue ); } - fprintf( Note, "SRC_USER % d\n", SrcTerms.User ); - fprintf( Note, "SRC_GPU_NPGROUP % d\n", SRC_GPU_NPGROUP ); + fprintf( Note, "SRC_DELEP_ENU % 14.7e\n", SrcTerms.Dlep_Enu ); + fprintf( Note, "SRC_DELEP_RHO1 % 14.7e\n", SrcTerms.Dlep_Rho1 ); + fprintf( Note, "SRC_DELEP_RHO2 % 14.7e\n", SrcTerms.Dlep_Rho2 ); + fprintf( Note, "SRC_DELEP_YE1 % 14.7e\n", SrcTerms.Dlep_Ye1 ); + fprintf( Note, "SRC_DELEP_YE2 % 14.7e\n", SrcTerms.Dlep_Ye2 ); + fprintf( Note, "SRC_DELEP_YEC % 14.7e\n", SrcTerms.Dlep_Yec ); } +# ifdef NEUTRINO_SCHEME +# if ( NEUTRINO_SCHEME == LIGHTBULB ) + fprintf( Note, "SRC_LIGHTBULB % d\n", SrcTerms.Lightbulb ); + fprintf( Note, "SRC_LIGHTBULB_LNUE % 14.7e\n", SrcTerms.Lightbulb_Lnue ); + fprintf( Note, "SRC_LIGHTBULB_TNUE % 14.7e\n", SrcTerms.Lightbulb_Tnue ); +# elif ( NEUTRINO_SCHEME == LEAKAGE ) + fprintf( Note, "SRC_LEAKAGE % d\n", SrcTerms.Leakage ); + fprintf( Note, "SRC_LEAKAGE_NRADIUS % d\n", SrcTerms.Leakage_NRadius ); + fprintf( Note, "SRC_LEAKAGE_NTHETA % d\n", SrcTerms.Leakage_NTheta ); + fprintf( Note, "SRC_LEAKAGE_NPHI % d\n", SrcTerms.Leakage_NPhi ); + fprintf( Note, "SRC_LEAKAGE_BINSIZE_RADIUS % 14.7e\n", SrcTerms.Leakage_BinSize_Radius ); + fprintf( Note, "SRC_LEAKAGE_RADIUSMAX % 14.7e\n", SrcTerms.Leakage_RadiusMax ); + fprintf( Note, "SRC_LEAKAGE_RADIUSMIN_LOG % 14.7e\n", SrcTerms.Leakage_RadiusMin_Log ); + fprintf( Note, "SRC_LEAKAGE_NUHEAT % d\n", SrcTerms.Leakage_NuHeat ); + fprintf( Note, "SRC_LEAKAGE_NUHEAT_FAC % 14.7e\n", SrcTerms.Leakage_NuHeat_Fac ); + fprintf( Note, "SRC_LEAKAGE_OPT_TEMP % d\n", SrcTerms.Leakage_Opt_Temp ); +# endif +# endif // #ifdef NEUTRINO_SCHEME + fprintf( Note, "SRC_USER % d\n", SrcTerms.User ); + fprintf( Note, "SRC_GPU_NPGROUP % d\n", SRC_GPU_NPGROUP ); fprintf( Note, "***********************************************************************************\n" ); fprintf( Note, "\n\n" ); diff --git a/src/GPU_API/CUAPI_MemAllocate_Fluid.cu b/src/GPU_API/CUAPI_MemAllocate_Fluid.cu index 43cbf709b..4e4698a34 100644 --- a/src/GPU_API/CUAPI_MemAllocate_Fluid.cu +++ b/src/GPU_API/CUAPI_MemAllocate_Fluid.cu @@ -43,6 +43,14 @@ extern real (*d_FC_Mag_Half)[NCOMP_MAG][ FLU_NXT_P1*SQR(FLU_NXT) ]; extern real (*d_EC_Ele )[NCOMP_MAG][ CUBE(N_EC_ELE) ]; #endif #endif // FLU_SCHEME +#if ( MODEL == HYDRO ) +extern real *d_SrcLeakage_Radius; +extern real *d_SrcLeakage_tau; +extern real *d_SrcLeakage_chi; +extern real *d_SrcLeakage_HeatFlux; +extern real *d_SrcLeakage_HeatERms; +extern real *d_SrcLeakage_HeatEAve; +#endif #if ( MODEL == ELBDM ) extern bool (*d_IsCompletelyRefined); @@ -110,6 +118,14 @@ int CUAPI_MemAllocate_Fluid( const int Flu_NPG, const int Pot_NPG, const int Src const long Flu_MemSize_S_In = sizeof(real )*Src_NP*FLU_NIN_S *CUBE(SRC_NXT); const long Flu_MemSize_S_Out = sizeof(real )*Src_NP*FLU_NOUT_S*CUBE(PS1); const long Corner_MemSize_S = sizeof(double)*Src_NP*3; +# if ( MODEL == HYDRO && NEUTRINO_SCHEME == LEAKAGE ) + const long Leak_Rad_MemSize = sizeof(real )* SrcTerms.Leakage_NRadius; + const long Leak_Tau_MemSize = sizeof(real )*(SrcTerms.Leakage_NRadius*SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3); + const long Leak_Chi_MemSize = sizeof(real )*(SrcTerms.Leakage_NRadius*SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3); + const long Leak_Flux_MemSize = sizeof(real )*(SrcTerms.Leakage_NRadius*SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3); + const long Leak_ERms_MemSize = sizeof(real )*( SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3); + const long Leak_EAve_MemSize = sizeof(real )*( SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3); +# endif // the size of the global memory arrays in different models # if ( FLU_SCHEME == MHM || FLU_SCHEME == MHM_RP || FLU_SCHEME == CTU ) @@ -203,6 +219,11 @@ int CUAPI_MemAllocate_Fluid( const int Flu_NPG, const int Pot_NPG, const int Src TotalSize += Corner_MemSize_S; } +# if ( MODEL == HYDRO && NEUTRINO_SCHEME == LEAKAGE ) + TotalSize += Leak_Rad_MemSize + Leak_Tau_MemSize + Leak_Chi_MemSize + + Leak_Flux_MemSize + Leak_ERms_MemSize + Leak_EAve_MemSize; +# endif + if ( MPI_Rank == 0 ) Aux_Message( stdout, "NOTE : total memory requirement in GPU fluid solver = %ld MB\n", TotalSize/(1<<20) ); @@ -263,6 +284,23 @@ int CUAPI_MemAllocate_Fluid( const int Flu_NPG, const int Pot_NPG, const int Src CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_Corner_Array_S, Corner_MemSize_S ) ); } +# if ( MODEL == HYDRO && NEUTRINO_SCHEME == LEAKAGE ) + CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_SrcLeakage_Radius, Leak_Rad_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_SrcLeakage_tau, Leak_Tau_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_SrcLeakage_chi, Leak_Chi_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_SrcLeakage_HeatFlux, Leak_Flux_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_SrcLeakage_HeatERms, Leak_ERms_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_SrcLeakage_HeatEAve, Leak_EAve_MemSize ) ); + +// store the device pointers in SrcTerms when using GPU + SrcTerms.Leakage_Radius_DevPtr = d_SrcLeakage_Radius; + SrcTerms.Leakage_tau_DevPtr = d_SrcLeakage_tau; + SrcTerms.Leakage_chi_DevPtr = d_SrcLeakage_chi; + SrcTerms.Leakage_HeatFlux_DevPtr = d_SrcLeakage_HeatFlux; + SrcTerms.Leakage_HeatERms_DevPtr = d_SrcLeakage_HeatERms; + SrcTerms.Leakage_HeatEAve_DevPtr = d_SrcLeakage_HeatEAve; +# endif + # if ( MODEL == ELBDM ) CUDA_CHECK_MALLOC( cudaMalloc( (void**) &d_IsCompletelyRefined, Flu_MemSize_IsCompletelyRefined ) ); @@ -337,6 +375,15 @@ int CUAPI_MemAllocate_Fluid( const int Flu_NPG, const int Pot_NPG, const int Src CUDA_CHECK_MALLOC( cudaMallocHost( (void**) &h_GramFE_TimeEvo, GramFE_TimeEvo_MemSize ) ); # endif +# if ( MODEL == HYDRO && NEUTRINO_SCHEME == LEAKAGE ) + CUDA_CHECK_MALLOC( cudaMallocHost( (void**) &h_SrcLeakage_Radius, Leak_Rad_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMallocHost( (void**) &h_SrcLeakage_tau, Leak_Tau_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMallocHost( (void**) &h_SrcLeakage_chi, Leak_Chi_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMallocHost( (void**) &h_SrcLeakage_HeatFlux, Leak_Flux_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMallocHost( (void**) &h_SrcLeakage_HeatERms, Leak_ERms_MemSize ) ); + CUDA_CHECK_MALLOC( cudaMallocHost( (void**) &h_SrcLeakage_HeatEAve, Leak_EAve_MemSize ) ); +# endif + // create streams Stream = new cudaStream_t [GPU_NStream]; for (int s=0; s as we still rely on these constants (e.g., DENS, DUAL) in the fluid solvers # if ( EOS == EOS_NUCLEAR ) # if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) - Idx_Temp_IG = AddField( "Temp_IG", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); - if ( Idx_Temp_IG != TEMP_IG ) Aux_Error( ERROR_INFO, "inconsistent Idx_Temp_IG (%d != %d) !!\n", Idx_Temp_IG, TEMP_IG ); + Idx_Temp_IG = AddField( "Temp_IG", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); + if ( Idx_Temp_IG != TEMP_IG ) Aux_Error( ERROR_INFO, "inconsistent Idx_Temp_IG (%d != %d) !!\n", Idx_Temp_IG, TEMP_IG ); # endif - Idx_dEdt_Nu = AddField( "dEdt_Nu", FIXUP_FLUX_NO, FIXUP_REST_NO, FLOOR_NO, NORMALIZE_NO, INTERP_FRAC_NO ); - if ( Idx_dEdt_Nu != DEDT_NU ) Aux_Error( ERROR_INFO, "inconsistent Idx_dEdt_Nu (%d != %d) !!\n", Idx_dEdt_Nu, DEDT_NU ); + Idx_dEdt_Nu = AddField( "dEdt_Nu", FIXUP_FLUX_NO, FIXUP_REST_NO, FLOOR_NO, NORMALIZE_NO, INTERP_FRAC_NO ); + if ( Idx_dEdt_Nu != DEDT_NU ) Aux_Error( ERROR_INFO, "inconsistent Idx_dEdt_Nu (%d != %d) !!\n", Idx_dEdt_Nu, DEDT_NU ); - Idx_Ye = AddField( "Ye" , FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_YES ); - if ( Idx_Ye != YE ) Aux_Error( ERROR_INFO, "inconsistent Idx_Ye (%d != %d) !!\n", Idx_Ye, YE ); +# if ( NEUTRINO_SCHEME == LEAKAGE ) + Idx_dYedt_Nu = AddField( "dYedt_Nu", FIXUP_FLUX_NO, FIXUP_REST_NO, FLOOR_NO, NORMALIZE_NO, INTERP_FRAC_NO ); + if ( Idx_dYedt_Nu != DYEDT_NU ) Aux_Error( ERROR_INFO, "inconsistent Idx_dYedt_Nu (%d != %d) !!\n", Idx_dYedt_Nu, DYEDT_NU ); # endif + Idx_Ye = AddField( "Ye", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_YES ); + if ( Idx_Ye != YE ) Aux_Error( ERROR_INFO, "inconsistent Idx_Ye (%d != %d) !!\n", Idx_Ye, YE ); +# endif // # if ( EOS == EOS_NUCLEAR ) + # ifdef COSMIC_RAY Idx_CRay = AddField( "CRay", FIXUP_FLUX_YES, FIXUP_REST_YES, FLOOR_YES, NORMALIZE_NO, INTERP_FRAC_NO ); if ( Idx_CRay != CRAY ) Aux_Error( ERROR_INFO, "inconsistent Idx_CRay (%d != %d) !!\n", Idx_CRay, CRAY ); diff --git a/src/Init/Init_GAMER.cpp b/src/Init/Init_GAMER.cpp index c4974cdd7..39b2fbe61 100644 --- a/src/Init/Init_GAMER.cpp +++ b/src/Init/Init_GAMER.cpp @@ -354,6 +354,12 @@ void Init_GAMER( int *argc, char ***argv ) if ( MPI_Rank == 0 ) Aux_Message( stdout, "%s ...\n", "Initializing source-term fields" ); +// disable source-term modules that take effects on active fields even dt=0.0 + const bool Backup_Deleptonization = SrcTerms.Deleptonization; + + SrcTerms.Deleptonization = false; + + for (int lv=0; lvAdd( "SRC_LIGHTBULB", &SrcTerms.Lightbulb, false, Useless_bool, Useless_bool ); ReadPara->Add( "SRC_LIGHTBULB_LNUE", &SrcTerms.Lightbulb_Lnue, 1.0e52, 0.0, NoMax_double ); ReadPara->Add( "SRC_LIGHTBULB_TNUE", &SrcTerms.Lightbulb_Tnue, 4.0, 0.0, NoMax_double ); + ReadPara->Add( "SRC_LEAKAGE", &SrcTerms.Leakage, false, Useless_bool, Useless_bool ); + ReadPara->Add( "SRC_LEAKAGE_NRADIUS", &SrcTerms.Leakage_NRadius, 768, 1, NoMax_int ); + ReadPara->Add( "SRC_LEAKAGE_NTHETA", &SrcTerms.Leakage_NTheta, 1, 1, NoMax_int ); + ReadPara->Add( "SRC_LEAKAGE_NPHI", &SrcTerms.Leakage_NPhi, 1, 1, NoMax_int ); + ReadPara->Add( "SRC_LEAKAGE_BINSIZE_RADIUS", &SrcTerms.Leakage_BinSize_Radius, 1.0e5, Eps_double, NoMax_double ); + ReadPara->Add( "SRC_LEAKAGE_RADIUSMAX", &SrcTerms.Leakage_RadiusMax, 3.0e8, Eps_double, NoMax_double ); + ReadPara->Add( "SRC_LEAKAGE_RADIUSMIN_LOG", &SrcTerms.Leakage_RadiusMin_Log, 3.0e7, Eps_double, NoMax_double ); + ReadPara->Add( "SRC_LEAKAGE_NUHEAT", &SrcTerms.Leakage_NuHeat, true, Useless_bool, Useless_bool ); + ReadPara->Add( "SRC_LEAKAGE_NUHEAT_FAC", &SrcTerms.Leakage_NuHeat_Fac, 1.0, 0.0, NoMax_double ); + ReadPara->Add( "SRC_LEAKAGE_OPT_TEMP", &SrcTerms.Leakage_Opt_Temp, LEAK_TEMP_BINDATA, LEAK_TEMP_INDCELL, LEAK_TEMP_BINDATA ); ReadPara->Add( "SRC_USER", &SrcTerms.User, false, Useless_bool, Useless_bool ); // do not check SRC_GPU_NPGROUP since it may be reset by either Init_ResetParameter() or CUAPI_SetMemSize() ReadPara->Add( "SRC_GPU_NPGROUP", &SRC_GPU_NPGROUP, -1, NoMin_int, NoMax_int ); diff --git a/src/Init/Init_MemAllocate_Fluid.cpp b/src/Init/Init_MemAllocate_Fluid.cpp index c629151dd..de397d6a6 100644 --- a/src/Init/Init_MemAllocate_Fluid.cpp +++ b/src/Init/Init_MemAllocate_Fluid.cpp @@ -103,6 +103,22 @@ void Init_MemAllocate_Fluid( const int Flu_NPatchGroup, const int Pot_NPatchGrou h_GramFE_TimeEvo = new gramfe_matmul_float [PS2][ 2*FLU_NXT ]; # endif +# if ( MODEL == HYDRO && NEUTRINO_SCHEME == LEAKAGE ) + h_SrcLeakage_Radius = new real [SrcTerms.Leakage_NRadius]; + h_SrcLeakage_tau = new real [SrcTerms.Leakage_NRadius*SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3]; + h_SrcLeakage_chi = new real [SrcTerms.Leakage_NRadius*SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3]; + h_SrcLeakage_HeatFlux = new real [SrcTerms.Leakage_NRadius*SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3]; + h_SrcLeakage_HeatERms = new real [ SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3]; + h_SrcLeakage_HeatEAve = new real [ SrcTerms.Leakage_NTheta*SrcTerms.Leakage_NPhi*3]; + + SrcTerms.Leakage_Radius_DevPtr = h_SrcLeakage_Radius; + SrcTerms.Leakage_tau_DevPtr = h_SrcLeakage_tau; + SrcTerms.Leakage_chi_DevPtr = h_SrcLeakage_chi; + SrcTerms.Leakage_HeatFlux_DevPtr = h_SrcLeakage_HeatFlux; + SrcTerms.Leakage_HeatERms_DevPtr = h_SrcLeakage_HeatERms; + SrcTerms.Leakage_HeatEAve_DevPtr = h_SrcLeakage_HeatEAve; +# endif + } // FUNCTION : Init_MemAllocate_Fluid diff --git a/src/Init/Init_ResetParameter.cpp b/src/Init/Init_ResetParameter.cpp index 5e3bd9b97..acbf77866 100644 --- a/src/Init/Init_ResetParameter.cpp +++ b/src/Init/Init_ResetParameter.cpp @@ -1250,6 +1250,24 @@ void Init_ResetParameter() # endif // #ifdef STAR_FORMATION +# if ( defined NEUTRINO_SCHEME && NEUTRINO_SCHEME == LEAKAGE ) + SrcTerms.Leakage_BinSize_Radius /= UNIT_L; + SrcTerms.Leakage_RadiusMax /= UNIT_L; + SrcTerms.Leakage_RadiusMin_Log /= UNIT_L; + + PRINT_RESET_PARA( SrcTerms.Leakage_BinSize_Radius, FORMAT_REAL, "to be consistent with the code units" ); + PRINT_RESET_PARA( SrcTerms.Leakage_RadiusMax, FORMAT_REAL, "to be consistent with the code units" ); + PRINT_RESET_PARA( SrcTerms.Leakage_RadiusMin_Log, FORMAT_REAL, "to be consistent with the code units" ); + + if ( Mis_CompareRealValue( SrcTerms.Leakage_RadiusMax, SrcTerms.Leakage_RadiusMin_Log, NULL, false ) ) + { + SrcTerms.Leakage_NRadius = int( SrcTerms.Leakage_RadiusMin_Log / SrcTerms.Leakage_BinSize_Radius ); + + PRINT_RESET_PARA( SrcTerms.Leakage_NRadius, FORMAT_INT, "to be consistent with the specified maximum radius" ); + } +# endif // if ( defined NEUTRINO_SCHEME && NEUTRINO_SCHEME == LEAKAGE ) + + // disable OPT__MINIMIZE_MPI_BARRIER in the serial mode # ifdef SERIAL if ( OPT__MINIMIZE_MPI_BARRIER ) @@ -1322,14 +1340,14 @@ void Init_ResetParameter() { GREP_MINBINSIZE = amr->dh[MAX_LEVEL]; - PRINT_RESET_PARA( GREP_MINBINSIZE, FORMAT_FLT, "" ); + PRINT_RESET_PARA( GREP_MINBINSIZE, FORMAT_REAL, "" ); } if ( GREP_MAXRADIUS <= 0.0 ) { GREP_MAXRADIUS = sqrt(3.0) * amr->BoxSize[0]; - PRINT_RESET_PARA( GREP_MAXRADIUS, FORMAT_FLT, "" ); + PRINT_RESET_PARA( GREP_MAXRADIUS, FORMAT_REAL, "" ); } # endif diff --git a/src/Main/Main.cpp b/src/Main/Main.cpp index 82902ff1a..861406b6a 100644 --- a/src/Main/Main.cpp +++ b/src/Main/Main.cpp @@ -340,6 +340,8 @@ double Src_Dlep_AuxArray_Flt [SRC_NAUX_DLEP ]; int Src_Dlep_AuxArray_Int [SRC_NAUX_DLEP ]; double Src_Lightbulb_AuxArray_Flt[SRC_NAUX_LIGHTBULB]; int Src_Lightbulb_AuxArray_Int[SRC_NAUX_LIGHTBULB]; +double Src_Leakage_AuxArray_Flt [SRC_NAUX_LEAKAGE ]; +int Src_Leakage_AuxArray_Int [SRC_NAUX_LEAKAGE ]; #endif double Src_User_AuxArray_Flt [SRC_NAUX_USER ]; int Src_User_AuxArray_Int [SRC_NAUX_USER ]; @@ -477,6 +479,14 @@ real (*h_Flu_Array_S_Out[2])[FLU_NOUT_S][ CUBE(PS1) ] = { NULL, NUL real (*h_Mag_Array_S_In [2])[NCOMP_MAG][ SRC_NXT_P1*SQR(SRC_NXT) ] = { NULL, NULL }; #endif double (*h_Corner_Array_S[2])[3] = { NULL, NULL }; +#if ( MODEL == HYDRO ) +real *h_SrcLeakage_Radius = NULL; +real *h_SrcLeakage_tau = NULL; +real *h_SrcLeakage_chi = NULL; +real *h_SrcLeakage_HeatFlux = NULL; +real *h_SrcLeakage_HeatERms = NULL; +real *h_SrcLeakage_HeatEAve = NULL; +#endif @@ -564,6 +574,14 @@ real (*d_Flu_Array_S_Out)[FLU_NOUT_S][ CUBE(PS1) ] = NULL; real (*d_Mag_Array_S_In)[NCOMP_MAG ][ SRC_NXT_P1*SQR(SRC_NXT) ] = NULL; #endif double (*d_Corner_Array_S)[3] = NULL; +#if ( MODEL == HYDRO ) +real *d_SrcLeakage_Radius = NULL; +real *d_SrcLeakage_tau = NULL; +real *d_SrcLeakage_chi = NULL; +real *d_SrcLeakage_HeatFlux = NULL; +real *d_SrcLeakage_HeatERms = NULL; +real *d_SrcLeakage_HeatEAve = NULL; +#endif #endif // #ifdef GPU diff --git a/src/Makefile_base b/src/Makefile_base index a8ff7da00..a5be58982 100644 --- a/src/Makefile_base +++ b/src/Makefile_base @@ -78,7 +78,8 @@ CPU_FILE += Aux_Check_Parameter.cpp Aux_Check_Conservation.cpp Aux_Check.cp Aux_GetMemInfo.cpp Aux_Message.cpp Aux_Record_PatchCount.cpp Aux_TakeNote.cpp Aux_Timing.cpp \ Aux_Check_MemFree.cpp Aux_Record_Performance.cpp Aux_CheckFileExist.cpp Aux_Array.cpp \ Aux_Record_User.cpp Aux_Record_CorrUnphy.cpp Aux_Record_Center.cpp Aux_SwapPointer.cpp Aux_Check_NormalizePassive.cpp \ - Aux_LoadTable.cpp Aux_IsFinite.cpp Aux_ComputeProfile.cpp Aux_FindExtrema.cpp Aux_FindWeightedAverageCenter.cpp Aux_PauseManually.cpp + Aux_LoadTable.cpp Aux_IsFinite.cpp Aux_ComputeProfile.cpp Aux_FindExtrema.cpp Aux_FindWeightedAverageCenter.cpp Aux_PauseManually.cpp \ + Aux_ComputeRay.cpp CPU_FILE += CPU_FluidSolver.cpp Flu_AdvanceDt.cpp Flu_Prepare.cpp Flu_Close.cpp Flu_FixUp_Flux.cpp \ Flu_FixUp_Restrict.cpp Flu_AllocateFluxArray.cpp Flu_BoundaryCondition_User.cpp Flu_ResetByUser.cpp \ @@ -311,16 +312,16 @@ endif # SUPPORT_LIBYT # local source terms source files # ------------------------------------------------------------------------------------ GPU_FILE += CUAPI_Asyn_SrcSolver.cu CUSRC_SrcSolver_IterateAllCells.cu CUSRC_Src_Deleptonization.cu \ - CUSRC_Src_User_Template.cu CUSRC_Src_Lightbulb.cu + CUSRC_Src_User_Template.cu CUSRC_Src_Lightbulb.cu CUSRC_Src_Leakage.cu CUSRC_Src_Leakage_ComputeLeak.cu CPU_FILE += CPU_SrcSolver.cpp CPU_SrcSolver_IterateAllCells.cpp CPU_Src_Deleptonization.cpp \ - CPU_Src_User_Template.cpp CPU_Src_Lightbulb.cpp + CPU_Src_User_Template.cpp CPU_Src_Lightbulb.cpp CPU_Src_Leakage.cpp CPU_Src_Leakage_ComputeLeak.cpp CPU_FILE += Src_AdvanceDt.cpp Src_Prepare.cpp Src_Close.cpp Src_Init.cpp Src_End.cpp \ Src_WorkBeforeMajorFunc.cpp -vpath %.cu SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization SourceTerms/LightBulb -vpath %.cpp SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization SourceTerms/LightBulb +vpath %.cu SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization SourceTerms/LightBulb SourceTerms/Leakage +vpath %.cpp SourceTerms SourceTerms/User_Template SourceTerms/Deleptonization SourceTerms/LightBulb SourceTerms/Leakage # Grackle source files (included only if "SUPPORT_GRACKLE" is turned on) diff --git a/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp b/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp index bac49f49c..b6ee22896 100644 --- a/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp +++ b/src/SourceTerms/CPU_SrcSolver_IterateAllCells.cpp @@ -117,9 +117,14 @@ void CPU_SrcSolver_IterateAllCells( if ( SrcTerms.Lightbulb ) SrcTerms.Lightbulb_FuncPtr( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, PassiveFloor, &EoS, SrcTerms.Lightbulb_AuxArrayDevPtr_Flt, SrcTerms.Lightbulb_AuxArrayDevPtr_Int ); -# endif -// (3) user-defined +// (3) leakage + if ( SrcTerms.Leakage ) + SrcTerms.Leakage_FuncPtr ( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, PassiveFloor, &EoS, + SrcTerms.Leakage_AuxArrayDevPtr_Flt, SrcTerms.Leakage_AuxArrayDevPtr_Int ); +# endif // if ( MODEL == HYDRO ) + +// (4) user-defined if ( SrcTerms.User ) SrcTerms.User_FuncPtr ( fluid, B, &SrcTerms, dt, dh, x, y, z, TimeNew, TimeOld, MinDens, MinPres, MinEint, PassiveFloor, &EoS, SrcTerms.User_AuxArrayDevPtr_Flt, SrcTerms.User_AuxArrayDevPtr_Int ); diff --git a/src/SourceTerms/Leakage/CPU_Src_Leakage.cpp b/src/SourceTerms/Leakage/CPU_Src_Leakage.cpp new file mode 100644 index 000000000..3fd973c03 --- /dev/null +++ b/src/SourceTerms/Leakage/CPU_Src_Leakage.cpp @@ -0,0 +1,1150 @@ +#include "CUFLU.h" +#include "Global.h" +#include "NuclearEoS.h" +#include "Src_Leakage.h" + +#if ( MODEL == HYDRO ) + + + +// external functions and GPU-related set-up +#ifdef __CUDACC__ + +#include "CUDA_CheckError.h" +#include "CUFLU_Shared_FluUtility.cu" +#include "CUDA_ConstMemory.h" +#include "CUSRC_Src_Leakage_ComputeLeak.cu" + +extern real *d_SrcLeakage_Radius; +extern real *d_SrcLeakage_tau; +extern real *d_SrcLeakage_chi; +extern real *d_SrcLeakage_HeatFlux; +extern real *d_SrcLeakage_HeatERms; +extern real *d_SrcLeakage_HeatEAve; + +#endif // #ifdef __CUDACC__ + + +// external variables, local and external function prototypes +#ifndef __CUDACC__ + +// construct the leakage profiles sequentially to minimize memory usage +// --> avoid bad_alloc errors when invoking Aux_AllocateArray3D() in Aux_ComputeRay() +//#define LEAKAGE_SEQUENTIAL_MAPPING + +static long Leakage_Step = -1; + +#if ( EOS == EOS_NUCLEAR ) +extern int g_nye; +#endif + +void Src_SetAuxArray_Leakage( double [], int [] ); +void Src_SetCPUFunc_Leakage( SrcFunc_t & ); +#ifdef GPU +void Src_SetGPUFunc_Leakage( SrcFunc_t & ); +#endif +void Src_SetConstMemory_Leakage( const double AuxArray_Flt[], const int AuxArray_Int[], + double *&DevPtr_Flt, int *&DevPtr_Int ); +void Src_PassData2GPU_Leakage(); +double Src_Leakage_ConstructSeries( const int NBin, const double xmin, const double xmax, const double dx ); + +extern void Src_Leakage_ComputeTau( Profile_t *Ray[], double *Edge, + const int NRadius, const int NTheta, const int NPhi, + real *tau_Ruff, real *chi_Ross, real *Heat_Flux, + real *Heat_ERms, real *Heat_Eave ); +extern void Src_Leakage_ComputeLeak( const real Dens_Code, const real Temp_Kelv, const real Ye, const real chi[], const real tau[], + const real Heat_Flux[], const real *Heat_ERms, const real *Heat_EAve, + real *dEdt, real *dYedt, real *Lum, real *Heat, real *NetHeat, + const bool NuHeat, const real NuHeat_Fac, const real UNIT_D, const EoS_t *EoS ); + +#endif // #ifndef __CUDACC__ + + +GPU_DEVICE static +int Src_Leakage_BinarySearch( const real Array[], int Min, int Max, const real Key ); + +GPU_DEVICE static +real Src_Leakage_LinearInterp( const real *array, const real *xs, const real x_in ); + +GPU_DEVICE static +real Src_Leakage_BilinearInterp( const real *array, const real *xs, const real *ys, + const real x_in, const real y_in ); + +GPU_DEVICE static +real Src_Leakage_TrilinearInterp( const real *array, const real *xs, const real *ys, const real *zs, + const real x_in, const real y_in, const real z_in ); + + + +/******************************************************** +1. Leakage source term + --> Enabled by the runtime option "SRC_LEAKAGE" + +2. This file is shared by both CPU and GPU + + CUSRC_Src_Leakage.cu -> CPU_Src_Leakage.cpp + +3. Four steps are required to implement a source term + + I. Set auxiliary arrays + II. Implement the source-term function + III. [Optional] Add the work to be done every time + before calling the major source-term function + IV. Set initialization functions + +4. The source-term function must be thread-safe and + not use any global variable +********************************************************/ + + + +// ======================= +// I. Set auxiliary arrays +// ======================= + +//------------------------------------------------------------------------------------------------------- +// Function : Src_SetAuxArray_Leakage +// Description : Set the auxiliary arrays AuxArray_Flt/Int[] +// +// Note : 1. Invoked by Src_Init_Leakage() +// 2. AuxArray_Flt/Int[] have the size of SRC_NAUX_LEAKAGE defined in Macro.h (default = 11) +// 3. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : AuxArray_Flt/Int : Floating-point/Integer arrays to be filled up +// +// Return : AuxArray_Flt/Int[] +//------------------------------------------------------------------------------------------------------- +#ifndef __CUDACC__ +void Src_SetAuxArray_Leakage( double AuxArray_Flt[], int AuxArray_Int[] ) +{ + + const int RayConfig = ( (SrcTerms.Leakage_NTheta == 1) << 1 ) + (SrcTerms.Leakage_NPhi == 1); + const int NRay = SrcTerms.Leakage_NTheta * SrcTerms.Leakage_NPhi; + const double Ye_Frac = 0.01; + + AuxArray_Flt[SRC_AUX_PNS_X ] = amr->BoxCenter[0]; + AuxArray_Flt[SRC_AUX_PNS_Y ] = amr->BoxCenter[1]; + AuxArray_Flt[SRC_AUX_PNS_Z ] = amr->BoxCenter[2]; + AuxArray_Flt[SRC_AUX_MAXRADIUS ] = SrcTerms.Leakage_RadiusMax; + AuxArray_Flt[SRC_AUX_RADMIN_LOG] = int( SrcTerms.Leakage_RadiusMin_Log / SrcTerms.Leakage_BinSize_Radius ) + * SrcTerms.Leakage_BinSize_Radius; + AuxArray_Flt[SRC_AUX_DRAD ] = SrcTerms.Leakage_BinSize_Radius; + AuxArray_Flt[SRC_AUX_DTHETA ] = M_PI / SrcTerms.Leakage_NTheta; + AuxArray_Flt[SRC_AUX_DPHI ] = 2.0 * M_PI / SrcTerms.Leakage_NPhi; +# if ( EOS == EOS_NUCLEAR ) + AuxArray_Flt[SRC_AUX_YEMIN ] = ( 1.0 + Ye_Frac ) * h_EoS_Table[NUC_TAB_YE][0]; + AuxArray_Flt[SRC_AUX_YEMAX ] = ( 1.0 - Ye_Frac ) * h_EoS_Table[NUC_TAB_YE][g_nye-1]; +# endif + AuxArray_Flt[SRC_AUX_VSQR2CODE ] = 1.0 / SQR( UNIT_V ); + + AuxArray_Int[SRC_AUX_RAYCONF ] = RayConfig; + AuxArray_Int[SRC_AUX_NRAD_LIN ] = int( SrcTerms.Leakage_RadiusMin_Log / SrcTerms.Leakage_BinSize_Radius ); + AuxArray_Int[SRC_AUX_STRIDE ] = NRay * NType_Neutrino; + AuxArray_Int[SRC_AUX_MODE ] = LEAK_MODE_EVOLVE; + +} // FUNCTION : Src_SetAuxArray_Leakage +#endif // #ifndef __CUDACC__ + + + +// ====================================== +// II. Implement the source-term function +// ====================================== + +//------------------------------------------------------------------------------------------------------- +// Function : Src_Leakage +// Description : Apply the leakage scheme for neutrino heating and cooling, and the change in electron fraction +// +// Note : 1. Invoked by CPU/GPU_SrcSolver_IterateAllCells() +// 2. See Src_SetAuxArray_Leakage() for the values stored in AuxArray_Flt/Int[] +// 3. Shared by both CPU and GPU +// +// Parameter : fluid : Fluid array storing both the input and updated values +// --> Including both active and passive variables +// B : Cell-centered magnetic field +// SrcTerms : Structure storing all source-term variables +// dt : Time interval to advance solution +// dh : Grid size +// x/y/z : Target physical coordinates +// TimeNew : Target physical time to reach +// TimeOld : Physical time before update +// --> This function updates physical time from TimeOld to TimeNew +// MinDens/Pres/Eint : Density, pressure, and internal energy floors +// EoS : EoS object +// AuxArray_* : Auxiliary arrays (see the Note above) +// +// Return : fluid[] +//----------------------------------------------------------------------------------------- +GPU_DEVICE_NOINLINE +static void Src_Leakage( real fluid[], const real B[], + const SrcTerms_t *SrcTerms, const real dt, const real dh, + const double x, const double y, const double z, + const double TimeNew, const double TimeOld, + const real MinDens, const real MinPres, const real MinEint, const long PassiveFloor, + const EoS_t *EoS, const double AuxArray_Flt[], const int AuxArray_Int[] ) +{ + +// check +# ifdef GAMER_DEBUG + if ( AuxArray_Flt == NULL ) printf( "ERROR : AuxArray_Flt == NULL in %s !!\n", __FUNCTION__ ); + if ( AuxArray_Int == NULL ) printf( "ERROR : AuxArray_Int == NULL in %s !!\n", __FUNCTION__ ); +# endif + + + const real PNS_x = AuxArray_Flt[SRC_AUX_PNS_X ]; + const real PNS_y = AuxArray_Flt[SRC_AUX_PNS_Y ]; + const real PNS_z = AuxArray_Flt[SRC_AUX_PNS_Z ]; + const real MaxRadius = AuxArray_Flt[SRC_AUX_MAXRADIUS ]; + const real RadMin_Log = AuxArray_Flt[SRC_AUX_RADMIN_LOG]; + const real dRad = AuxArray_Flt[SRC_AUX_DRAD ]; + const real dTht = AuxArray_Flt[SRC_AUX_DTHETA ]; + const real dPhi = AuxArray_Flt[SRC_AUX_DPHI ]; + const real YeMin = AuxArray_Flt[SRC_AUX_YEMIN ]; + const real YeMax = AuxArray_Flt[SRC_AUX_YEMAX ]; + const real sEint2Code = AuxArray_Flt[SRC_AUX_VSQR2CODE ]; + + const int RayConfig = AuxArray_Int[SRC_AUX_RAYCONF ]; + const int NRad_Lin = AuxArray_Int[SRC_AUX_NRAD_LIN ]; + const int Stride = AuxArray_Int[SRC_AUX_STRIDE ]; + const int Mode = AuxArray_Int[SRC_AUX_MODE ]; + + + const int NRadius = SrcTerms->Leakage_NRadius; + const int NTheta = SrcTerms->Leakage_NTheta; + const int NPhi = SrcTerms->Leakage_NPhi; +# ifdef __CUDACC__ + const real *Radius = SrcTerms->Leakage_Radius_DevPtr; +# else + const real *Radius = h_SrcLeakage_Radius; +# endif + + const double x0 = x - PNS_x; + const double y0 = y - PNS_y; + const double z0 = z - PNS_z; + const double Rad = sqrt( SQR(x0) + SQR(y0) + SQR(z0) ); + double Tht, Phi; + real xs[2], ys[2], zs[2]; + int Idx_Rad, Idx_Tht, Idx_Phi; + int NPhi_half = NPhi>>1; + + +// do nothing if the cell is beyond the sampled rays + if ( Rad > MaxRadius ) + { + if ( Mode == LEAK_MODE_RECORD ) + { + fluid[DENS] = (real)0.0; + fluid[MOMX] = (real)0.0; + fluid[MOMY] = (real)0.0; + fluid[MOMZ] = (real)0.0; + fluid[ENGY] = (real)0.0; +# ifdef YE + fluid[YE ] = (real)0.0; +# endif + } + +# ifdef DYEDT_NU + fluid[DEDT_NU ] = (real)0.0; + fluid[DYEDT_NU] = (real)0.0; +# endif + + return; + } + + +// (1) set up the data index and coordinate +// (1-1) spherical radius: +// --> use data at boundary for cells outside of the radius range +// --> for Rad <= RadMin_Log, add MIN( ..., NRadius-2 ) to deal with the special case RadMin_Log = MaxRadius + Idx_Rad = ( Rad <= RadMin_Log ) + ? MIN( MAX( int( Rad / dRad - (real)0.5 ), 0 ), NRadius-2 ) + : Src_Leakage_BinarySearch( Radius, NRad_Lin-1, NRadius-2, Rad ); + + xs[0] = Radius[Idx_Rad ]; + xs[1] = Radius[Idx_Rad+1]; + +// (1-2) polar angle, theta: +// --> for cells close to the pole (Tht < 0.5 * dTht or Tht > PI - 0.5 * dTht) +// (a) NPhi > 1, use data at Phi + PI for interpolation +// --> Idx_Tht = -1 if Tht < 0.5 * dTht +// NTheta - 1 if Tht > PI - 0.5 * dTht +// int( Tht / dTht - 0.5 ) otherwise +// (b) NPhi = 1, use data at boundary +// --> Idx_Tht = 0 if Tht < 0.5 * dTht +// NTheta - 2 if Tht > PI - 0.5 * dTht +// int( Tht / dTht - 0.5 ) otherwise + if ( NTheta > 1 ) + { +// set theta = 0.5 * PI arbitrarily if the cell is at the reference center + Tht = ( Rad == 0.0 ) + ? 0.5 * M_PI + : acos( z0 / Rad ); + Idx_Tht = ( NPhi > 1 ) + ? int( Tht / dTht + 0.5 ) - 1 + : MIN( MAX( int( Tht / dTht - 0.5 ), 0 ), NTheta-2 ); + + ys[0] = ( (real)0.5 + (real)Idx_Tht ) * dTht; + ys[1] = ( (real)1.5 + (real)Idx_Tht ) * dTht; + } + +// (1-3) azimuthal angle, phi: +// --> use periodic boundary condition +// --> Idx_Phi = NPhi - 1 if Phi < 0.5 * dPhi or Phi > 2 * PI - 0.5 * dPhi +// int( Phi / dPhi - 0.5 ) otherwise + if ( NPhi > 1 ) + { + Phi = ( ( x0 == 0.0 ) && ( y0 == 0.0 ) ) + ? 0.0 + : ( ( y0 >= 0.0 ) ? atan2( y0, x0 ) : atan2( y0, x0 ) + 2.0 * M_PI ); + Idx_Phi = int( Phi / dPhi + 0.5 ) - 1; + +// deal with the case Phi < 0.5 * dPhi + if ( Idx_Phi == -1 ) { Idx_Phi += NPhi; Phi += 2.0 * M_PI; } + + zs[0] = ( (real)0.5 + (real)Idx_Phi ) * dPhi; + zs[1] = ( (real)1.5 + (real)Idx_Phi ) * dPhi; + } + + +// (2-1) prepare the leakage data using linear interpolation +# ifdef __CUDACC__ + real *Ray_tau = SrcTerms->Leakage_tau_DevPtr; + real *Ray_chi = SrcTerms->Leakage_chi_DevPtr; + real *Ray_Flux = SrcTerms->Leakage_HeatFlux_DevPtr; + real *Ray_ERms = SrcTerms->Leakage_HeatERms_DevPtr; + real *Ray_EAve = SrcTerms->Leakage_HeatEAve_DevPtr; +# else + real *Ray_tau = h_SrcLeakage_tau; + real *Ray_chi = h_SrcLeakage_chi; + real *Ray_Flux = h_SrcLeakage_HeatFlux; + real *Ray_ERms = h_SrcLeakage_HeatERms; + real *Ray_EAve = h_SrcLeakage_HeatEAve; +# endif + + real tau[NType_Neutrino], chi[NType_Neutrino], Heat_Flux[NType_Neutrino]; + real Heat_EAve[NType_Neutrino], Heat_ERms[NType_Neutrino]; + + switch ( RayConfig ) + { + case 0: // NTheta != 1, NPhi != 1 + { + real tau_tmp[8], chi_tmp[8], Flux_tmp[8], ERms_tmp[4], EAve_tmp[4]; + + for (int n=0; n= 0 ) + ? ( j == NTheta ? NTheta-1 : j ) + : 0; + const int kk = ( j == -1 || j == NTheta ) + ? k + NPhi_half + : k; + + const int IRay = jj + NTheta * ( kk % NPhi ); + const int Idx_In_2D = n + IRay * NType_Neutrino; + + for (int i=Idx_Rad; i<=Idx_Rad+1; i++, Idx_Out_3D++) + { + const int Idx_In_3D = Idx_In_2D + i * Stride; + + tau_tmp [Idx_Out_3D] = Ray_tau [Idx_In_3D]; + chi_tmp [Idx_Out_3D] = Ray_chi [Idx_In_3D]; + Flux_tmp[Idx_Out_3D] = Ray_Flux[Idx_In_3D]; + } + + ERms_tmp[Idx_Out_2D] = Ray_ERms[Idx_In_2D]; + EAve_tmp[Idx_Out_2D] = Ray_EAve[Idx_In_2D]; + } + + tau [n] = Src_Leakage_TrilinearInterp( tau_tmp , xs, ys, zs, Rad, Tht, Phi ); + chi [n] = Src_Leakage_TrilinearInterp( chi_tmp , xs, ys, zs, Rad, Tht, Phi ); + Heat_Flux[n] = Src_Leakage_TrilinearInterp( Flux_tmp, xs, ys, zs, Rad, Tht, Phi ); + Heat_ERms[n] = Src_Leakage_BilinearInterp ( ERms_tmp, ys, zs, Tht, Phi ); + Heat_EAve[n] = Src_Leakage_BilinearInterp ( EAve_tmp, ys, zs, Tht, Phi ); + } // for (int n=0; nGuessHTilde_FuncPtr, + EoS->HTilde2Temp_FuncPtr, EoS->AuxArrayDevPtr_Flt, + EoS->AuxArrayDevPtr_Int, EoS->Table ); +# else + const real Eint_Code = Hydro_Con2Eint( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], fluid[ENGY], + true, MinEint, PassiveFloor, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, + EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); +# endif + +# ifdef YE + const real Ye = fluid[YE] / fluid[DENS]; +# else + const real Ye = NULL_REAL; +# endif + +# ifdef TEMP_IG + const real Temp_IG_Kelv = fluid[TEMP_IG]; +# else + const real Temp_IG_Kelv = NULL_REAL; +# endif + +# if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) + const int NTarget = 0; +# else + const int NTarget = 1; +# endif + int In_Int[NTarget+1]; + real In_Flt[4], Out[NTarget+1]; + + In_Flt[0] = Dens_Code; + In_Flt[1] = Eint_Code; + In_Flt[2] = Ye; + In_Flt[3] = Temp_IG_Kelv; + + In_Int[0] = NTarget; +# if ( NUC_TABLE_MODE == NUC_TABLE_MODE_ENGY ) + In_Int[1] = NUC_VAR_IDX_EORT; +# endif + +# ifdef __CUDACC__ + EoS->General_FuncPtr( NUC_MODE_ENGY, Out, In_Flt, In_Int, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); +# else + EoS_General_CPUPtr ( NUC_MODE_ENGY, Out, In_Flt, In_Int, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); +# endif + + const real Temp_Kelv = Out[0]; + + +// (3) compute the local heating rates for internal energy and electron fraction + const bool NuHeat = SrcTerms->Leakage_NuHeat; + const real NuHeat_Fac = SrcTerms->Leakage_NuHeat_Fac; + const real Unit_D = SrcTerms->Unit_D; + const real Unit_T = SrcTerms->Unit_T; + real Lum[NType_Neutrino], Heat[NType_Neutrino], NetHeat[NType_Neutrino]; + real dEdt_CGS, dYedt; + + Src_Leakage_ComputeLeak( Dens_Code, Temp_Kelv, Ye, chi, tau, Heat_Flux, Heat_ERms, Heat_EAve, + &dEdt_CGS, &dYedt, Lum, Heat, NetHeat, NuHeat, NuHeat_Fac, Unit_D, EoS ); + +// convert the heating rates to code unit +// --> the returned dEdt_CGS is in erg/g/s and dYedt is in 1/s + const real dEdt_Code = dEdt_CGS * sEint2Code * Unit_T * Dens_Code; + const real dYedt_Code = dYedt * Unit_T; + + +// (4) for recording mode, store the dEdt, luminosity, heating rate, and net heating rate, then return + if ( Mode == LEAK_MODE_RECORD ) + { + fluid[DENS ] = dEdt_CGS * Dens_Code * Unit_D; + fluid[MOMX ] = Lum [0]; + fluid[MOMY ] = Lum [1]; + fluid[MOMZ ] = Lum [2]; + fluid[ENGY ] = Heat[0]; +# ifdef YE + fluid[YE ] = Heat[1]; +# endif +# ifdef DYEDT_NU + fluid[DEDT_NU ] = NetHeat[0]; + fluid[DYEDT_NU] = NetHeat[1]; +# endif + + return; + } + + +// (5-1) store the heating rate +# ifdef DYEDT_NU + fluid[DEDT_NU ] = dEdt_Code; + fluid[DYEDT_NU] = dYedt_Code; +# endif + +// (5-2) check whether the updated Ye is at least 1% away from the table boundary +// --> skip this cell if it is not + const real dYe = dYedt_Code * dt; + + if ( ( Ye + dYe < YeMin ) || ( Ye + dYe > YeMax ) ) + { +# ifdef GAMER_DEBUG + printf( "Invalid heating rate of Ye: Ye=%13.7e, dYedt_Code=%13.7e, dt_Code=%13.7e, YeMin=%13.7e, YeMax=%13.7e\n", + Ye, dYedt_Code, dt, YeMin, YeMax ); +# endif + + return; + } + +// (5-3) verify whether the updated internal energy density and Ye are valid for the nuclear EoS table + const real Eint_Update = Eint_Code + dEdt_Code * dt; + const real Ye_Update = Ye + dYedt_Code * dt; + + In_Flt[1] = Eint_Update; + In_Flt[2] = Ye_Update; + +# ifdef __CUDACC__ + EoS->General_FuncPtr( NUC_MODE_ENGY, Out, In_Flt, In_Int, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); +# else + EoS_General_CPUPtr ( NUC_MODE_ENGY, Out, In_Flt, In_Int, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); +# endif + + const real Temp_Chk = Out[0]; + +// (5-4) update the internal energy density and Ye if the EoS solver succeeds (Temp_Chk != NAN) + if ( Temp_Chk == Temp_Chk ) + { + fluid[ENGY] = Hydro_ConEint2Etot( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], Eint_Update, Emag ); +# ifdef YE + fluid[YE ] = Ye_Update * fluid[DENS]; +# endif + } + + +// (6) final check +# ifdef GAMER_DEBUG + if ( Hydro_IsUnphysical_Single( Eint_Update, "output internal energy density", + (real)0.0, __FLT_MAX__, ERROR_INFO, UNPHY_VERBOSE ) ) + { + printf( " Dens=%13.7e code units, Eint=%13.7e code units, Ye=%13.7e\n", Dens_Code, Eint_Code, Ye ); + printf( " Radius=%13.7e cm, Temp=%13.7e Kelvin, dEdt=%13.7e, dYedt=%13.7e\n", Rad * SrcTerms->Unit_L, Temp_Kelv, dEdt_Code, dYedt ); + + for (int n=0; n The major source-term function will update the system from TimeOld to TimeNew +// dt : Time interval to advance solution +// --> Physical coordinates : TimeNew - TimeOld == dt +// Comoving coordinates : TimeNew - TimeOld == delta(scale factor) != dt +// AuxArray_Flt/Int : Auxiliary arrays +// --> Can be used and/or modified here +// --> Must call Src_SetConstMemory_Leakage() after modification +// +// Return : AuxArray_Flt/Int[] +//------------------------------------------------------------------------------------------------------- +#ifndef __CUDACC__ +void Src_WorkBeforeMajorFunc_Leakage( const int lv, const double TimeNew, const double TimeOld, const double dt, + double AuxArray_Flt[], int AuxArray_Int[] ) +{ + + const int NRad = SrcTerms.Leakage_NRadius; + const int NTheta = SrcTerms.Leakage_NTheta; + const int NPhi = SrcTerms.Leakage_NPhi; + const int NRay = NTheta * NPhi; + + +// (0) do nothing if the leakage profile for current step is already constructed + if ( Step == Leakage_Step ) return; + + +// (1) set the leakage center to the GREP center, if available + double Leakage_Center[3]; + +# ifdef GRAVITY + for (int i=0; i<3; i++) Leakage_Center[i] = GREP_Center[i]; +# else + for (int i=0; i<3; i++) Leakage_Center[i] = amr->BoxCenter[i]; +# endif + + +// (2) set up the coordinate of leakage rays +// --> in linear scale for r <= RadiusMin_Log +// log RadiusMin_Log < r < MaxRadius + const double BinSize_Linear = AuxArray_Flt[SRC_AUX_DRAD ]; + const double MaxRadius = AuxArray_Flt[SRC_AUX_MAXRADIUS ]; + const int NRad_Linear = AuxArray_Int[SRC_AUX_NRAD_LIN ]; + const double RadiusMin_Log = AuxArray_Flt[SRC_AUX_RADMIN_LOG]; + + const double Factor = Src_Leakage_ConstructSeries( NRad - NRad_Linear, RadiusMin_Log, MaxRadius, BinSize_Linear ); + double BinSize_Log = BinSize_Linear; + double Edge[NRad+1]; + + Edge[0] = 0.0; + + for (int i=1; i<=NRad_Linear; i++) + Edge[i] = Edge[i-1] + BinSize_Linear; + + for (int i=NRad_Linear+1; i<=NRad; i++) + { + BinSize_Log *= Factor; + Edge[i] = Edge[i-1] + BinSize_Log; + } + +// compute and store radius in the host array with typecasting + for (int i=0; i prepare data at TimeOld because data on higher level at TimeNew are not available +# ifndef _YE + const long _YE = NULL_INT; +# endif + const int NProf = 3; + const int NData = NRad * NTheta * NPhi; + long TVar[] = { _DENS, _NONE, _YE }; + + Profile_t Ray_Dens_Code, Ray_Temp_Kelv, Ray_Ye; + Profile_t *Leakage_Ray[] = { &Ray_Dens_Code, &Ray_Temp_Kelv, &Ray_Ye }; + + switch ( SrcTerms.Leakage_Opt_Temp ) + { + case LEAK_TEMP_INDCELL: TVar[1] = _TEMP; break; + case LEAK_TEMP_BINDATA: TVar[1] = _EINT; break; + default: + Aux_Error( ERROR_INFO, "unsupported temperature computation scheme %s = %d !!\n", + "SRC_LEAKAGE_OPT_TEMP", SrcTerms.Leakage_Opt_Temp ); + } // switch ( SrcTerms.Leakage_Opt_Temp ) + + +# ifdef LEAKAGE_SEQUENTIAL_MAPPING + for (int i=0; iFreeMemory(); + + +// (8) update the time stamp + Leakage_Step = Step; + +} // FUNCTION : Src_WorkBeforeMajorFunc_Leakage +#endif // #ifndef __CUDACC__ + + + +#ifdef __CUDACC__ +//------------------------------------------------------------------------------------------------------- +// Function : Src_PassData2GPU_Leakage +// Description : Transfer data to GPU +// +// Note : 1. Invoked by Src_WorkBeforeMajorFunc_Leakage() +// 2. Use synchronous transfer +// +// Parameter : None +// +// Return : None +//------------------------------------------------------------------------------------------------------- +void Src_PassData2GPU_Leakage() +{ + + const int NRadius = SrcTerms.Leakage_NRadius; + const int NTheta = SrcTerms.Leakage_NTheta; + const int NPhi = SrcTerms.Leakage_NPhi; + + const long Size_Radius = sizeof(real)*NRadius; + const long Size_tau = sizeof(real)*NRadius*NTheta*NPhi*NType_Neutrino; + const long Size_chi = sizeof(real)*NRadius*NTheta*NPhi*NType_Neutrino; + const long Size_HeatFlux = sizeof(real)*NRadius*NTheta*NPhi*NType_Neutrino; + const long Size_HeatERms = sizeof(real)* NTheta*NPhi*NType_Neutrino; + const long Size_HeatEAve = sizeof(real)* NTheta*NPhi*NType_Neutrino; + +// use synchronous transfer + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcLeakage_Radius, h_SrcLeakage_Radius, Size_Radius, cudaMemcpyHostToDevice ) ); + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcLeakage_tau, h_SrcLeakage_tau, Size_tau, cudaMemcpyHostToDevice ) ); + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcLeakage_chi, h_SrcLeakage_chi, Size_chi, cudaMemcpyHostToDevice ) ); + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcLeakage_HeatFlux, h_SrcLeakage_HeatFlux, Size_HeatFlux, cudaMemcpyHostToDevice ) ); + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcLeakage_HeatERms, h_SrcLeakage_HeatERms, Size_HeatERms, cudaMemcpyHostToDevice ) ); + CUDA_CHECK_ERROR( cudaMemcpy( d_SrcLeakage_HeatEAve, h_SrcLeakage_HeatEAve, Size_HeatEAve, cudaMemcpyHostToDevice ) ); + +} // FUNCTION : Src_PassData2GPU_Leakage +#endif // #ifdef __CUDACC__ + + + +// ================================ +// IV. Set initialization functions +// ================================ + +#ifdef __CUDACC__ +# define FUNC_SPACE __device__ static +#else +# define FUNC_SPACE static +#endif + +FUNC_SPACE SrcFunc_t SrcFunc_Ptr = Src_Leakage; + +//----------------------------------------------------------------------------------------- +// Function : Src_SetCPU/GPUFunc_Leakage +// Description : Return the function pointer of the CPU/GPU source-term function +// +// Note : 1. Invoked by Src_Init_Leakage() +// 2. Call-by-reference +// +// Parameter : SrcFunc_CPU/GPUPtr : CPU/GPU function pointer to be set +// +// Return : SrcFunc_CPU/GPUPtr +//----------------------------------------------------------------------------------------- +#ifdef __CUDACC__ +__host__ +void Src_SetGPUFunc_Leakage( SrcFunc_t &SrcFunc_GPUPtr ) +{ + CUDA_CHECK_ERROR( cudaMemcpyFromSymbol( &SrcFunc_GPUPtr, SrcFunc_Ptr, sizeof(SrcFunc_t) ) ); +} + +#else + +void Src_SetCPUFunc_Leakage( SrcFunc_t &SrcFunc_CPUPtr ) +{ + SrcFunc_CPUPtr = SrcFunc_Ptr; +} + +#endif // #ifdef __CUDACC__ ... else ... + + + +#ifdef __CUDACC__ +//------------------------------------------------------------------------------------------------------- +// Function : Src_SetConstMemory_Leakage +// Description : Set the constant memory variables on GPU +// +// Note : 1. Adopt the suggested approach for CUDA version >= 5.0 +// 2. Invoked by Src_Init_Leakage() and, if necessary, Src_WorkBeforeMajorFunc_Leakage() +// 3. SRC_NAUX_LEAKAGE is defined in Macro.h +// +// Parameter : AuxArray_Flt/Int : Auxiliary arrays to be copied to the constant memory +// DevPtr_Flt/Int : Pointers to store the addresses of constant memory arrays +// +// Return : c_Src_Leakage_AuxArray_Flt[], c_Src_Leakage_AuxArray_Int[], DevPtr_Flt, DevPtr_Int +//--------------------------------------------------------------------------------------------------- +void Src_SetConstMemory_Leakage( const double AuxArray_Flt[], const int AuxArray_Int[], + double *&DevPtr_Flt, int *&DevPtr_Int ) +{ + +// copy data to constant memory + CUDA_CHECK_ERROR( cudaMemcpyToSymbol( c_Src_Leakage_AuxArray_Flt, AuxArray_Flt, SRC_NAUX_LEAKAGE*sizeof(double) ) ); + CUDA_CHECK_ERROR( cudaMemcpyToSymbol( c_Src_Leakage_AuxArray_Int, AuxArray_Int, SRC_NAUX_LEAKAGE*sizeof(int ) ) ); + +// obtain the constant-memory pointers + CUDA_CHECK_ERROR( cudaGetSymbolAddress( (void **)&DevPtr_Flt, c_Src_Leakage_AuxArray_Flt ) ); + CUDA_CHECK_ERROR( cudaGetSymbolAddress( (void **)&DevPtr_Int, c_Src_Leakage_AuxArray_Int ) ); + +} // FUNCTION : Src_SetConstMemory_Leakage +#endif // #ifdef __CUDACC__ + + + +#ifndef __CUDACC__ + +//----------------------------------------------------------------------------------------- +// Function : Src_Init_Leakage +// Description : Initialize the leakage source term +// +// Note : 1. Set auxiliary arrays by invoking Src_SetAuxArray_*() +// --> Copy to the GPU constant memory and store the associated addresses +// 2. Set the source-term function by invoking Src_SetCPU/GPUFunc_*() +// 3. Invoked by Src_Init() +// 4. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : None +// +// Return : None +//----------------------------------------------------------------------------------------- +void Src_Init_Leakage() +{ + +// set the auxiliary arrays + Src_SetAuxArray_Leakage( Src_Leakage_AuxArray_Flt, Src_Leakage_AuxArray_Int ); + +// copy the auxiliary arrays to the GPU constant memory and store the associated addresses +# ifdef GPU + Src_SetConstMemory_Leakage( Src_Leakage_AuxArray_Flt, Src_Leakage_AuxArray_Int, + SrcTerms.Leakage_AuxArrayDevPtr_Flt, SrcTerms.Leakage_AuxArrayDevPtr_Int ); +# else + SrcTerms.Leakage_AuxArrayDevPtr_Flt = Src_Leakage_AuxArray_Flt; + SrcTerms.Leakage_AuxArrayDevPtr_Int = Src_Leakage_AuxArray_Int; +# endif + +// set the major source-term function + Src_SetCPUFunc_Leakage( SrcTerms.Leakage_CPUPtr ); + +# ifdef GPU + Src_SetGPUFunc_Leakage( SrcTerms.Leakage_GPUPtr ); + SrcTerms.Leakage_FuncPtr = SrcTerms.Leakage_GPUPtr; +# else + SrcTerms.Leakage_FuncPtr = SrcTerms.Leakage_CPUPtr; +# endif + +} // FUNCTION : Src_Init_Leakage + + + +//----------------------------------------------------------------------------------------- +// Function : Src_End_Leakage +// Description : Release the resources used by the leakage source term +// +// Note : 1. Invoked by Src_End() +// 2. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : None +// +// Return : None +//----------------------------------------------------------------------------------------- +void Src_End_Leakage() +{ + +// not used by this source term + +} // FUNCTION : Src_End_Leakage + +#endif // #ifndef __CUDACC__ + + + +//------------------------------------------------------------------------------------------------------- +// Function : Src_Leakage_ConstructSeries +// Description : Construct a series with NBin+1 points in which the spacing increases logarithmically +// --> dx * Factor, dx * Factor^2, ..., dx * Factor^NBin +// +// Note : 1. Invoked by Src_WorkBeforeMajorFunc_Leakage() +// 2. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// +// Parameter : NBin : Number of bin +// xmin : The minimum value of constructed series +// xmax : The maximum value of constructed series +// dx : Initial value of spacing +// +// Return : Factor +//------------------------------------------------------------------------------------------------------- +#ifndef __CUDACC__ +double Src_Leakage_ConstructSeries( const int NBin, const double xmin, const double xmax, const double dx ) +{ + + const double Tolerance = 1.0e-10; + const double Target = ( xmax - xmin ) / dx; + const double NBin_double = double(NBin); + int NIter = 100; + double Factor, Correction, Sum, Derivative; + + +// estimate the factor + Factor = exp( log(Target) / NBin_double ); + +// solve ( xmax - xmin ) / dx = Factor + Factor^2 + ... + Factor^NBin + while ( NIter-- ) + { +// analytical formula of Sum = Factor + Factor^2 + ... + Factor^NBin +// Derivative = d(Sum) / d(Factor) + Sum = ( pow( Factor, NBin_double + 1.0 ) - Factor ) / ( Factor - 1.0 ); + Derivative = ( NBin_double + 1.0 ) * Sum / Factor + ( NBin_double - Sum ) / ( Factor - 1.0 ); + Correction = ( Sum - Target ) / Derivative; + + if ( fabs( Correction / Factor ) < Tolerance ) break; + + Factor -= Correction; + } + + + return Factor; + +} // FUNCTION : Src_Leakage_ConstructSeries +#endif + + + +//------------------------------------------------------------------------------------------------------- +// Function : Src_Leakage_BinarySearch +// Description : Use binary search to find the proper array index "Idx" in the input "Array" satisfying +// +// Array[Idx] <= Key < Array[Idx+1] +// +// Note : 1. A variant function of Mis_BinarySearch_Real(), which is shared by both CPU and GPU +// 2. Return "Min" instead if "Key < Array[Min]" +// "Max" instead if "Key > Array[Max]" +// +// Parameter : Array : Sorted look-up array (in ascending numerical order) +// Min : Minimum array index for searching +// Max : Maximum array index for searching +// Key : Target value to search for +// +// Return : Idx : if target is found +// Min : if Key < Array[Min] +// Max : if Key >= Array[Max] +//------------------------------------------------------------------------------------------------------- +GPU_DEVICE static +int Src_Leakage_BinarySearch( const real Array[], int Min, int Max, const real Key ) +{ + +// check whether the input key lies outside the target range + if ( Key < Array[Min] ) return Min; + if ( Key >= Array[Max] ) return Max; + +// binary search + int Idx = -2; + + while ( ( Idx=(Min+Max)/2 ) != Min ) + { + if ( Array[Idx] > Key ) Max = Idx; + else Min = Idx; + } + + + return Idx; + +} // FUNCTION : Src_Leakage_BinarySearch + + + +//------------------------------------------------------------------------------------- +// Function : Src_Leakage_LinearInterp +// Description : Linear interpolation +// +// Note : 1. Return array[0] instead if x_in < xs[0] +// array[1] instead if x_in > xs[1] +//------------------------------------------------------------------------------------- +GPU_DEVICE static +real Src_Leakage_LinearInterp( const real *array, const real *xs, const real x_in ) +{ + + if ( x_in < xs[0] ) return array[0]; + if ( x_in > xs[1] ) return array[1]; + + const real frac_x = ( x_in - xs[0] ) / ( xs[1] - xs[0] ); + + + return array[0] + ( array[1] - array[0] ) * frac_x; + +} // FUNCTION : Src_Leakage_LinearInterp + + + +//------------------------------------------------------------------------------------- +// Function : Src_Leakage_BilinearInterp +// Description : Bilinear interpolation +//------------------------------------------------------------------------------------- +GPU_DEVICE static +real Src_Leakage_BilinearInterp( const real *array, const real *xs, const real *ys, + const real x_in, const real y_in ) +{ + + real v[2] = { Src_Leakage_LinearInterp( array, xs, x_in ), + Src_Leakage_LinearInterp( array+2, xs, x_in ) }; + + + return Src_Leakage_LinearInterp( v, ys, y_in ); + +} // FUNCTION : Src_Leakage_BilinearInterp + + + +//------------------------------------------------------------------------------------- +// Function : Src_Leakage_TrilinearInterp +// Description : Trilinear interpolation +//------------------------------------------------------------------------------------- +GPU_DEVICE static +real Src_Leakage_TrilinearInterp( const real *array, const real *xs, const real *ys, const real *zs, + const real x_in, const real y_in, const real z_in ) +{ + + real v[2] = { Src_Leakage_BilinearInterp( array, xs, ys, x_in, y_in ), + Src_Leakage_BilinearInterp( array+4, xs, ys, x_in, y_in ) }; + + + return Src_Leakage_LinearInterp( v, zs, z_in ); + +} // FUNCTION : Src_Leakage_TrilinearInterp + + + +#endif // #if ( MODEL == HYDRO ) diff --git a/src/SourceTerms/Leakage/CPU_Src_Leakage_ComputeLeak.cpp b/src/SourceTerms/Leakage/CPU_Src_Leakage_ComputeLeak.cpp new file mode 100644 index 000000000..1bffdbf73 --- /dev/null +++ b/src/SourceTerms/Leakage/CPU_Src_Leakage_ComputeLeak.cpp @@ -0,0 +1,1222 @@ +#include "CUFLU.h" +#include "NuclearEoS.h" +#include "Src_Leakage.h" + + + +static const real ONETHIRD = 1.0 / 3.0; +static const real ONESIXTH = 1.0 / 6.0; +#ifdef FLOAT8 +static const real EXP_OVERFLOW_THRESHOLD = 709.782712893384; +#else +static const real EXP_OVERFLOW_THRESHOLD = 88.722839; +#endif + + +#ifndef __CUDACC__ + +// not apply any correction to the leakage profiles for the stability issue +//#define LEAKAGE_PROF_CORR + +extern double CCSN_Leakage_EAve [NType_Neutrino]; +extern double CCSN_Leakage_RadNS[NType_Neutrino]; + +#endif + + +GPU_DEVICE +void Src_Leakage_ComputeLeak( const real Dens_Code, const real Temp_Kelv, const real Ye, const real chi[], const real tau[], + const real Heat_Flux[], const real *Heat_ERms, const real *Heat_EAve, + real *dEdt, real *dYedt, real *Lum, real *Heat, real *NetHeat, + const bool NuHeat, const real NuHeat_Fac, const real UNIT_D, const EoS_t *EoS ); + +GPU_DEVICE static +real Compute_FermiIntegral( const int Order, const real eta ); + + + + +#ifndef __CUDACC__ +//------------------------------------------------------------------------------------------------------- +// Function : Src_Leakage_ComputeTau +// Description : Compute the neutrino optical depth using the Ruffert and Rosswog approaches, +// as well as the neutrino flux and mean/RMS neutrino energy for neutrino heating +// +// Note : 1. Invoked by Src_WorkBeforeMajorFunc_Leakage() +// 2. The stored Heat_Flux is scaled by Const_hc_MeVcm_CUBE to avoid overflow +// 3. Add "#ifndef __CUDACC__" since this routine is only useful on CPU +// 4. Ref: M. Rufferrt, H.-Th. Janka, G. Schaefer, 1996, A&A, 311, 532 (arXiv: 9509006) +// S. Rosswog & M. Liebendoerfer, 2003, MNRAS, 342, 673 (arXiv: 0302301) +// E. O'Connor & C. D. Ott, 2010, Class. Quantum Grav., 27, 114103 (arXiv: 0912.2393) +// +// Parameter : Ray : Profile_t object array stores the profiles of all rays +// --> Density : in code unit +// --> Temperature : in Kelvin +// --> Ye : dimensionless +// Edge : Edge of bins in the radial direction in code unit +// NRadius : Number of bins in the radial direction +// NTheta : Number of bins in the polar direction +// NPhi : Number of bins in the azimuthal direction +// tau_Ruff : Array to store the opacity calculated via the Ruffert scheme +// chi_Ross : Array to store the energy-independent optical depth calculated via the Rosswog scheme +// Heat_Flux : Array to store the neutrino flux +// Heat_ERms : Array to store the rms neutrino energy at neutrino sphere +// Heat_EAve : Array to store the mean neutrino energy at neutrino sphere +// +// Return : tau_Ruff, chi_Ross, Heat_Flux, Heat_ERms, Heat_EAve +//------------------------------------------------------------------------------------------------------- +void Src_Leakage_ComputeTau( Profile_t *Ray[], double *Edge, + const int NRadius, const int NTheta, const int NPhi, + real *tau_Ruff, real *chi_Ross, real *Heat_Flux, + real *Heat_ERms, real *Heat_EAve ) +{ + +# ifdef FLOAT8 + const double Tolerance_Leak = 1.0e-10; +# else + const double Tolerance_Leak = 1.0e-5; +# endif + +# if ( NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) + const real EoS_TempMin_Kelv = POW( (real)10.0, h_EoS_Table[NUC_TAB_TORE ][0] ) / Kelvin2MeV; +# else + const real EoS_TempMin_Kelv = POW( (real)10.0, h_EoS_Table[NUC_TAB_EORT_MODE][0] ) / Kelvin2MeV; +# endif + +# ifdef LEAKAGE_PROF_CORR + const real Dens_CorrThresh_Code = 1.0e3 / UNIT_D; +# endif + + const bool NuHeat = SrcTerms.Leakage_NuHeat; + const real NuHeat_Fac = SrcTerms.Leakage_NuHeat_Fac; + + +// prepare the line element, surface area at bin center, and bin volume + double Edge_CGS[NRadius+1], BinWidth[NRadius], BinSurfArea[NRadius], BinVol[NRadius]; + + for (int i=0; i= IRay_start && j < IRay_end ) continue; + + for (int k=0; kNCell[Idx]; + Dens_Code[TID][i] = Ray[0]->Data [Idx]; + Temp_Kelv[TID][i] = Ray[1]->Data [Idx]; + Ye [TID][i] = Ray[2]->Data [Idx]; + + Dens_CGS [TID][i] = Dens_Code[TID][i] * UNIT_D; + Temp_MeV [TID][i] = Temp_Kelv[TID][i] * Kelvin2MeV; + + if ( NCell[TID][i] == 0L ) + Aux_Error( ERROR_INFO, "empty bin in Ray (%d) in %s !!\n", j, __FUNCTION__ ); + } + +# ifdef LEAKAGE_PROF_CORR +// (1-2) fix potential undershoots near shock, adopted from FLASH + for (int i=1; i 0.53 ) ) + Ye[TID][i] = 0.5 * ( Ye[TID][i-1] + Ye[TID][i+1] ); + } +# endif + +// (1-3) unit conversion + for (int i=0; i not need to include mass difference term (Q/T) +// since we have rest masses in the chemical potentials + eta_nu[TID][i][0] = eta_e[TID][i] - eta_n[TID][i] + eta_p[TID][i]; // (A5) in Ruffert et. al (1996) + eta_nu[TID][i][1] = -eta_nu[TID][i][0]; // (A5) in Ruffert et. al (1996) + eta_nu[TID][i][2] = 0.0; + } // for (int i=0; i assume the optical depth is zero at the right edge of outermost bin + for (int k=0; k=0; i--) + tau[TID][i][k] = tau[TID][i+1][k] + + 0.5 * kappa_tot[TID][i ][k] * BinWidth[i ] + + 0.5 * kappa_tot[TID][i+1][k] * BinWidth[i+1]; + } + +// use new optical depth to update opacity + for (int i=0; i the eta^0_nu is set to 0.0 + eta_nu_loc[TID][i][0] = eta_nu[TID][i][0] * ( 1.0 - exp( -tau[TID][i][0] ) ); // (A3) + eta_nu_loc[TID][i][1] = eta_nu[TID][i][1] * ( 1.0 - exp( -tau[TID][i][1] ) ); // (A4) + eta_nu_loc[TID][i][2] = 0.0; // (A2) + +// number fraction with Pauli blocking effects (Y_NN), assumed completely dissociated + Yn = ( 1.0 - Ye[TID][i] ) / ( 1.0 + TwoThirds * fmax( 0.0, eta_n[TID][i] ) ); // (A8) + Yp = Ye[TID][i] / ( 1.0 + TwoThirds * fmax( 0.0, eta_p[TID][i] ) ); // (A8) + +// number fraction with Fermion blocking effects + if ( x_h[TID][i] < 0.5 ) + { + fac1 = exp( -eta_hat[TID][i] ); + + Ynp = ( 2.0 * Ye[TID][i] - 1.0 ) / ( fac1 - 1.0 ); // (A13) + Ypn = fac1 * Ynp; // (A14) + } + + else + { + Ynp = x_n[TID][i]; + Ypn = x_p[TID][i]; + } + + Ynp = fmax( 0.0, Ynp ); + Ypn = fmax( 0.0, Ypn ); + +// update opacity +// --> electron neutrino (nu_e) + fac1 = Compute_FermiIntegral( 53, eta_nu_loc[TID][i][0] ); + fac2 = 1.0 + exp( eta_e[TID][i] - Compute_FermiIntegral( 54, eta_nu_loc[TID][i][0] ) ); // (A15) + + kappa_scat_n[0] = kappa_scat_n_fac[TID][i] * Yn * fac1; // (A6) + kappa_scat_p[0] = kappa_scat_p_fac[TID][i] * Yp * fac1; // (A6) + kappa_abs_n = kappa_abs_fac [TID][i] * Ynp * fac1 / fac2; // (A11) + +// --> electron anti-neutrino (nu_a) + fac1 = Compute_FermiIntegral( 53, eta_nu_loc[TID][i][1] ); + fac2 = 1.0 + exp( -eta_e[TID][i] - Compute_FermiIntegral( 54, eta_nu_loc[TID][i][1] ) ); // (A16) + + kappa_scat_n[1] = kappa_scat_n_fac[TID][i] * Yn * fac1; // (A6) + kappa_scat_p[1] = kappa_scat_p_fac[TID][i] * Yp * fac1; // (A6) + kappa_abs_p = kappa_abs_fac [TID][i] * Ypn * fac1 / fac2; // (A11) + +// --> heavy-lepton neutrino (nu_x) + fac1 = Compute_FermiIntegral( 53, eta_nu_loc[TID][i][2] ); + + kappa_scat_n[2] = kappa_scat_n_fac[TID][i] * Yn * fac1; // (A6) + kappa_scat_p[2] = kappa_scat_p_fac[TID][i] * Yp * fac1; // (A6) + +// sum each contribution + kappa_tot[TID][i][0] = kappa_scat_p[0] + kappa_scat_n[0] + kappa_abs_n; // nu_e; (A17) + kappa_tot[TID][i][1] = kappa_scat_p[1] + kappa_scat_n[1] + kappa_abs_p; // nu_a; (A18) + kappa_tot[TID][i][2] = kappa_scat_p[2] + kappa_scat_n[2]; // nu_x; (A19) + } // for (int i=0; i Tolerance_Leak ) IsConverged = false; + }} + + NIter++; + } // while ( ! IsConverged && NIter < NIter_Max ) + + if ( ! IsConverged ) + { + Aux_Message( stderr, "\n#%6s %8s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n", + "Index", "NCell", + "kappa_old[0]", "kappa_old[1]", "kappa_old[2]", + "kappa_new[0]", "kappa_new[1]", "kappa_new[2]", + "tau[0]", "tau[1]", "tau[2]" ); + + for (int i=0; i use the neutrino degeneracy obtained from the Ruffert scheme + double blocking_factor_nue, blocking_factor_nua, eta_pn, eta_np; + double kappa_scat_fac_Rosswog, kappa_abs_fac_Rosswog; + double kappa_tilde_scat_n[NType_Neutrino], kappa_tilde_scat_p[NType_Neutrino], kappa_tilde_scat_x[NType_Neutrino]; + double kappa_tilde_abs_n [NType_Neutrino], kappa_tilde_abs_p [NType_Neutrino], kappa_tilde_abs_x [NType_Neutrino]; + +// (4-1) construct the energy-independent opacity (zeta) + for (int i=0; i the power of Abar is 1 because the opacity multiples the number fraction, not mass fraction + kappa_scat_fac_Rosswog *= 0.25 * abar[TID][i] * SQR( 1.0 - zbar[TID][i] / abar[TID][i] ); + + for (int k=0; k the factor Const_NA is moved from eta_pn/eta_np to kappa_abs_fac_Rosswog + kappa_abs_fac_Rosswog = Const_Rosswog_kappa_a; + blocking_factor_nue = 1.0 + exp( eta_e[TID][i] - Compute_FermiIntegral( 54, eta_nu[TID][i][0] ) ); + blocking_factor_nua = 1.0 + exp( -eta_e[TID][i] - Compute_FermiIntegral( 54, eta_nu[TID][i][1] ) ); + + if ( Dens_CGS[TID][i] < 1.0e11 ) + { +// non-degenerate state, use mass fraction as chemical potential + eta_pn = Dens_CGS[TID][i] * x_p[TID][i]; + eta_np = Dens_CGS[TID][i] * x_n[TID][i]; + } + + else + { + fac1 = Dens_CGS[TID][i] * ( x_n[TID][i] - x_p[TID][i] ); + + eta_pn = fac1 / ( exp( eta_hat[TID][i] ) - 1.0 ); // (A9) + eta_np = -fac1 / ( exp( -eta_hat[TID][i] ) - 1.0 ); // (A9) + } + + eta_pn = fmax( 0.0, eta_pn ); + eta_np = fmax( 0.0, eta_np ); + + kappa_tilde_abs_n[0] = eta_np * kappa_abs_fac_Rosswog / blocking_factor_nue; + kappa_tilde_abs_n[1] = 0.0; // no absorption of a-type neutrinos on neutrons + kappa_tilde_abs_n[2] = 0.0; // no absorption of x-type neutrinos on neutrons + kappa_tilde_abs_p[0] = 0.0; // no absorption of e-type neutrinos on protons + kappa_tilde_abs_p[1] = eta_pn * kappa_abs_fac_Rosswog / blocking_factor_nua; + kappa_tilde_abs_p[2] = 0.0; // no absorption of x-type neutrinos on protons + kappa_tilde_abs_x[0] = 0.0; // no absorption on nuclei + kappa_tilde_abs_x[1] = 0.0; // no absorption on nuclei + kappa_tilde_abs_x[2] = 0.0; // no absorption on nuclei + +// sum each contribution to get zeta; (A21) + for (int k=0; k assume chi is zero at the right edge of outermost bin + for (int k=0; k=0; i--) + chi[TID][i][k] = chi[TID][i+1][k] + + 0.5 * zeta[TID][i ][k] * BinWidth[i ] + + 0.5 * zeta[TID][i+1][k] * BinWidth[i+1]; + } + +// store chi + for (int i=0; i the returned Lum represents the luminosity per unit volume + Src_Leakage_ComputeLeak( Dens_Code[TID][i], Temp_Kelv[TID][i], Ye[TID][i], chi_Ross_3D[i][j], tau_Ruff_3D[i][j], + Heat_Flux_3D[i][j], Heat_ERms_2D[j], Heat_EAve_2D[j], + dEdt[TID]+i, dYedt[TID]+i, Lum, Heat, NetHeat, + NuHeat, NuHeat_Fac, UNIT_D, &EoS ); + + for (int k=0; kGeneral_FuncPtr( NUC_MODE_TEMP, Out, In_Flt, In_Int, EoS->AuxArrayDevPtr_Flt, EoS->AuxArrayDevPtr_Int, EoS->Table ); +# else + EoS_General_CPUPtr ( NUC_MODE_TEMP, Out, In_Flt, In_Int, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); +# endif + + const real eta_e = Out[0] / Temp_MeV; + const real eta_p = Out[1] / Temp_MeV; + const real eta_n = Out[2] / Temp_MeV; + const real x_h = Out[3]; + const real x_n = Out[4]; + const real x_p = Out[5]; + const real abar = Out[6]; + const real zbar = Out[7]; + const real eta_hat = eta_n - eta_p - Const_Qnp / Temp_MeV; + real eta_nu[NType_Neutrino]; + +// not need to include mass difference term (Q/T) since we have rest masses in the chemical potentials + eta_nu[0] = eta_e - eta_n + eta_p; + eta_nu[1] = -eta_nu[0]; + eta_nu[2] = (real)0.0; + +// (1-2) do nothing if outside the shock + if ( ( x_h > (real)0.5 ) && ( Dens_CGS < (real)1.0e13 ) ) + { + *dEdt = (real)0.0; + *dYedt = (real)0.0; + + for (int k=0; k the power of Abar is 1 because the opacity multiples the number fraction, not mass fraction + kappa_scat_fac_Rosswog *= (real)0.25 * abar * SQR( (real)1.0 - zbar / abar ); + + for (int k=0; k the factor Const_NA is moved from eta_pn/eta_np to kappa_abs_fac_Rosswog + kappa_abs_fac_Rosswog = (real)Const_Rosswog_kappa_a; + blocking_factor[0] = (real)1.0 + EXP( eta_e - Compute_FermiIntegral( 54, eta_nu[0] ) ); + blocking_factor[1] = (real)1.0 + EXP( -eta_e - Compute_FermiIntegral( 54, eta_nu[1] ) ); + + if ( Dens_CGS < (real)1.0e11 ) + { +// non-degenerate state, use mass fraction as chemical potential + eta_pn = Dens_CGS * x_p; + eta_np = Dens_CGS * x_n; + } + + else + { + const real factor = Dens_CGS * ( x_n - x_p ); + + eta_pn = factor / ( EXP( eta_hat ) - (real)1.0 ); // (A9) + eta_np = -factor / ( EXP( -eta_hat ) - (real)1.0 ); // (A9) + } + + eta_pn = FMAX( (real)0.0, eta_pn ); + eta_np = FMAX( (real)0.0, eta_np ); + + kappa_tilde_abs_n[0] = eta_np * kappa_abs_fac_Rosswog / blocking_factor[0]; + kappa_tilde_abs_n[1] = (real)0.0; // no absorption of a-type neutrinos on neutrons + kappa_tilde_abs_n[2] = (real)0.0; // no absorption of x-type neutrinos on neutrons + kappa_tilde_abs_p[0] = (real)0.0; // no absorption of e-type neutrinos on protons + kappa_tilde_abs_p[1] = eta_pn * kappa_abs_fac_Rosswog / blocking_factor[1]; + kappa_tilde_abs_p[2] = (real)0.0; // no absorption of x-type neutrinos on protons + kappa_tilde_abs_x[0] = (real)0.0; // no absorption on nuclei + kappa_tilde_abs_x[1] = (real)0.0; // no absorption on nuclei + kappa_tilde_abs_x[2] = (real)0.0; // no absorption on nuclei + +// sum each contribution to get zeta; (A21) + for (int k=0; k the diffusion time is increased by a factor of 2, as suggested by O'Connor & Ott (2010) + real R_diff [NType_Neutrino]; // number diffusion rate per volume + real Q_diff [NType_Neutrino]; // energy diffusion rate per volume + real weight [NType_Neutrino] = { (real)1.0, (real)1.0, (real)4.0 }; // numerical multiplicity factor, g_nu + + for (int k=0; k the factor Const_NA is moved from eta_pn/eta_np to beta + const real beta = (real)Const_leakage_beta; // (A7) + + R_loc[0] = beta * eta_pn * Temp_MeV_Quintic * FermiInte_p[1]; // (A6) + Q_loc[0] = beta * eta_pn * Temp_MeV_Sextic * FermiInte_p[2]; // (A10) + R_loc[1] = beta * eta_np * Temp_MeV_Quintic * FermiInte_m[1]; // (A12) + Q_loc[1] = beta * eta_np * Temp_MeV_Sextic * FermiInte_m[2]; // (A13) + R_loc[2] = (real)0.0; + Q_loc[2] = (real)0.0; + +// (4-2) electron-positron pair annihilation in the Ruffert scheme + const real epsilon_m = Temp_MeV_Quartic * FermiInte_p[0]; // (B5) + const real epsilon_p = Temp_MeV_Quartic * FermiInte_m[0]; // (B5) + const real pair_R_factor = epsilon_m * epsilon_p; // for (B8) and (B10), factored out the constants + const real factor_pair = (real)0.5 * ( FermiInte_p[3] + FermiInte_m[3] ); // (B16) + real R_pair, Q_pair; + + for (int k=0; k * exp(-2 * tau_nu) + + for (int k=0; k (real)1.0e-3 ) + { + const real eta_sqr = SQR(eta); + + switch ( Order ) + { + case 0 : + integral = ( eta > EXP_OVERFLOW_THRESHOLD ) + ? eta + : LOG( (real)1.0 + EXP(eta) ); + break; + + case 1 : + integral = (real)0.5 * eta_sqr + (real)1.6449; + integral /= (real)1.0 + EXP( (real)-1.6855 * eta ); + break; + + case 2 : + integral = ( ONETHIRD * eta_sqr + (real)3.2899 ) * eta; + integral /= (real)1.0 - EXP( (real)-1.8246 * eta ); + break; + + case 3 : + integral = ( (real)0.25 * eta_sqr + (real)4.9348 ) * eta_sqr + (real)11.3644; + integral /= (real)1.0 + EXP( (real)-1.9039 * eta ); + break; + + case 4 : + integral = ( ( (real)0.2 * eta_sqr + (real)6.5797 ) * eta_sqr + (real)45.4576 ) * eta; + integral /= (real)1.0 - EXP( (real)-1.9484 * eta ); + break; + + case 5 : + integral = ( ( ONESIXTH * eta_sqr + (real)8.2247 ) * eta_sqr + (real)113.6439 ) * eta_sqr + (real)236.5323; + integral /= (real)1.0 + EXP( (real)-1.9727 * eta ); + break; + + case 43: // case 4 / case 3 + integral = ( ( (real)0.2 * eta_sqr + (real)6.5797 ) * eta_sqr + (real) 45.4576 ) * eta; + integral /= ( (real)0.25 * eta_sqr + (real)4.9348 ) * eta_sqr + (real) 11.3644; + integral *= ( (real)1.0 + EXP( (real)-1.9039 * eta ) ) + / ( (real)1.0 - EXP( (real)-1.9484 * eta ) ); + break; + + case 53: // case 5 / case 3 + integral = ( ( ONESIXTH * eta_sqr + (real)8.2247 ) * eta_sqr + (real)113.6439 ) * eta_sqr + (real)236.5323; + integral /= ( (real)0.25 * eta_sqr + (real)4.9348 ) * eta_sqr + (real) 11.3644; + integral *= ( (real)1.0 + EXP( (real)-1.9039 * eta ) ) + / ( (real)1.0 + EXP( (real)-1.9727 * eta ) ); + break; + + case 54: // case 5 / case 4 + integral = ( ( ONESIXTH * eta_sqr + (real)8.2247 ) * eta_sqr + (real)113.6439 ) * eta_sqr + (real)236.5323; + integral /= ( ( (real)0.2 * eta_sqr + (real)6.5797 ) * eta_sqr + (real) 45.4576 ) * eta; + integral *= ( (real)1.0 - EXP( (real)-1.9484 * eta ) ) + / ( (real)1.0 + EXP( (real)-1.9727 * eta ) ); + break; + } // switch ( Order ) + } + + else + { + switch ( Order ) + { + case 0 : + integral = LOG( (real)1.0 + EXP(eta) ); + break; + + case 1 : + integral = EXP(eta); + integral /= (real)1.0 + (real)0.2159 * EXP( (real)0.8857 * eta ); + break; + + case 2 : + integral = (real)2.0 * EXP(eta); + integral /= (real)1.0 + (real)0.1092 * EXP( (real)0.8908 * eta ); + break; + + case 3 : + integral = (real)6.0 * EXP(eta); + integral /= (real)1.0 + (real)0.0559 * EXP( (real)0.9069 * eta ); + break; + + case 4 : + integral = (real)24.0 * EXP(eta); + integral /= (real)1.0 + (real)0.0287 * EXP( (real)0.9257 * eta ); + break; + + case 5 : + integral = (real)120.0 * EXP(eta); + integral /= (real)1.0 + (real)0.0147 * EXP( (real)0.9431 * eta ); + break; + + case 43: // case 4 / case 3 + integral = (real)4.0 * ( (real)1.0 + (real)0.0559 * EXP( (real)0.9069 * eta ) ) + / ( (real)1.0 + (real)0.0287 * EXP( (real)0.9257 * eta ) ); + break; + + case 53: // case 5 / case 3 + integral = (real)20.0 * ( (real)1.0 + (real)0.0559 * EXP( (real)0.9069 * eta ) ) + / ( (real)1.0 + (real)0.0147 * EXP( (real)0.9431 * eta ) ); + break; + + case 54: // case 5 / case 4 + integral = (real)5.0 * ( (real)1.0 + (real)0.0287 * EXP( (real)0.9257 * eta ) ) + / ( (real)1.0 + (real)0.0147 * EXP( (real)0.9431 * eta ) ); + break; + } // switch ( Order ) + } // if ( eta > (real)1.0e-3 ) ... else ... + + + return integral; + +} // FUNCTION : Compute_FermiIntegral diff --git a/src/SourceTerms/Leakage/CUSRC_Src_Leakage.cu b/src/SourceTerms/Leakage/CUSRC_Src_Leakage.cu new file mode 120000 index 000000000..308a86af5 --- /dev/null +++ b/src/SourceTerms/Leakage/CUSRC_Src_Leakage.cu @@ -0,0 +1 @@ +CPU_Src_Leakage.cpp \ No newline at end of file diff --git a/src/SourceTerms/Leakage/CUSRC_Src_Leakage_ComputeLeak.cu b/src/SourceTerms/Leakage/CUSRC_Src_Leakage_ComputeLeak.cu new file mode 120000 index 000000000..76cb1bfd2 --- /dev/null +++ b/src/SourceTerms/Leakage/CUSRC_Src_Leakage_ComputeLeak.cu @@ -0,0 +1 @@ +CPU_Src_Leakage_ComputeLeak.cpp \ No newline at end of file diff --git a/src/SourceTerms/Leakage/Src_Leakage.h b/src/SourceTerms/Leakage/Src_Leakage.h new file mode 100644 index 000000000..5cfed5f1d --- /dev/null +++ b/src/SourceTerms/Leakage/Src_Leakage.h @@ -0,0 +1,87 @@ +#ifndef __SRC_LEAKAGE_H__ +#define __SRC_LEAKAGE_H__ + + + +#include "Macro.h" +#include "PhysicalConstant.h" + + +// auxiliary array indices +#define SRC_AUX_PNS_X 0 // AuxArray_Flt: x-coordinate of proto-neutron star +#define SRC_AUX_PNS_Y 1 // AuxArray_Flt: y-coordinate of proto-neutron star +#define SRC_AUX_PNS_Z 2 // AuxArray_Flt: z-coordinate of proto-neutron star +#define SRC_AUX_MAXRADIUS 3 // AuxArray_Flt: maximum radius +#define SRC_AUX_RADMIN_LOG 4 // AuxArray_Flt: minimum radius of logarithmic bins +#define SRC_AUX_DRAD 5 // AuxArray_Flt: spacing linear bins in the radial direction +#define SRC_AUX_DTHETA 6 // AuxArray_Flt: spacing in the polar direction +#define SRC_AUX_DPHI 7 // AuxArray_Flt: spacing in the azimuthal direction +#define SRC_AUX_YEMIN 8 // AuxArray_Flt: minimum allowed Ye after update in the leakage scheme +#define SRC_AUX_YEMAX 9 // AuxArray_Flt: maximum allowed Ye after update in the leakage scheme +#define SRC_AUX_VSQR2CODE 10 // AuxArray_Flt: convert velocity^2 to code unit + +#define SRC_AUX_RAYCONF 0 // AuxArray_Int: configuration of the leakage rays +#define SRC_AUX_NRAD_LIN 1 // AuxArray_Int: number of linear bin in the radial direction +#define SRC_AUX_STRIDE 2 // AuxArray_Int: stride for indexing +#define SRC_AUX_MODE 3 // AuxArray_Int: control the behavior of the leakage scheme + + +// additional physical constants for the leakage scheme +const double Const_Qnp = 1.293333; // m_n - m_p, MeV +const double Const_Cv = 0.5 + 2.0 * 0.23; // vector coupling +const double Const_Ca = 0.5; // axial coupling +const double Const_alpha = 1.23; // g_A +const double Const_me_MeV = 5.10998950e-1; // electron mass in MeV +const double Const_sigma0 = 1.76e-44; // reference weak-interaction cross section in cm^2 +const double Const_mn = 1.67492749804e-24; // neutron mass in gram +const double Const_fsc = 7.2973525693e-3; // fine-structure constant + + +// derived constant +const double Const_gamma_0 = 5.56515284698977e-2; // 2.0 * sqrt( Const_fsc / (3.0*M_PI) ) +const double Const_hc_MeVcm = 1.0e-6 * 2.0 * M_PI * Const_Planck_eV * Const_c; // h*c in MeV*cm +const double Const_alpha_SQR = SQR ( Const_alpha ); +const double Const_me_MeV_SQR = SQR ( Const_me_MeV ); +const double Const_hc_MeVcm_CUBE = CUBE( Const_hc_MeVcm ); + + +// constant factor in the leakage scheme +// --> Ruffert scheme +const double Const_sn_0 = ( 1.0 + 5.0 * Const_alpha_SQR ) / 24.0; // cross section coefficients +const double Const_sp_0 = ( 4.0 * SQR(Const_Cv - 1.0) + 5.0 * Const_alpha_SQR ) / 24.0; // cross section coefficients +const double Const_Ruffert_kappa_sn = Const_sn_0 * Const_NA * Const_sigma0 / Const_me_MeV_SQR; // constant coefficient of scattering process +const double Const_Ruffert_kappa_sp = Const_sp_0 * Const_NA * Const_sigma0 / Const_me_MeV_SQR; // constant coefficient of scattering process +const double Const_Ruffert_kappa_a = 0.25 * ( 1.0 + 3.0 * Const_alpha_SQR ) + * Const_NA * Const_sigma0 / Const_me_MeV_SQR; // constant coefficient of absorption process + +// --> Rosswog scheme +const double Const_Rosswog_kappa_s = 0.25 * Const_NA * Const_sigma0 / Const_me_MeV_SQR; // constant coefficient of scattering process +const double Const_Rosswog_kappa_a = Const_Ruffert_kappa_a; // constant coefficient of absorption process + +// --> leakage scheme (factored out 1.0 / Const_hc_MeVcm_CUBE to avoid overflow) +const double Const_leakage_diff = 4.0 * M_PI * Const_c / 6.0; +const double Const_leakage_beta = M_PI * Const_c * ( 1.0 + 3.0 * Const_alpha_SQR ) + * Const_NA * Const_sigma0 / Const_me_MeV_SQR; +const double Const_leakage_pair = 64.0 * SQR( M_PI ) * Const_sigma0 * Const_c + / ( Const_hc_MeVcm_CUBE * Const_me_MeV_SQR ); +const double Const_leakage_pair_ea = ( SQR( Const_Cv - Const_Ca ) + SQR( Const_Cv + Const_Ca ) ) + * Const_leakage_pair / 36.0; +const double Const_leakage_pair_x = ( SQR( Const_Cv - Const_Ca ) + SQR( Const_Cv + Const_Ca - 2.0 ) ) + * Const_leakage_pair / 9.0; +const double Const_leakage_gamma = CUBE( M_PI ) * Const_sigma0 * Const_c + / ( 3.0 * Const_fsc * Const_me_MeV_SQR * Const_hc_MeVcm_CUBE ); +const double Const_leakage_gamma_ea = Const_leakage_gamma * SQR( Const_Cv ); +const double Const_leakage_gamma_x = Const_leakage_gamma * SQR( Const_Cv - 1.0 ) * 4.0; +const double Const_leakage_brem = 0.5 * 0.231 * ( 2.0778e2 / Const_MeV ) * Const_hc_MeVcm_CUBE; +const double Const_leakage_heat = 0.25 * ( 1.0 + 3.0 * Const_alpha_SQR ) * Const_sigma0 + / ( Const_me_MeV_SQR * Const_mn ); + + +// miscellaneous +const int NType_Neutrino = 3; // electron neutrino (e), electron anti-neutrino (a), heavy-lepton neutrino (x) +const double TwoThirds = 2.0 / 3.0; +const double Kelvin2MeV = 1.0e-6 * Const_kB_eV; // Kelvin to MeV +const double Erg2MeV = 1.0 / Const_MeV; // erg to MeV + + +#endif // #ifndef __SRC_LEAKAGE_H__ diff --git a/src/SourceTerms/LightBulb/CPU_Src_Lightbulb.cpp b/src/SourceTerms/LightBulb/CPU_Src_Lightbulb.cpp index 7d7db5cbf..9a1994ecd 100644 --- a/src/SourceTerms/LightBulb/CPU_Src_Lightbulb.cpp +++ b/src/SourceTerms/LightBulb/CPU_Src_Lightbulb.cpp @@ -213,7 +213,7 @@ static void Src_Lightbulb( real fluid[], const real B[], fluid[ENGY] = Hydro_ConEint2Etot( fluid[DENS], fluid[MOMX], fluid[MOMY], fluid[MOMZ], Eint_Update, Emag ); # ifdef DEDT_NU - fluid[DEDT_NU] = FABS( rate_Code * Dens_Code ); + fluid[DEDT_NU] = rate_Code * Dens_Code; # endif diff --git a/src/SourceTerms/Src_AdvanceDt.cpp b/src/SourceTerms/Src_AdvanceDt.cpp index 0f65e5a4a..272497950 100644 --- a/src/SourceTerms/Src_AdvanceDt.cpp +++ b/src/SourceTerms/Src_AdvanceDt.cpp @@ -2,11 +2,6 @@ -// flag for checking whether the dEdt_Nu field is initialized -bool IsInit_dEdt_Nu = false; - - - //------------------------------------------------------------------------------------------------------- // Function : Src_AdvanceDt @@ -45,7 +40,4 @@ void Src_AdvanceDt( const int lv, const double TimeNew, const double TimeOld, co InvokeSolver( SRC_SOLVER, lv, TimeNew, TimeOld, dt, NULL_REAL, SaveSg_Flu, SaveSg_Mag, NULL_INT, OverlapMPI, Overlap_Sync ); - - if ( SrcTerms.Lightbulb ) IsInit_dEdt_Nu = true; - } // FUNCTION : Src_AdvanceDt diff --git a/src/SourceTerms/Src_End.cpp b/src/SourceTerms/Src_End.cpp index 25ed4ed54..fcd73ba24 100644 --- a/src/SourceTerms/Src_End.cpp +++ b/src/SourceTerms/Src_End.cpp @@ -6,6 +6,7 @@ #if ( MODEL == HYDRO ) void Src_End_Deleptonization(); void Src_End_Lightbulb(); +void Src_End_Leakage(); #endif // this function pointer can be set by a test problem initializer for a non-built-in source term @@ -38,7 +39,10 @@ void Src_End() if ( SrcTerms.Lightbulb ) Src_End_Lightbulb(); -# endif + + if ( SrcTerms.Leakage ) + Src_End_Leakage(); +# endif // if ( MODEL == HYDRO ) // users may not define Src_End_User_Ptr if ( SrcTerms.User && Src_End_User_Ptr ) diff --git a/src/SourceTerms/Src_Init.cpp b/src/SourceTerms/Src_Init.cpp index fb67b0c7a..2bea72810 100644 --- a/src/SourceTerms/Src_Init.cpp +++ b/src/SourceTerms/Src_Init.cpp @@ -6,6 +6,7 @@ #if ( MODEL == HYDRO ) void Src_Init_Deleptonization(); void Src_Init_Lightbulb(); +void Src_Init_Leakage(); #endif // this function pointer can be set by a test problem initializer for a user-specified source term @@ -37,6 +38,7 @@ void Src_Init() # if ( MODEL == HYDRO ) SrcTerms.Deleptonization || SrcTerms.Lightbulb || + SrcTerms.Leakage || # endif SrcTerms.User ) @@ -77,7 +79,16 @@ void Src_Init() # endif SrcTerms.Lightbulb_AuxArrayDevPtr_Flt = NULL; SrcTerms.Lightbulb_AuxArrayDevPtr_Int = NULL; + + SrcTerms.Leakage_FuncPtr = NULL; + SrcTerms.Leakage_CPUPtr = NULL; +# ifdef GPU + SrcTerms.Leakage_GPUPtr = NULL; # endif + SrcTerms.Leakage_AuxArrayDevPtr_Flt = NULL; + SrcTerms.Leakage_AuxArrayDevPtr_Int = NULL; +# endif // #if ( MODEL == HYDRO ) + SrcTerms.User_FuncPtr = NULL; SrcTerms.User_CPUPtr = NULL; # ifdef GPU @@ -114,9 +125,22 @@ void Src_Init() if ( SrcTerms.Lightbulb_GPUPtr == NULL ) Aux_Error( ERROR_INFO, "SrcTerms.Lightbulb_GPUPtr == NULL !!\n" ); # endif } -# endif -// (3) user-specified source term +// (3) leakage + if ( SrcTerms.Leakage ) + { + Src_Init_Leakage(); + +// check if the source-term function is set properly + if ( SrcTerms.Leakage_FuncPtr == NULL ) Aux_Error( ERROR_INFO, "SrcTerms.Leakage_FuncPtr == NULL !!\n" ); + if ( SrcTerms.Leakage_CPUPtr == NULL ) Aux_Error( ERROR_INFO, "SrcTerms.Leakage_CPUPtr == NULL !!\n" ); +# ifdef GPU + if ( SrcTerms.Leakage_GPUPtr == NULL ) Aux_Error( ERROR_INFO, "SrcTerms.Leakage_GPUPtr == NULL !!\n" ); +# endif + } +# endif // if ( MODEL == HYDRO ) + +// (4) user-specified source term if ( SrcTerms.User ) { if ( Src_Init_User_Ptr == NULL ) Aux_Error( ERROR_INFO, "Src_Init_User_Ptr == NULL !!\n" ); diff --git a/src/SourceTerms/Src_WorkBeforeMajorFunc.cpp b/src/SourceTerms/Src_WorkBeforeMajorFunc.cpp index 13bb5bc1a..384c95911 100644 --- a/src/SourceTerms/Src_WorkBeforeMajorFunc.cpp +++ b/src/SourceTerms/Src_WorkBeforeMajorFunc.cpp @@ -9,6 +9,9 @@ void Src_WorkBeforeMajorFunc_Deleptonization( const int lv, const double TimeNew void Src_WorkBeforeMajorFunc_Lightbulb( const int lv, const double TimeNew, const double TimeOld, const double dt, double AuxArray_Flt[], int AuxArray_Int[] ); + +void Src_WorkBeforeMajorFunc_Leakage( const int lv, const double TimeNew, const double TimeOld, const double dt, + double AuxArray_Flt[], int AuxArray_Int[] ); #endif // this function pointer can be set by a test problem initializer for a user-specified source term @@ -46,13 +49,18 @@ void Src_WorkBeforeMajorFunc( const int lv, const double TimeNew, const double T Src_WorkBeforeMajorFunc_Deleptonization( lv, TimeNew, TimeOld, dt, Src_Dlep_AuxArray_Flt, Src_Dlep_AuxArray_Int ); -// (2) Lightbulb +// (2) lightbulb if ( SrcTerms.Lightbulb ) - Src_WorkBeforeMajorFunc_Lightbulb( lv, TimeNew, TimeOld, dt, - Src_Lightbulb_AuxArray_Flt, Src_Lightbulb_AuxArray_Int ); -# endif + Src_WorkBeforeMajorFunc_Lightbulb ( lv, TimeNew, TimeOld, dt, + Src_Lightbulb_AuxArray_Flt, Src_Lightbulb_AuxArray_Int ); + +// (3) leakage + if ( SrcTerms.Leakage && lv == 0 ) + Src_WorkBeforeMajorFunc_Leakage ( lv, TimeNew, TimeOld, dt, + Src_Leakage_AuxArray_Flt, Src_Leakage_AuxArray_Int ); +# endif // #if ( MODEL == HYDRO ) -// (3) user-specified source term +// (4) user-specified source term // --> users may not define Src_WorkBeforeMajorFunc_User_Ptr if ( SrcTerms.User && Src_WorkBeforeMajorFunc_User_Ptr != NULL ) Src_WorkBeforeMajorFunc_User_Ptr ( lv, TimeNew, TimeOld, dt, diff --git a/src/TestProblem/Hydro/CCSN/Flag_CCSN.cpp b/src/TestProblem/Hydro/CCSN/Flag_CCSN.cpp index 7a9c8a88d..00654f104 100644 --- a/src/TestProblem/Hydro/CCSN/Flag_CCSN.cpp +++ b/src/TestProblem/Hydro/CCSN/Flag_CCSN.cpp @@ -9,6 +9,11 @@ extern double CCSN_CC_MaxRefine_Dens1; extern double CCSN_CC_MaxRefine_Dens2; extern double CCSN_CentralDens; +extern bool CCSN_Is_PostBounce; +extern double CCSN_REF_RBase; +extern double CCSN_Rsh_Max; +extern double CCSN_Rsh_Ave; + extern double CCSN_MaxRefine_Rad; extern double CCSN_AngRes_Min; extern double CCSN_AngRes_Max; @@ -85,7 +90,7 @@ bool Flag_CoreCollapse( const int i, const int j, const int k, const int lv, con //------------------------------------------------------------------------------------------------------- -// Function : Flag_Lightbulb +// Function : Flag_PostBounce // Description : Check if the element (i,j,k) of the input data satisfies the user-defined flag criteria // // Note : 1. Invoked by "Flag_Check" using the function pointer "Flag_User_Ptr" @@ -104,7 +109,7 @@ bool Flag_CoreCollapse( const int i, const int j, const int k, const int lv, con // Return : "true" if the flag criteria are satisfied // "false" if the flag criteria are not satisfied //------------------------------------------------------------------------------------------------------- -bool Flag_Lightbulb( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold ) +bool Flag_PostBounce( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold ) { bool Flag = false; @@ -137,7 +142,7 @@ bool Flag_Lightbulb( const int i, const int j, const int k, const int lv, const return Flag; -} // FUNCTION : Flag_Lightbulb +} // FUNCTION : Flag_PostBounce @@ -172,13 +177,32 @@ bool Flag_Region_CCSN( const int i, const int j, const int k, const int lv, cons # else const double dR [3] = { Pos[0]-amr->BoxCenter[0], Pos[1]-amr->BoxCenter[1], Pos[2]-amr->BoxCenter[2] }; # endif - const double R = sqrt( SQR(dR[0]) + SQR(dR[1]) + SQR(dR[2]) ); + const double R = sqrt( SQR(dR[0]) + SQR(dR[1]) + SQR(dR[2]) ); + +// must check CCSN_MaxRefine_Rad before evaluating other criteria + if ( R * UNIT_L <= CCSN_MaxRefine_Rad ) + return true; // check the maximum allowed refinement level based on angular resolution if ( CCSN_AngRes_Max > 0.0 && 2.0 * R * CCSN_AngRes_Max > dh ) Within = false; + if ( !Within ) return Within; + + +// check allowed maximum refine level based on distance to the center + if ( CCSN_REF_RBase > 0.0 ) + { + const double R_base = CCSN_REF_RBase; + int ratio = (int) ( R / R_base ); + int dlv = 0; + + while ( ratio ) { dlv += 1; ratio >>= 1; } + + if ( lv + dlv >= MAX_LEVEL ) Within = false; + } + return Within; diff --git a/src/TestProblem/Hydro/CCSN/Init_TestProb_Hydro_CCSN.cpp b/src/TestProblem/Hydro/CCSN/Init_TestProb_Hydro_CCSN.cpp index 069d77792..4f91b47c4 100644 --- a/src/TestProblem/Hydro/CCSN/Init_TestProb_Hydro_CCSN.cpp +++ b/src/TestProblem/Hydro/CCSN/Init_TestProb_Hydro_CCSN.cpp @@ -56,13 +56,15 @@ static int CCSN_Eint_Mode; // Mode of obtaining internal double CCSN_MaxRefine_Rad; // radius in cm within which to refine to the maximum allowed level double CCSN_CC_CentralDensFac; // factor that reduces the dt constrained by the central density (in cgs) during the core collapse double CCSN_CC_Red_DT; // reduced time step (in s) when the central density exceeds CCSN_CC_CentralDensFac before bounce - double CCSN_LB_TimeFac; // factor that scales the dt constrained by lightbulb scheme + double CCSN_NuHeat_TimeFac; // factor that scales the dt constrained by the lightbulb/leakage scheme int CCSN_CC_Rot; // mode for rotational profile (0:off, 1:analytical, 2:table) // --> analytical formula: Omega(r)=Omega_0*[R_0^2/(r^2+R_0^2)], where r is the spherical radius double CCSN_CC_Rot_R0; // characteristic radius R_0 (in cm) in the analytical rotational profile double CCSN_CC_Rot_Omega0; // central angular frequency Omega_0 (in rad/s) in the analytical rotational profile double CCSN_CC_Rot_Fac; // multiplication factor for the tabular rotational profile + double CCSN_REF_RBase; // reference distance for determining a maximum refinement level based on distance from the box center (in cm) + bool CCSN_Is_PostBounce = false; // boolean that indicates whether core bounce has occurred double CCSN_AngRes_Min; // minimum angular resolution in degree @@ -70,19 +72,22 @@ static int CCSN_Eint_Mode; // Mode of obtaining internal double CCSN_Shock_ThresFac_Pres; // pressure threshold factor for detecting postbounce shock double CCSN_Shock_ThresFac_Vel; // velocity threshold facotr for detecting postbounce shock int CCSN_Shock_Weight; // weighting of each cell for detecting postbounce shock (1:volume, 2:1/volume) + + int CCSN_DT_YE; // dt criterion on Ye (1=Ye-Ye_min/Ye_max, 2=Ye, 3=none) [1] // ======================================================================================= // problem-specific function prototypes void Record_CCSN_CentralQuant(); +void Record_CCSN_Leakage(); void Record_CCSN_GWSignal(); void Detect_CoreBounce(); void Detect_Shock(); -double Mis_GetTimeStep_Lightbulb( const int lv, const double dTime_dt ); double Mis_GetTimeStep_CoreCollapse( const int lv, const double dTime_dt ); +double Mis_GetTimeStep_PostBounce( const int lv, const double dTime_dt ); bool Flag_Region_CCSN( const int i, const int j, const int k, const int lv, const int PID ); bool Flag_CoreCollapse( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold ); -bool Flag_Lightbulb( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold ); +bool Flag_PostBounce( const int i, const int j, const int k, const int lv, const int PID, const double *Threshold ); @@ -178,18 +183,20 @@ void LoadInputTestProb( const LoadParaMode_t load_mode, ReadPara_t *ReadPara, HD LOAD_PARA( load_mode, "CCSN_CC_MaxRefine_Dens2", &CCSN_CC_MaxRefine_Dens2, 1.0e12, 0.0, NoMax_double ); LOAD_PARA( load_mode, "CCSN_CC_CentralDensFac", &CCSN_CC_CentralDensFac, 1.0e13, Eps_double, NoMax_double ); LOAD_PARA( load_mode, "CCSN_CC_Red_DT", &CCSN_CC_Red_DT, 1.0e-5, Eps_double, NoMax_double ); - LOAD_PARA( load_mode, "CCSN_LB_TimeFac", &CCSN_LB_TimeFac, 0.1, Eps_double, 1.0 ); + LOAD_PARA( load_mode, "CCSN_NuHeat_TimeFac", &CCSN_NuHeat_TimeFac, 0.1, Eps_double, 1.0 ); LOAD_PARA( load_mode, "CCSN_CC_Rot", &CCSN_CC_Rot, 2, 0, 2 ); LOAD_PARA( load_mode, "CCSN_CC_Rot_R0", &CCSN_CC_Rot_R0, 2.0e8, Eps_double, NoMax_double ); LOAD_PARA( load_mode, "CCSN_CC_Rot_Omega0", &CCSN_CC_Rot_Omega0, 0.5, 0.0, NoMax_double ); LOAD_PARA( load_mode, "CCSN_CC_Rot_Fac", &CCSN_CC_Rot_Fac, -1.0, NoMin_double, NoMax_double ); + LOAD_PARA( load_mode, "CCSN_REF_RBase", &CCSN_REF_RBase, 1.25e7, NoMin_double, NoMax_double ); LOAD_PARA( load_mode, "CCSN_Is_PostBounce", &CCSN_Is_PostBounce, false, Useless_bool, Useless_bool ); LOAD_PARA( load_mode, "CCSN_MaxRefine_Rad", &CCSN_MaxRefine_Rad, 3.0e6, Eps_double, NoMax_double ); LOAD_PARA( load_mode, "CCSN_AngRes_Min", &CCSN_AngRes_Min, -1.0, NoMin_double, NoMax_double ); LOAD_PARA( load_mode, "CCSN_AngRes_Max", &CCSN_AngRes_Max, -1.0, NoMin_double, NoMax_double ); LOAD_PARA( load_mode, "CCSN_Shock_ThresFac_Pres", &CCSN_Shock_ThresFac_Pres, 0.5, Eps_double, NoMax_double ); - LOAD_PARA( load_mode, "CCSN_Shock_ThresFac_Vel" , &CCSN_Shock_ThresFac_Vel, 0.1, Eps_double, NoMax_double ); - LOAD_PARA( load_mode, "CCSN_Shock_Weight" , &CCSN_Shock_Weight, 2, 1, 2 ); + LOAD_PARA( load_mode, "CCSN_Shock_ThresFac_Vel", &CCSN_Shock_ThresFac_Vel, 0.1, Eps_double, NoMax_double ); + LOAD_PARA( load_mode, "CCSN_Shock_Weight", &CCSN_Shock_Weight, 2, 1, 2 ); + LOAD_PARA( load_mode, "CCSN_DT_YE", &CCSN_DT_YE, 1, 1, 3 ); } // FUNCITON : LoadInputTestProb @@ -325,7 +332,7 @@ void SetParameter() if ( CCSN_AngRes_Max > 0.0 ) { CCSN_AngRes_Max *= Deg2Rad; - PRINT_RESET_PARA( CCSN_AngRes_Max, FORMAT_DOUBLE, "" ); + PRINT_RESET_PARA( CCSN_AngRes_Max, FORMAT_REAL, "" ); if ( !OPT__FLAG_REGION ) Aux_Error( ERROR_INFO, "%s is disabled for %s = %13.7e !!\n", "OPT__FLAG_REGION", "CCSN_AngRes_Max", CCSN_AngRes_Max ); @@ -333,7 +340,7 @@ void SetParameter() if ( CCSN_AngRes_Min > 0.0 ) { CCSN_AngRes_Min *= Deg2Rad; - PRINT_RESET_PARA( CCSN_AngRes_Min, FORMAT_DOUBLE, "" ); + PRINT_RESET_PARA( CCSN_AngRes_Min, FORMAT_REAL, "" ); } if ( CCSN_AngRes_Min > 0.0 && CCSN_AngRes_Max > 0.0 && @@ -342,6 +349,8 @@ void SetParameter() // (2) set the problem-specific derived parameters +// convert runtime parameters to the code unit + CCSN_REF_RBase /= UNIT_L; // (3) reset other general-purpose parameters @@ -377,7 +386,7 @@ void SetParameter() Aux_Message( stdout, " sampling interval of GW signals = %13.7e\n", CCSN_GW_DT ); Aux_Message( stdout, " mode for obtaining internal energy = %d\n", CCSN_Eint_Mode ); if ( CCSN_Prob != Migration_Test ) { - Aux_Message( stdout, " scaling factor for lightbulb dt = %13.7e\n", CCSN_LB_TimeFac ); + Aux_Message( stdout, " scaling factor for lightbulb/leakage dt = %13.7e\n", CCSN_NuHeat_TimeFac ); Aux_Message( stdout, " has core bounce occurred = %d\n", CCSN_Is_PostBounce ); Aux_Message( stdout, " pressure threshold factor for detecting shock = %13.7e\n", CCSN_Shock_ThresFac_Pres ); Aux_Message( stdout, " velocity threshold factor for detecting shock = %13.7e\n", CCSN_Shock_ThresFac_Vel ); @@ -400,6 +409,8 @@ void SetParameter() Aux_Message( stdout, " radius within which to refine to the maximum allowed level (in cm) = %13.7e\n", CCSN_MaxRefine_Rad ); Aux_Message( stdout, " minimum angular resolution (in degrees) = %13.7e\n", CCSN_AngRes_Min/Deg2Rad ); Aux_Message( stdout, " maximum angular resolution (in degrees) = %13.7e\n", CCSN_AngRes_Max/Deg2Rad ); + Aux_Message( stdout, " reference distance for the maximum refinement level = %13.7e\n", CCSN_REF_RBase ); + Aux_Message( stdout, " dt criterion on Ye = %d\n", CCSN_DT_YE ); Aux_Message( stdout, "=======================================================================================\n" ); } @@ -518,14 +529,17 @@ void SetGridIC( real fluid[], const double x, const double y, const double z, co # if ( EOS == EOS_NUCLEAR ) real *Passive = new real [NCOMP_PASSIVE]; - Passive[ YE - NCOMP_FLUID ] = Ye*Dens; - Passive[ DEDT_NU - NCOMP_FLUID ] = TINY_NUMBER; + Passive[ YE - NCOMP_FLUID ] = Ye*Dens; + Passive[ DEDT_NU - NCOMP_FLUID ] = 0.0; +# ifdef DYEDT_NU + Passive[ DYEDT_NU - NCOMP_FLUID ] = 0.0; +# endif # ifdef TEMP_IG - Passive[ TEMP_IG - NCOMP_FLUID ] = Temp; + Passive[ TEMP_IG - NCOMP_FLUID ] = Temp; # endif # else real *Passive = NULL; -# endif +# endif // #if ( EOS == EOS_NUCLEAR ) ... else ... if ( CCSN_Eint_Mode == 1 ) // Temperature Mode { @@ -810,16 +824,62 @@ void Load_IC_Prof_CCSN() void Record_CCSN() { -// (1) shock detection +// (1) check whether the core bounce occurs + if ( !CCSN_Is_PostBounce ) + { + Detect_CoreBounce(); + + if ( CCSN_Is_PostBounce ) + { +// dump the bounce time in standard output + if ( MPI_Rank == 0 ) Aux_Message( stdout, "Bounce time = %13.7e seconds !!\n", Time[0] * UNIT_T ); + +// disable the deleptonization scheme, and enable the lightbulb/leakage scheme + SrcTerms.Deleptonization = false; + +# ifdef NEUTRINO_SCHEME +# if ( NEUTRINO_SCHEME == LIGHTBULB ) + SrcTerms.Lightbulb = true; + if ( MPI_Rank == 0 ) Aux_Message( stdout, "Enable the lightbulb scheme !!\n" ); +# elif ( NEUTRINO_SCHEME == LEAKAGE ) + SrcTerms.Leakage = true; + if ( MPI_Rank == 0 ) Aux_Message( stdout, "Enable the leakage scheme !!\n" ); +# else + if ( MPI_Rank == 0 ) Aux_Message( stdout, "No NEUTRINO_SCHEME specified !!\n" ); +# endif +# endif + + Src_Init(); + +// initialize the dEdt_Nu field + for (int lv=0; lvFluSg[lv], amr->MagSg[lv], false, false ); + + Buf_GetBufferData( lv, amr->FluSg[lv], amr->MagSg[lv], NULL_INT, DATA_GENERAL, _TOTAL, _MAG, Flu_ParaBuf, USELB_YES ); + } + +// forced output data at core bounce + Output_DumpData( 2 ); + + } + } // if ( !CCSN_Is_PostBounce ) + + +// (2) shock detection if ( CCSN_Prob != Migration_Test && CCSN_Is_PostBounce ) Detect_Shock(); -// (2) record quantities at the center +// (3) record quantities at the center +# if ( defined NEUTRINO_SCHEME && NEUTRINO_SCHEME == LEAKAGE ) + if ( CCSN_Is_PostBounce ) Record_CCSN_Leakage(); +# endif + Record_CCSN_CentralQuant(); -// (3) GW signal +// (4) GW signal # ifdef GRAVITY if ( CCSN_GW_OUTPUT ) { @@ -852,29 +912,6 @@ void Record_CCSN() } // if ( CCSN_GW_OUTPUT ) # endif - -// (4) check whether the core bounce occurs - if ( !CCSN_Is_PostBounce ) - { - Detect_CoreBounce(); - - if ( CCSN_Is_PostBounce ) - { -// dump the bounce time in standard output - if ( MPI_Rank == 0 ) Aux_Message( stdout, "Bounce time = %13.7e seconds !!\n", Time[0] * UNIT_T ); - -// disable the deleptonization scheme, and enable the lightbulb scheme - SrcTerms.Deleptonization = false; - SrcTerms.Lightbulb = true; - - Src_Init(); - -// forced output data at core bounce - Output_DumpData( 2 ); - - } - } // if ( !CCSN_Is_PostBounce ) - } // FUNCTION : Record_CCSN() @@ -882,21 +919,26 @@ void Record_CCSN() //------------------------------------------------------------------------------------------------------- // Function : Mis_GetTimeStep_CCSN // Description : Interface for invoking several functions for estimating the evolution time-step +// for CCSN simulations //------------------------------------------------------------------------------------------------------- double Mis_GetTimeStep_CCSN( const int lv, const double dTime_dt ) { double dt_CCSN = HUGE_NUMBER; - if ( SrcTerms.Lightbulb ) + if ( ! CCSN_Is_PostBounce ) { - const double dt_LB = Mis_GetTimeStep_Lightbulb( lv, dTime_dt ); - - dt_CCSN = fmin( dt_CCSN, dt_LB ); + if ( SrcTerms.Deleptonization ) + dt_CCSN = fmin( dt_CCSN, Mis_GetTimeStep_CoreCollapse( lv, dTime_dt ) ); } - if ( !CCSN_Is_PostBounce && SrcTerms.Deleptonization ) - dt_CCSN = fmin( dt_CCSN, Mis_GetTimeStep_CoreCollapse( lv, dTime_dt ) ); + else + { +# if ( NEUTRINO_SCHEME == LIGHTBULB || NEUTRINO_SCHEME == LEAKAGE ) + if ( SrcTerms.Lightbulb || SrcTerms.Leakage ) + dt_CCSN = fmin( dt_CCSN, Mis_GetTimeStep_PostBounce ( lv, dTime_dt ) ); +# endif + } return dt_CCSN; @@ -935,9 +977,10 @@ bool Flag_CCSN( const int i, const int j, const int k, const int lv, const int P if ( Flag ) return Flag; } - if ( ( CCSN_Prob == Post_Bounce ) || SrcTerms.Lightbulb ) + if ( CCSN_Is_PostBounce && + ( SrcTerms.Lightbulb || SrcTerms.Leakage ) ) { - Flag |= Flag_Lightbulb( i, j, k, lv, PID, Threshold ); + Flag |= Flag_PostBounce( i, j, k, lv, PID, Threshold ); if ( Flag ) return Flag; } @@ -994,11 +1037,9 @@ void Init_TestProb_Hydro_CCSN() // set the function pointers of various problem-specific routines Init_Function_User_Ptr = SetGridIC; - Flag_User_Ptr = Flag_CCSN; Flag_Region_Ptr = Flag_Region_CCSN; Aux_Record_User_Ptr = Record_CCSN; End_User_Ptr = End_CCSN; - Mis_GetTimeStep_User_Ptr = Mis_GetTimeStep_CCSN; # if ( EOS == EOS_NUCLEAR && NUC_TABLE_MODE == NUC_TABLE_MODE_TEMP ) Flu_ResetByUser_Func_Ptr = Flu_ResetByUser_CCSN; @@ -1008,6 +1049,12 @@ void Init_TestProb_Hydro_CCSN() Output_HDF5_InputTest_Ptr = LoadInputTestProb; # endif + if ( CCSN_Prob != Migration_Test ) + { + Flag_User_Ptr = Flag_CCSN; + Mis_GetTimeStep_User_Ptr = Mis_GetTimeStep_CCSN; + } + # ifdef MHD switch ( CCSN_Mag ) { diff --git a/src/TestProblem/Hydro/CCSN/Record_CCSN.cpp b/src/TestProblem/Hydro/CCSN/Record_CCSN.cpp index a33b1dded..a1d7a3b70 100644 --- a/src/TestProblem/Hydro/CCSN/Record_CCSN.cpp +++ b/src/TestProblem/Hydro/CCSN/Record_CCSN.cpp @@ -2,14 +2,28 @@ double CCSN_CentralDens; - double CCSN_Rsh_Min = 0.0; - double CCSN_Rsh_Max = 0.0; - double CCSN_Rsh_Ave = 0.0; + + double CCSN_Rsh_Min = 0.0; + double CCSN_Rsh_Max = 0.0; + double CCSN_Rsh_Ave = 0.0; + double CCSN_Rsh_Ave_V = 0.0; + double CCSN_Rsh_Ave_Vinv = 0.0; + + double CCSN_Leakage_NetHeatGain = 0.0; + double CCSN_Leakage_Lum [3] = { 0.0 }; + double CCSN_Leakage_Heat [3] = { 0.0 }; + double CCSN_Leakage_NetHeat[3] = { 0.0 }; + double CCSN_Leakage_EAve [3] = { 0.0 }; + double CCSN_Leakage_RadNS [3] = { 0.0 }; + extern bool CCSN_Is_PostBounce; extern double CCSN_Shock_ThresFac_Pres; extern double CCSN_Shock_ThresFac_Vel; extern int CCSN_Shock_Weight; +extern void Src_WorkBeforeMajorFunc_Leakage( const int lv, const double TimeNew, const double TimeOld, const double dt, + double AuxArray_Flt[], int AuxArray_Int[] ); + //------------------------------------------------------------------------------------------------------- @@ -123,30 +137,70 @@ void Record_CCSN_CentralQuant() // write to the file "Record__CentralQuant" by the MPI process which has the target patch +# if ( NEUTRINO_SCHEME == LEAKAGE ) + const int NColumn = 28; +# else + const int NColumn = 14; +# endif + if ( MPI_Rank == Data_Int[0] ) { static bool FirstTime = true; - char filename_central_quant[2*MAX_STRING]; - sprintf( filename_central_quant, "%s/Record__CentralQuant", OUTPUT_DIR ); + char FileName[2*MAX_STRING]; + sprintf( FileName, "%s/Record__CentralQuant", OUTPUT_DIR ); // file header if ( FirstTime ) { - if ( Aux_CheckFileExist(filename_central_quant) ) - { - Aux_Message( stderr, "WARNING : file \"%s\" already exists !!\n", filename_central_quant ); - } + if ( Aux_CheckFileExist(FileName) ) + Aux_Message( stderr, "WARNING : file \"%s\" already exists !!\n", FileName ); else { - FILE *file_cent_quant = fopen( filename_central_quant, "w" ); - fprintf( file_cent_quant, "#%14s %12s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s %16s\n", - "1_Time [sec]", "2_Step", "3_PosX [cm]", "4_PosY [cm]", "5_PosZ [cm]", - "6_Dens [g/cm^3]", "7_Ye", "8_Rsh_Min [cm]", "9_Rsh_Ave [cm]", "10_Rsh_Max [cm]", - "11_GREP_PosX [cm]", "12_GREP_PosY [cm]", "13_GREP_PosZ [cm]" ); - fclose( file_cent_quant ); + FILE *File = fopen( FileName, "w" ); + +// column index + Aux_Message( File, "#%14s %8s", "[ 1]", "[ 2]" ); + for (int c=2; cFluSgTime[lv][ amr->FluSg[lv] ]; + + Src_WorkBeforeMajorFunc_Leakage( lv, TimeNew, TimeNew, 0.0, + Src_Leakage_AuxArray_Flt, Src_Leakage_AuxArray_Int ); + +// enable the record mode + Src_Leakage_AuxArray_Int[3] = LEAK_MODE_RECORD; + +// allocate memory for per-thread arrays +# ifdef OPENMP + const int NT = OMP_NTHREAD; +# else + const int NT = 1; +# endif + + double NetHeatGain = 0.0; // net heating rate in gain region + double Lum [NType_Neutrino] = { 0.0 }; + double Heat [NType_Neutrino] = { 0.0 }; + double NetHeat[NType_Neutrino] = { 0.0 }; + + double *OMP_NetHeatGain = new double [NT]; + double **OMP_Lum, **OMP_Heat, **OMP_NetHeat; + + Aux_AllocateArray2D( OMP_Lum, NT, NType_Neutrino ); + Aux_AllocateArray2D( OMP_Heat, NT, NType_Neutrino ); + Aux_AllocateArray2D( OMP_NetHeat, NT, NType_Neutrino ); + + +# pragma omp parallel + { +# ifdef OPENMP + const int TID = omp_get_thread_num(); +# else + const int TID = 0; +# endif + +// initialize arrays + OMP_NetHeatGain[TID] = 0.0; + + for (int n=0; ndh[lv]; + const double dv_CGS = CUBE( dh * UNIT_L ); + + +# pragma omp for schedule( static ) + for (int PID=0; PIDNPatchComma[lv][1]; PID++) + { +// skip non-leaf patches + if ( amr->patch[0][lv][PID]->son != -1 ) continue; + + const double z0 = amr->patch[0][lv][PID]->EdgeL[2]; + const double y0 = amr->patch[0][lv][PID]->EdgeL[1]; + const double x0 = amr->patch[0][lv][PID]->EdgeL[0]; + + for (int k=0; kpatch[ amr->FluSg[lv] ][lv][PID]->fluid[v][k][j][i]; + +# ifdef MHD + real B[NCOMP_MAG]; + + MHD_GetCellCenteredBFieldInPatch( B, lv, PID, i, j, k, amr->MagSg[lv] ); +# else + real *B = NULL; +# endif + + SrcTerms.Leakage_CPUPtr( fluid, B, &SrcTerms, 0.0, NULL_REAL, x, y, z, NULL_REAL, NULL_REAL, + MIN_DENS, MIN_PRES, MIN_EINT, PassiveFloorMask, &EoS, + Src_Leakage_AuxArray_Flt, Src_Leakage_AuxArray_Int ); + + if ( fluid[DENS] > 0.0 ) + OMP_NetHeatGain[TID] += fluid[DENS] * dv_CGS; + + OMP_Lum [TID][0] += fluid[MOMX] * dv_CGS; + OMP_Lum [TID][1] += fluid[MOMY] * dv_CGS; + OMP_Lum [TID][2] += fluid[MOMZ] * dv_CGS; + OMP_Heat [TID][0] += fluid[ENGY] * dv_CGS; +# ifdef YE + OMP_Heat [TID][1] += fluid[YE ] * dv_CGS; +# endif +# ifdef DYEDT_NU + OMP_NetHeat[TID][0] += fluid[DEDT_NU ] * dv_CGS; + OMP_NetHeat[TID][1] += fluid[DYEDT_NU] * dv_CGS; +# endif + + }}} // i,j,k + } // for (int PID=0; PIDNPatchComma[lv][1]; PID++) + } // for (int lv=0; lv<=MAX_LEVEL; lv++) + } // OpenMP parallel region + + +// sum over all OpenMP threads + for (int t=0; tBoxCenter[0], amr->BoxCenter[1], amr->BoxCenter[2] }; + +# ifndef MHD + const int OPT__MAG_INT_SCHEME = INT_NONE; +# endif + # ifdef OPENMP const int NT = OMP_NTHREAD; # else const int NT = 1; # endif - const int NGhost = 1; - const int PS1P2 = PS1 + 2*NGhost; - const int NPG_Max = FLU_GPU_NPGROUP; - const int TVarCC = _DENS | _PASSIVE | _VELX | _VELY | _VELZ | _PRES; - const int Stride1 = CUBE( PS1P2 ); - const int Stride2 = NCOMP_TOTAL * Stride1; - const int VELX = NCOMP_PASSIVE + 1; - const int VELY = NCOMP_PASSIVE + 2; - const int VELZ = NCOMP_PASSIVE + 3; - const int PRES = NCOMP_PASSIVE + 4; - double Center[3]; - double Shock_Min = HUGE_NUMBER; - double Shock_Max = -HUGE_NUMBER; - double Shock_Ave = 0.0; - double Shock_Weight = 0.0; - int Shock_Found = false; - - double OMP_Shock_Min [NT]; - double OMP_Shock_Max [NT]; - double OMP_Shock_Ave [NT]; - double OMP_Shock_Weight[NT]; - int OMP_Shock_Found [NT]; - - real *OMP_Fluid = new real [ 8*NPG_Max*Stride2 ]; - real *OMP_Pres_Min = new real [ 8*NPG_Max*Stride1 ]; - real *OMP_Cs_Min = new real [ 8*NPG_Max*Stride1 ]; + double Shock_Min = HUGE_NUMBER; + double Shock_Max = -HUGE_NUMBER; + double Shock_Ave_V = 0.0; + double Shock_Ave_Vinv = 0.0; + double Shock_Weight_V = 0.0; + double Shock_Weight_Vinv = 0.0; + int Shock_Found = false; + + double OMP_Shock_Min [NT]; + double OMP_Shock_Max [NT]; + double OMP_Shock_Ave_V [NT]; + double OMP_Shock_Ave_Vinv [NT]; + double OMP_Shock_Weight_V [NT]; + double OMP_Shock_Weight_Vinv[NT]; + int OMP_Shock_Found [NT]; + + real *OMP_Fluid = new real [ 8*NPG_Max*NCOMP_TOTAL*CUBE(SHK_NXT) ]; + real *OMP_Pres_Min = new real [ 8*NPG_Max* CUBE(SHK_NXT) ]; + real *OMP_Cs_Min = new real [ 8*NPG_Max* CUBE(SHK_NXT) ]; +# ifdef MHD + real *OMP_Mag = new real [ 8*NPG_Max*NCOMP_MAG * SQR(SHK_NXT)*SHK_NXT_P1 ]; +# else + real *OMP_Mag = NULL; +# endif + for (int t=0; tdh[lv]; const double dv = CUBE( dh ); + const double _dv = 1.0 / dv; const int NTotal = amr->NPatchComma[lv][1] / 8; int *PID0_List = new int [NTotal]; +// obsolete!! to be removed in the future. double Weight; switch ( CCSN_Shock_Weight ) { @@ -562,11 +812,11 @@ void Detect_Shock() { const int NPG = MIN( NPG_Max, NTotal-Disp ); -// (1-a) prepare the primitive and passive variables -// note that the data are prepared in order of density, Passive, velx, vely, velz, and pressure - Prepare_PatchData( lv, Time[lv], OMP_Fluid, NULL, - NGhost, NPG, PID0_List+Disp, TVarCC, _NONE, - OPT__FLU_INT_SCHEME, INT_NONE, UNIT_PATCH, NSIDE_26, false, +// (1) prepare the data +// --> prepare all fields to ensure monotonicity + Prepare_PatchData( lv, Time[lv], OMP_Fluid, OMP_Mag, + SHK_GHOST_SIZE, NPG, PID0_List+Disp, _TOTAL, _MAG, + OPT__FLU_INT_SCHEME, OPT__MAG_INT_SCHEME, UNIT_PATCH, NSIDE_26, false, OPT__BC_FLU, BC_POT_NONE, -1.0, -1.0, -1.0, -1.0, false ); @@ -583,41 +833,66 @@ void Detect_Shock() if ( amr->patch[0][lv][PID]->son != -1 ) continue; - real (*Fluid )[PS1P2][PS1P2][PS1P2] = ( real(*)[PS1P2][PS1P2][PS1P2] ) ( OMP_Fluid + PID_IDX*Stride2 ); - real (*Pres_Min) [PS1P2][PS1P2] = ( real(*) [PS1P2][PS1P2] ) ( OMP_Pres_Min + PID_IDX*Stride1 ); - real (*Cs_Min ) [PS1P2][PS1P2] = ( real(*) [PS1P2][PS1P2] ) ( OMP_Cs_Min + PID_IDX*Stride1 ); -// (1-b) compute the sound speed and store in the density field - for (int k=0; k the sound speed and pressure are stored in the density and energy fields, respectively + for (int k=0; k=(real)0.0), MIN_PRES, PassiveFloorMask, Emag, + EoS_DensEint2Pres_CPUPtr, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, + EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table, NULL ); + + const real CSqr = EoS_DensPres2CSqr_CPUPtr( FluidForEoS[DENS], Pres, FluidForEoS+NCOMP_FLUID, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); - Fluid[DENS][k][j][i] = SQRT( CSqr ); +// store pressure and sound speed, and convert momentum to velocity + Fluid[DENS][k][j][i] = SQRT( CSqr ); + Fluid[MOMX][k][j][i] /= FluidForEoS[DENS]; + Fluid[MOMY][k][j][i] /= FluidForEoS[DENS]; + Fluid[MOMZ][k][j][i] /= FluidForEoS[DENS]; + Fluid[ENGY][k][j][i] = Pres; }}} // i, j, k -// (2) use a naive method to find the minimum of pressure and sound speed in the local 3x3 subarray - for (int k=NGhost; kpatch[0][lv][PID]->EdgeL[2] + (k+0.5)*dh; const int kk = k+NGhost; - for (int j=0; jpatch[0][lv][PID]->EdgeL[1] + (j+0.5)*dh; const int jj = j+NGhost; - for (int i=0; ipatch[0][lv][PID]->EdgeL[0] + (i+0.5)*dh; const int ii = i+NGhost; + for (int k=0; kpatch[0][lv][PID]->EdgeL[2] + (k+0.5)*dh; const int kk = k+SHK_GHOST_SIZE; + for (int j=0; jpatch[0][lv][PID]->EdgeL[1] + (j+0.5)*dh; const int jj = j+SHK_GHOST_SIZE; + for (int i=0; ipatch[0][lv][PID]->EdgeL[0] + (i+0.5)*dh; const int ii = i+SHK_GHOST_SIZE; const double dx = x - Center[0]; const double dy = y - Center[1]; @@ -637,23 +912,25 @@ void Detect_Shock() const double r = sqrt( SQR( dx ) + SQR( dy ) + SQR( dz ) ); // (3) evaluate the undivided gradient of pressure and the undivided divergence of velocity - real GradP = (real)0.5 * ( FABS( Fluid[PRES][kk+1][jj ][ii ] - Fluid[PRES][kk-1][jj ][ii ] ) - + FABS( Fluid[PRES][kk ][jj+1][ii ] - Fluid[PRES][kk ][jj-1][ii ] ) - + FABS( Fluid[PRES][kk ][jj ][ii+1] - Fluid[PRES][kk ][jj ][ii-1] ) ); + real GradP = (real)0.5 * ( FABS( Fluid[ENGY][kk+1][jj ][ii ] - Fluid[ENGY][kk-1][jj ][ii ] ) + + FABS( Fluid[ENGY][kk ][jj+1][ii ] - Fluid[ENGY][kk ][jj-1][ii ] ) + + FABS( Fluid[ENGY][kk ][jj ][ii+1] - Fluid[ENGY][kk ][jj ][ii-1] ) ); - real DivV = (real)0.5 * ( ( Fluid[VELZ][kk+1][jj ][ii ] - Fluid[VELZ][kk-1][jj ][ii ] ) - + ( Fluid[VELY][kk ][jj+1][ii ] - Fluid[VELY][kk ][jj-1][ii ] ) - + ( Fluid[VELX][kk ][jj ][ii+1] - Fluid[VELX][kk ][jj ][ii-1] ) ); + real DivV = (real)0.5 * ( ( Fluid[MOMZ][kk+1][jj ][ii ] - Fluid[MOMZ][kk-1][jj ][ii ] ) + + ( Fluid[MOMY][kk ][jj+1][ii ] - Fluid[MOMY][kk ][jj-1][ii ] ) + + ( Fluid[MOMX][kk ][jj ][ii+1] - Fluid[MOMX][kk ][jj ][ii-1] ) ); // (4) examine the criteria for detecting strong shock if ( ( GradP >= CCSN_Shock_ThresFac_Pres * Pres_Min[kk][jj][ii] ) && ( DivV <= -CCSN_Shock_ThresFac_Vel * Cs_Min [kk][jj][ii] ) ) { - OMP_Shock_Min [TID] = fmin( r, OMP_Shock_Min[TID] ); - OMP_Shock_Max [TID] = fmax( r, OMP_Shock_Max[TID] ); - OMP_Shock_Ave [TID] += r * Weight; - OMP_Shock_Weight[TID] += Weight; - OMP_Shock_Found [TID] = true; + OMP_Shock_Min [TID] = fmin( r, OMP_Shock_Min[TID] ); + OMP_Shock_Max [TID] = fmax( r, OMP_Shock_Max[TID] ); + OMP_Shock_Ave_V [TID] += r * dv; + OMP_Shock_Ave_Vinv [TID] += r * _dv; + OMP_Shock_Weight_V [TID] += dv; + OMP_Shock_Weight_Vinv[TID] += _dv; + OMP_Shock_Found [TID] = true; } }}} // i,j,k @@ -668,32 +945,41 @@ void Detect_Shock() delete [] OMP_Fluid; delete [] OMP_Pres_Min; delete [] OMP_Cs_Min; +# ifdef MHD + delete [] OMP_Mag; +# endif // collect data over all OpenMP threads for (int t=0; tdh[lv]; # pragma omp for schedule( runtime ) for (int PID=0; PIDNPatchComma[lv][1]; PID++) { + + if ( amr->patch[0][lv][PID]->son != -1 ) continue; + for (int k=0; kpatch[ amr->FluSg[lv] ][lv][PID]->fluid[MOMY][k][j][i]; const real Momz = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid[MOMZ][k][j][i]; const real Engy = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid[ENGY][k][j][i]; +# if ( defined YE && NEUTRINO_SCHEME == LEAKAGE ) + const real Ye = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid[YE ][k][j][i] / Dens; +# else + const real Ye = NULL_REAL; +# endif # ifdef MHD - real B[NCOMP_MAG]; + real B[NCOMP_MAG]; MHD_GetCellCenteredBFieldInPatch( B, lv, PID, i, j, k, amr->MagSg[lv] ); @@ -84,47 +95,40 @@ double Mis_GetTimeStep_Lightbulb( const int lv, const double dTime_dt ) const real Emag = NULL_REAL; # endif // ifdef MHD ... else ... - const real Eint_Code = Hydro_Con2Eint( Dens, Momx, Momy, Momz, Engy, true, MIN_EINT, PassiveFloorMask, Emag, + const real Eint_Code = Hydro_Con2Eint( Dens, Momx, Momy, Momz, Engy, false, MIN_EINT, PassiveFloorMask, Emag, EoS_GuessHTilde_CPUPtr, EoS_HTilde2Temp_CPUPtr, EoS_AuxArray_Flt, EoS_AuxArray_Int, h_EoS_Table ); - real dEint_Code = NULL_REAL; +# ifdef DEDT_NU + const real dEint_Code = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid[DEDT_NU][k][j][i]; +# else + const real dEint_Code = NULL_REAL; +# endif // ifdef DEDT_NU ... else ... - if ( IsInit_dEdt_Nu ) - { -// use the stored neutrino heating/cooling rate -# ifdef DEDT_NU - dEint_Code = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid[DEDT_NU][k][j][i]; -# endif - } +# ifdef DYEDT_NU + const real dYedt = amr->patch[ amr->FluSg[lv] ][lv][PID]->fluid[DYEDT_NU][k][j][i]; + real dYe; - else + switch ( CCSN_DT_YE ) { -// call Src_Lightbulb() to compute the neutrino heating/cooling rate if not initialized yet - const double z = amr->patch[0][lv][PID]->EdgeL[2] + (k+0.5)*dh; - const double y = amr->patch[0][lv][PID]->EdgeL[1] + (j+0.5)*dh; - const double x = amr->patch[0][lv][PID]->EdgeL[0] + (i+0.5)*dh; - -// get the input arrays - real fluid[FLU_NIN_S]; - - for (int v=0; vpatch[ amr->FluSg[lv] ][lv][PID]->fluid[v][k][j][i]; - - - SrcTerms.Lightbulb_CPUPtr( fluid, B, &SrcTerms, 0.0, NULL_REAL, x, y, z, NULL_REAL, NULL_REAL, - MIN_DENS, MIN_PRES, MIN_EINT, PassiveFloorMask, &EoS, - Src_Lightbulb_AuxArray_Flt, Src_Lightbulb_AuxArray_Int ); + case 1 : dYe = ( dYedt > 0.0 ) ? ( YeMax - Ye ) : ( YeMin - Ye ); break; + case 2 : dYe = Ye; break; + case 3 : dYe = __FLT_MAX__; break; + } +# else + const real dYedt = NULL_REAL; + real dYe = NULL_REAL; +# endif // ifdef DYEDT_NU ... else ... -# ifdef DEDT_NU - dEint_Code = fluid[DEDT_NU]; -# endif - } // if ( IsInit_dEdt_Nu ) ... else ... + double dtInv_NuHeat_ThisCell = FABS( dEint_Code / Eint_Code ); - const double dt_LB_Inv_ThisCell = FABS( dEint_Code / Eint_Code ); +# if ( NEUTRINO_SCHEME == LEAKAGE ) + dtInv_NuHeat_ThisCell = FMAX( dYedt / dYe, dtInv_NuHeat_ThisCell ); +# endif // compare the inverse of ratio to avoid zero division, and store the maximum value - OMP_dt_LB_Inv[TID] = FMAX( OMP_dt_LB_Inv[TID], dt_LB_Inv_ThisCell ); + OMP_dtInv_NuHeat[TID] = FMAX( OMP_dtInv_NuHeat[TID], dtInv_NuHeat_ThisCell ); }}} // i,j,k } // for (int PID=0; PIDNPatchComma[lv][1]; PID++) @@ -132,23 +136,23 @@ double Mis_GetTimeStep_Lightbulb( const int lv, const double dTime_dt ) // find the maximum over all OpenMP threads - for (int TID=0; TID