aboutsummaryrefslogtreecommitdiffstats
path: root/fg21sim/extragalactic/clusters/solver.py
diff options
context:
space:
mode:
authorAaron LI <aly@aaronly.me>2017-07-25 20:38:32 +0800
committerAaron LI <aly@aaronly.me>2017-07-25 20:38:32 +0800
commite6b727785fab339103a103194317ffdc6ba37684 (patch)
tree2c71d378ba32af5fe976e18090d9decfbccc8763 /fg21sim/extragalactic/clusters/solver.py
parente58c402015658b19576070fd3d10566d5f1e039c (diff)
downloadfg21sim-e6b727785fab339103a103194317ffdc6ba37684.tar.bz2
clusters/solver.py: Use multiple data points to fit the power-law
Signed-off-by: Aaron LI <aly@aaronly.me>
Diffstat (limited to 'fg21sim/extragalactic/clusters/solver.py')
-rw-r--r--fg21sim/extragalactic/clusters/solver.py36
1 files changed, 19 insertions, 17 deletions
diff --git a/fg21sim/extragalactic/clusters/solver.py b/fg21sim/extragalactic/clusters/solver.py
index d7660b0..67bda58 100644
--- a/fg21sim/extragalactic/clusters/solver.py
+++ b/fg21sim/extragalactic/clusters/solver.py
@@ -322,30 +322,32 @@ class FokkerPlanckSolver:
def fix_boundary(self, uc):
"""
- Truncate the lower end (i.e., near the lower boundary) of the
- distribution/spectrum and then extrapolate as a power law, in order
- to avoid the unphysical pile-up of electrons at the lower regime.
-
- TODO:
- Fit a power law to the same number (``buffer_np``) of data points,
- then extrapolate it to fix the lower buffer region.
+ Due to the no-flux boundary condition adopted, particles may
+ unphysically pile up near the lower boundary. Therefore, a
+ buffer region spanning ``self.buffer_np`` cells is chosen, within
+ which the densities are replaced by extrapolating from the upper
+ density distribution as a power law, and the power-law index
+ 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)
"""
if self.buffer_np is None:
return uc
+ if (uc <= 0.0).sum() > 0:
+ logger.warning("solved density has negative values!")
+ return uc
- uc = np.asarray(uc)
x = self.x
- # Calculate the power-law index
- xa = x[self.buffer_np]
- xb = x[self.buffer_np+1]
- ya = uc[self.buffer_np]
- yb = uc[self.buffer_np+1]
- if ya > 0 and yb > 0:
- # Truncate and extrapolate as a power law
- s = np.log(yb/ya) / np.log(xb/xa)
- uc[:self.buffer_np] = ya * (x[:self.buffer_np] / xa) ** s
+ # Calculate the power-law index from ``self.buffer_np`` 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)]
+ # Power-law fit
+ pfit = np.polyfit(np.log10(xp), np.log10(yp), deg=1)
+ xbuf = x[:self.buffer_np]
+ ybuf = 10 ** np.polyval(pfit, np.log10(xbuf))
+ uc[:self.buffer_np] = ybuf
return uc
def time_step(self):