Support the leakage scheme for CCSN

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

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

example/test_problem/Hydro/CCSN/Input__TestProb.CoreCollapse

@@ -29,7 +29,7 @@ CCSN_CC_MaxRefine_Dens2   1.0e12                   # central density threshold t
 CCSN_MaxRefine_RadFac     0.15                     # factor that determines the maximum refinement level based on distance from the box center  [0.15]
 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

example/test_problem/Hydro/CCSN/Input__TestProb.Leakage

@@ -0,0 +1,40 @@
+# 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   0                                 # output GW signal (0=off, 1=on)                    [0]
+CCSN_GW_DT       1.0                               # 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_MaxRefine_RadFac     0.15                     # factor that determines the maximum refinement level based on distance from the box center  [0.15]
+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.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
+CCSN_CC_Rot_R0            1.03376e8                # characteristic radius R_0 (in cm) in the analytical rotational profile             [2.0e8]
+CCSN_CC_Rot_Omega0        2.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_Is_PostBounce        0                        # boolean that indicates whether core bounce has occurred  [0]

example/test_problem/Hydro/CCSN/Input__TestProb.Lightbulb

@@ -21,6 +21,6 @@ CCSN_Eint_Mode   1                                 # Mode of obtaining internal
                                                    #   2=Pres Mode: Eint(dens, pres, [Ye]) )
 
 CCSN_MaxRefine_RadFac  0.15                        # factor that determines the maximum refinement level based on distance from the box center  [0.15]
-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     1                           # boolean that indicates whether core bounce has occurred   [0]

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

example/test_problem/Hydro/CCSN/clean.sh

@@ -3,4 +3,4 @@ rm -f Record__Note Record__Timing Record__TimeStep Record__PatchCount Record__Du
       Diag* 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
+      GRACKLE_INFO Record__DivB Record__CentralQuant Record__QuadMom_2nd Leakage_Lum_*

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     ] );

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 );

include/Global.h

@@ -289,6 +289,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     ];
@@ -373,6 +375,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

include/Macro.h

@@ -94,6 +94,7 @@
 #define LIGHTBULB    1
 #define IDSA         2
 #define M1           3
+#define LEAKAGE      4
 
 
 // Poisson solvers
@@ -165,13 +166,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
@@ -275,12 +285,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 )
@@ -323,16 +338,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 )
@@ -358,10 +378,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
 
@@ -395,10 +418,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
 
@@ -776,9 +802,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[]
 

include/Prototype.h

@@ -40,6 +40,10 @@ void Aux_ComputeProfile( Profile_t *Prof[], const double Center[], const double
                          const double PrepTime );
 void Aux_FindExtrema( Extrema_t *Extrema, const ExtremaMode_t Mode, const int MinLv, const int MaxLv,
                       const PatchType_t PatchType );
+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 PrepTime );
 #ifndef SERIAL
 void Aux_Record_BoundaryPatch( const int lv, int *NList, int **IDList, int **PosList );
 #endif

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,30 @@ 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;
+   real     *Leakage_Radius_DevPtr;
+   real     *Leakage_tau_DevPtr;
+   real     *Leakage_chi_DevPtr;
+   real     *Leakage_Heat_Flux_DevPtr;
+   real     *Leakage_HeatE_Rms_DevPtr;
+   real     *Leakage_HeatE_Ave_DevPtr;
+#  endif // if ( MODEL == HYDRO )
 
 // user-specified source term
    SrcFunc_t User_FuncPtr;

src/Auxiliary/Aux_Check_Parameter.cpp

@@ -754,8 +754,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 )
@@ -1477,8 +1477,44 @@ 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" );
+
+   if ( SrcTerms.Lightbulb  &&  SrcTerms.Leakage )
+      Aux_Error( ERROR_INFO, "SRC_LIGHTBULB and SRC_LEAKAGE cannot be enabled simultaneously !!\n" );
 #  endif
 
+   if ( SrcTerms.Lightbulb )
+   {
+#     if ( !defined NEUTRINO_SCHEME  ||  NEUTRINO_SCHEME != LIGHTBULB )
+         Aux_Error( ERROR_INFO, "NEUTRINO_SCHEME != LIGHTBULB for the lightbulb scheme !!\n" );
+#     endif
+   }
+
+   if ( SrcTerms.Leakage )
+   {
+#     if ( !defined NEUTRINO_SCHEME  ||  NEUTRINO_SCHEME != LEAKAGE )
+         Aux_Error( ERROR_INFO, "NEUTRINO_SCHEME != LEAKAGE for the leakage scheme !!\n" );
+#     endif
+
+      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, "For %s != 1, the value (%d) must be 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, "SRC_LEAKAGE_RADIUSMIN_LOG (%d) / SRC_LEAKAGE_BINSIZE_RADIUS (%d) > SRC_LEAKAGE_NRADIUS (%d) !!\n",
+                    SrcTerms.Leakage_RadiusMin_Log, SrcTerms.Leakage_BinSize_Radius, SrcTerms.Leakage_NRadius );
+   }
+
 // warning
 // ------------------------------
    if ( MPI_Rank == 0 ) {