Skip to content

RuntimeError: Integration called with impossible bounds #2

@vijaymahatma

Description

@vijaymahatma

I am getting this error when running a particular model:

lsoda-- at current t (=r1), mxstep (=i1) steps
taken on this call before reaching tout
in above message, i1 = 500
in above message, r1 = 0.3655406222012D+16
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 3.450883e+58
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 4.955897e+59
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 8.983564e+60
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 1.149010e+61
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 8.583934e+62
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 3.620057e+63
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 8.448133e+63
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 8.515429e+64
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 1.965629e+65
AccuracyWarning)
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.py:709: AccuracyWarning: divmax (10) exceeded. Latest difference = 1.254024e+65
AccuracyWarning)
limit is nan: 0 4.1663449e-316 2.3262218789397e-310

RuntimeError Traceback (most recent call last)
/beegfs/general/mahatmav/analytic/4C39.24/4C39.py in ()
68 for Q in q_list: #min and max jet powers based on B field limits and dynamical ages
69
---> 70 env.solve(Q,tv)
71 env.findb()
72 env.findsynch(408e6) #frequency of spw7 map used to measure radio luminosity

/beegfs/general/mahatmav/analytic/4C39.24/solver.pyc in solve(self, Q, tv, tstop)
465 vt=self.vtot(self.R[i],self.Rp[i])
466 self.vt.append(vt)
--> 467 speeds=self.iter_dLdt([self.R[i],self.Rp[i]],tv[i])
468 vl=self.tempvl # stored here by dLdt
469 self.vl.append(vl)

/beegfs/general/mahatmav/analytic/4C39.24/solver.pyc in iter_dLdt(self, L, t, iterlimit)
343 '''
344 v0=np.array([self.cs,self.cs])
--> 345 r0=self.dL_dt_vest(L,t,vcomb(v0,L))
346 v2=r0
347 r2=self.dL_dt_vest(L,t,vcomb(v2,L))

/beegfs/general/mahatmav/analytic/4C39.24/solver.pyc in dL_dt_vest(self, L, t, v_est)
313 Rp=L[1]
314
--> 315 vl=self.vlobe(R,Rp,v_est,t)
316 self.tempvl=vl
317 if t<=self.tstop:

/beegfs/general/mahatmav/analytic/4C39.24/solver.pyc in vlobe(self, R, Rp, v, t)
271 N=self.ndict[(R,Rp)]
272 except:
--> 273 N=self.intn(R,Rp)
274 self.ndict[(R,Rp)]=N
275 if t>self.tstop:

/beegfs/general/mahatmav/analytic/4C39.24/solver.pyc in intn(self, R, Rp)
255 return 0
256 # factor 2 here since we integrate only one lobe
--> 257 result=2.0*romberg(self._outer_int,0,R,args=(R,Rp))
258 return result
259
/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.pyc in romberg(function, a, b, args, tol, rtol, show, divmax, vec_func)
685 interval = [a, b]
686 intrange = b - a
--> 687 ordsum = _difftrap(vfunc, interval, n)
688 result = intrange * ordsum
689 resmat = [[result]]

/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.pyc in _difftrap(function, interval, numtraps)
558 raise ValueError("numtraps must be > 0 in difftrap().")
559 elif numtraps == 1:
--> 560 return 0.5*(function(interval[0])+function(interval[1]))
561 else:
562 numtosum = numtraps/2

/home/mahatmav/soft/anaconda2/lib/python2.7/site-packages/scipy/integrate/quadrature.pyc in vfunc(x)
117 def vfunc(x):
118 if np.isscalar(x):
--> 119 return func(x, *args)
120 x = np.asarray(x)
121 # call with first point to get output type

/beegfs/general/mahatmav/analytic/4C39.24/solver.pyc in _outer_int(self, z, R, Rp)
247 if np.isnan(limit):
248 print 'limit is nan:',z,R,Rp
--> 249 raise RuntimeError('Integration called with impossible bounds')
250 return romberg(self._inner_int,0,limit,args=(z,R,Rp),divmax=20)
251

RuntimeError: Integration called with impossible bounds

This happens when I run the following model:
kt=2e3
p0=1e-11
rc=10
beta=1.5
q=2.22
z=1.88
zeta=0.01
tmin=0
tmax=200
tv=np.logspace(-6,np.log10(tmax),100)*Myr

env=Evolve_RG('beta',kT=kteV,p0=p0,rc=rckpc,beta=beta,do_adiabatic=True,q=q,z=z,verbose=False,zeta=zeta)
Q=3.1622776601683792e+38
env.solve(Q,tv)

When inspecting the error, I noticed that env.R and env.Rp is not correctly being determined towards the end of the time steps, e.g:
In [26]: env.R
Out[26]:
array([9.45393960e+015, 1.14673115e+016, 1.39094641e+016, 1.68717133e+016,
2.04648221e+016, 2.48231426e+016, 3.01096391e+016, 3.65219820e+016,
4.42999388e+016, 5.37343393e+016, 6.51779507e+016, 7.90586673e+016,
9.58955109e+016, 1.16318037e+017, 1.41089876e+017, 1.71137286e+017,
2.07583786e+017, 2.51792167e+017, 3.05415451e+017, 3.70458695e+017,
4.49353967e+017, 5.45051285e+017, 6.61128921e+017, 8.01610783e+017,
9.63767961e+017, 1.14713836e+018, 1.35386495e+018, 1.58615967e+018,
1.84632736e+018, 2.13679121e+018, 2.46012024e+018, 2.81905825e+018,
3.21655461e+018, 3.65579659e+018, 4.14024412e+018, 4.67366713e+018,
5.26018638e+018, 5.90431854e+018, 6.61102673e+018, 7.38577715e+018,
8.23460360e+018, 9.16418057e+018, 1.01819071e+019, 1.12960028e+019,
1.25156175e+019, 1.38509583e+019, 1.53134343e+019, 1.69158248e+019,
1.86724712e+019, 2.05995010e+019, 2.27150848e+019, 2.50397368e+019,
2.75966624e+019, 3.04121678e+019, 3.35161395e+019, 3.69426151e+019,
4.07304633e+019, 4.49242058e+019, 4.95750169e+019, 5.47419588e+019,
6.04935203e+019, 6.69095673e+019, 7.40838422e+019, 8.21272220e+019,
9.11720216e+019, 1.01377775e+020, 1.12939116e+020, 1.26096710e+020,
1.41152655e+020, 1.58492565e+020, 1.78617810e+020, 2.02193426e+020,
2.30120664e+020, 2.63648850e+020, 3.04550824e+020, 3.55401731e+020,
4.20024613e+020, 5.04196275e+020, 6.16726697e+020, 7.70987065e+020,
9.86825596e+020, 1.29281036e+021, 1.72959617e+021, 2.35683652e+021,
3.26636464e+021, 4.59648234e+021, 6.55179824e+021, 9.44149798e+021,
1.37314563e+022, 2.01223064e+022, 2.96658900e+022, 4.39385697e+022,
6.52985721e+022, 9.72640769e+022, 1.45063975e+023, 2.16609533e+023,
3.24746035e+023, 3.48810809e+023, 1.12465777e-312, 1.06099790e-312])

Any ideas?

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions