From 480ae1b0caa2b5d8260df21dc7f541e6c79fdec6 Mon Sep 17 00:00:00 2001 From: Aaron LI Date: Sun, 8 Jan 2017 14:05:28 +0800 Subject: solver.py: Avoid possible overflow when w is too large --- fg21sim/extragalactic/clusters/solver.py | 10 ++++++++++ 1 file changed, 10 insertions(+) (limited to 'fg21sim/extragalactic/clusters/solver.py') diff --git a/fg21sim/extragalactic/clusters/solver.py b/fg21sim/extragalactic/clusters/solver.py index 9d4a781..21b9a68 100644 --- a/fg21sim/extragalactic/clusters/solver.py +++ b/fg21sim/extragalactic/clusters/solver.py @@ -216,6 +216,16 @@ class FokkerPlanckSolver: C_mhalf = self.X_mhalf(C) w_phalf = dx_phalf * B_phalf / C_phalf w_mhalf = dx_mhalf * B_mhalf / C_mhalf + # Avoid overflow when w is too large + w_max = 300 + with np.errstate(invalid="ignore"): + mask_phalf = (np.abs(w_phalf) > w_max) + mask_mhalf = (np.abs(w_mhalf) > w_max) + w_phalf[mask_phalf] = w_max * (np.sign(w_phalf[mask_phalf])) + w_mhalf[mask_mhalf] = w_max * (np.sign(w_mhalf[mask_mhalf])) + logger.warning("Fixed too large w values: " + + "%d+%d" % (len(mask_phalf), len(mask_mhalf))) + # W_phalf = self.W(w_phalf) W_mhalf = self.W(w_mhalf) Wplus_phalf = W_phalf * np.exp(w_phalf/2) -- cgit v1.2.2