From 9d53b17cd020dbc2cf750632d07b0f90a8887bf8 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Fri, 6 Oct 2017 13:16:04 +0800 Subject: clusters/solver.py: Update comments/references/descriptions etc. --- fg21sim/extragalactic/clusters/solver.py | 56 ++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 25 deletions(-) diff --git a/fg21sim/extragalactic/clusters/solver.py b/fg21sim/extragalactic/clusters/solver.py index fb7c52b..a26a9e8 100644 --- a/fg21sim/extragalactic/clusters/solver.py +++ b/fg21sim/extragalactic/clusters/solver.py @@ -78,7 +78,10 @@ class FokkerPlanckSolver: Q(x,t) : injection coefficient (>=0) T(x,t) : escape coefficient - NOTE: The no-flux boundary condition is used. + NOTE + ---- + The no-flux boundary condition is used, and optional boundary fix + may be applied. Parameters ---------- @@ -94,7 +97,7 @@ class FokkerPlanckSolver: Function f(x,t) to calculate the diffusion coefficient C(x,t) f_injection : function Function f(x,t) to calculate the injection coefficient Q(x,t) - f_escape : optional + f_escape : function, optional Function f(x,t) to calculate the escape coefficient T(x,t) buffer_np : int, optional Number of grid points taking as the buffer region near the lower @@ -118,10 +121,12 @@ class FokkerPlanckSolver: References ---------- - [1] Park & Petrosian 1996, ApJS, 103, 255 - http://adsabs.harvard.edu/abs/1996ApJS..103..255P - [2] Donnert & Brunetti 2014, MNRAS, 443, 3564 - http://adsabs.harvard.edu/abs/2014MNRAS.443.3564D + .. [park1996] + Park & Petrosian 1996, ApJS, 103, 255 + http://adsabs.harvard.edu/abs/1996ApJS..103..255P + .. [donnert2014] + Donnert & Brunetti 2014, MNRAS, 443, 3564 + http://adsabs.harvard.edu/abs/2014MNRAS.443.3564D """ def __init__(self, xmin, xmax, x_np, tstep, @@ -153,15 +158,16 @@ class FokkerPlanckSolver: dx[i] = (x[i+1] - x[i-1]) / 2 - NOTE: - Extrapolate the x grid by 1 point beyond each side, therefore + NOTE + ---- + Extrapolate the X grid by 1 point beyond each side, therefore avoid NaN for the first and last element of dx[i]. - Otherwise, the following calculation of tridiagonal coefficients - may be invalid on the boundary elements. + Otherwise, the subsequent calculation of tridiagonal coefficients + may be invalid for the boundary elements. - References: Ref.[1],Eq.(8) + References: Ref.[park1996],Eq.(8) """ - x = self.x + x = self.x # log scale # Extrapolate the x grid by 1 point beyond each side x2 = np.concatenate([ [x[0]**2/x[1]], @@ -179,7 +185,7 @@ class FokkerPlanckSolver: dx[i+1/2] = x[i+1] - x[i] Thus the last element is NaN. - References: Ref.[1],Eq.(8) + References: Ref.[park1996],Eq.(8) """ x = self.x dx_ = x[1:] - x[:-1] @@ -207,7 +213,7 @@ class FokkerPlanckSolver: X[i+1/2] = (X[i] + X[i+1]) / 2 Thus the last element is NaN. - References: Ref.[1],Eq.(10) + References: Ref.[park1996],Eq.(10) """ Xmid = (X[1:] + X[:-1]) / 2 return np.concatenate([Xmid, [np.nan]]) @@ -225,7 +231,7 @@ class FokkerPlanckSolver: @staticmethod def W(w): - # References: Ref.[1],Eqs.(27,35) + # References: Ref.[park1996],Eqs.(27,35) with np.errstate(invalid="ignore"): # Ignore NaN's w = np.abs(w) @@ -239,10 +245,8 @@ class FokkerPlanckSolver: @staticmethod def bound_w(w, wmin=1e-8, wmax=1e3): """ - Bound the absolute values of w within [wmin, wmax]. - - To avoid the underflow/overflow during later W/Wplus/Wminus - calculations. + Bound the absolute values of ``w`` within [wmin, wmax], to avoid + the underflow/overflow during later W/Wplus/Wminus calculations. """ with np.errstate(invalid="ignore"): # Ignore NaN's @@ -254,14 +258,14 @@ class FokkerPlanckSolver: return ww def Wplus(self, w): - # References: Ref.[1],Eq.(32) + # References: Ref.[park1996],Eq.(32) ww = self.bound_w(w) W = self.W(ww) Wplus = W * np.exp(ww/2) return Wplus def Wminus(self, w): - # References: Ref.[1],Eq.(32) + # References: Ref.[park1996],Eq.(32) ww = self.bound_w(w) W = self.W(ww) Wminus = W * np.exp(-ww/2) @@ -281,7 +285,7 @@ class FokkerPlanckSolver: invalid. Therefore, b[0] and b[N-1] should be alternatively calculated with (e.g., no-flux) boundary condition considered. - References: Ref.[1],Eqs.(16,18,34) + References: Ref.[park1996],Eqs.(16,18,34) """ x = self.x dx = self.dx @@ -330,16 +334,18 @@ class FokkerPlanckSolver: is determined by fitting to the data points of ``self.buffer_np`` cells on the upper side of the buffer region. - References: Ref.[2],Sec.(3.3) + TODO: also fix the upper boundary in the same way? + + References: Ref.[donnert2014],Sec.(3.3) """ if self.buffer_np is None: return uc if (uc <= 0.0).sum() > 0: - logger.warning("solved density has negative values!") + logger.warning("solved density has zero/negative values!") return uc x = self.x - # Calculate the power-law index from ``self.buffer_np`` data + # Calculate the power-law index for ``self.buffer_np`` from data # points just right of the buffer region. xp = x[self.buffer_np:(self.buffer_np*2)] yp = uc[self.buffer_np:(self.buffer_np*2)] -- cgit v1.2.2