I have a code that I am trying to copy from IDL into Python. I have most of it copied and it works, except that when I try to run the function I wrote (called broyfunc) through the broyden1 non-linear solver, I get the error message that I can't have an array with infs or NaNs. I don't know how to go about checking this either. I've included my error messages and my code. Thank you!
import numpy as np
import scipy
import sympy
import sympy.mpmath as mpmath
# Inputs needed for AN_RC_MOD
p0 = 1.5 # reference pressure (surface pressure) [bar]
T0 = 94 # temperature at reference pressure [K]
Ttoa = 175 # temperature at top of atmosphere (assumed as stratopause) [K]
F1 = 1.5 # flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2]
F2 = 1.1 # flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2]
F20 = 0.75 # flux absorbed in channel '2' down to reference level p0 [W/m**2]
Fi = 0 # internal heat flux [W/m**2]
gam = 1.4 # ratio of specific heats, cp/cv (see Eq. 8 of RC12)
alpha = 0.77 # empirical adjustment to dry adiabatic lapse rate (see Eq. 10 of RC12)
n = 4/3 # scaling parameter that relates tau ~ p^n, where tau is gray thermal optical depth (see Eq.6 of RC12)
beta = alpha*((gam-1)/gam)
# Outputs:
# taurc - gray thermal optical depth of radiative-convective boundary
# tau0 - gray thermal optical depth down to reference pressure p0
# k1 - attenuation parameter for channel '1' [dimensionless]
# k2 - attenuation parameter for channel '2' [dimensionless]
# Important constants
sigma = 5.67e-8 # Stefan-Boltzmann constant [J/s/m**2/K**4]
D = 1.66 # diffusivity factory, (see, e.g., Rodgers & Walshaw 1966)
# Need inputs for: (numbers are for Titan, from Example_w_fluxes.pro)
#F1C = F1 ;flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2]
#F20C = F20 ;flux absorbed in channel '2' down to reference level p0 [W/m**2]
#F2C = F2 ;flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2]
#FiC = Fi ;internal heat flux [W/m**2]
#T0C = T0 ;temperature at reference pressure [K]
#p0C = p0 ;reference pressure (e.g., surface pressure) [bar]
#TtoaC = Ttoa ;temperature at top of atmosphere [K]
#nC = n ;scaling power between gray optical depth and pressure
#betaC = alpha*(gam - 1.)/gam ;beta defined in Eq. 11 of RC12
# Need inputs for: (numbers are for Titan, from Example_w_fluxes.pro)
#F1C = F1 ;flux absorbed in channel '1' (assumed to be "stratosphere") [W/m**2]
#F20C = F20 ;flux absorbed in channel '2' down to reference level p0 [W/m**2]
#F2C = F2 ;flux absorbed in channel '2' (assumed to be "troposphere") [W/m**2]
#FiC = Fi ;internal heat flux [W/m**2]
#T0C = T0 ;temperature at reference pressure [K]
#p0C = p0 ;reference pressure (e.g., surface pressure) [bar]
#TtoaC = Ttoa ;temperature at top of atmosphere [K]
#nC = n ;scaling power between gray optical depth and pressure
#betaC = alpha*(gam - 1.)/gam ;beta defined in Eq. 11 of RC12
# x[0] = optical depth of Radiative-Convective boundary, taurc
# x[1] = total thermal optical depth down to reference pressure p0 (tau0)
# initial guess that R-C boundary is at D*taurc~1, and tau0~10*p0^n
x = [1/D, 10*p0**n] # initial guesses of taurc and tau0, respectively
print x
# Define function F for Broyden optimization
def broyfunc(x):
# Variables
taurc = x[0] # current value of optical depth of r-c boundary
tau0 = x[1] # current value of total thermal optical depth down to reference pressure p0
# Compute fluxes and temperatures at r-c boundary ;;;
# first check if solver is attempting to use negative values for the
# optical depths (which is unphysical), and, if so, return large numbers
# to the solver to prevent it from trying such values as solutions
if x[0] < 0 or x[1] < 0:
x = [1.e3,1.e3]
print x
else:
x = x
# Shortwave optical depth in channel "2" is given by
# LOG( F2/(F2-F20) ), and k2 is defined as the ratio of
# this to the thermal optical depth at p0 (i.e., tau0)
if F2 > 0:
k2 = math.log(F2/(F2-F20))/tau0
else:
# if the flux in channel "2" is zero,
# then no attenuation in this channel
k2 = 0
# k1 is computed by requiring the radiative temperature profile
# to equal the top-of-atmosphere temperature (Ttoa, set by user),
# which is accomplished by evaluating Eq. 18 in RC12 at tau=0
if F1 > 0:
k1 = D*( 2/F1*( sigma*Ttoa**4 - Fi/2 - F2/2*(1 + k2/D) ) - 1 )
else:
# if the flux in channel "1" is zero,
# then no attenuation in this channel
k1 = 0
# return values of k1, k2 to user
k1C = k1 # attenuation parameter for channel '1' (assumed as "stratosphere")
k2C = k2 # attenuation parameter for channel '2' (assumed as "troposphere")
print 'k1C', k1C
print 'k2C', k2C
# four cases: (1) both k1, k2 are greater than zero
# (2) k1 greater than zero, k2 equal to (or less than) zero
# (3) k2 greater than zero, k1 equal to (or less than) zero
# (4) both k1, k2, are equal to (or less than) zero
# case 1: both k1, k2 are greater than zero #
if k1 > 0 and k2 > 0:
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = (F1/2)*(1 + (D/k1) + (1 - D/k1)*math.exp(-k1*taurc)) + (F2/2)*(1 + (D/k2) + (1 - D/k2)*math.exp(-k2*taurc)) + (Fi/2)*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D/k1 + (k1/D - D/k1)*math.exp(-k1*taurc)) + F2/2*(1 + D/k2 + (k2/D - D/k2)*math.exp(-k2*taurc)) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# case 2: k1 greater than zero, k2 equal to (or less than) zero #
if k1 > 0 and k2 < 0:
k2C = 0
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = F1/2*(1 + D/k1 + (1 - D/k1)*math.exp(-k1*taurc)) + F2/2*(2 + D*taurc) + Fi/2*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D/k1 + (k1/D - D/k1)*math.exp(-k1*taurc)) + F2/2*(1 + D*taurc) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# case 3: k2 greater than zero, k1 equal to (or less than) zero #
if k1 < 0 and k2 > 0:
k1C = 0
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = F1/2*(2 + D*taurc) + F2/2*(1 + D/k2 + (1 - D/k2)*math.exp(-k2*taurc)) + Fi/2*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D*taurc) + F2/2*(1 + D/k2 + (k2/D - D/k2)*math.exp(-k2*taurc)) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# case 4: both k1, k2, are equal to (or less than) zero
if k1 < 0 and k2 < 0:
k1C = 0
k2C = 0
# upwelling thermal flux at R-C boundary from radiative equilibrium (RC12 Eq. 19)
Fup_r = F1/2*(2 + D*taurc) + F2/2*(2 + D*taurc) + Fi/2*(2 + D*taurc)
# temperature at R-C boundary from radiative equilibrium (RC12 Eq. 18)
sigmaT_r = F1/2*(1 + D*taurc) + F2/2*(1 + D*taurc) + Fi/2*(1 + D*taurc)
T_r = (sigmaT_r/sigma)**0.25
# Upwelling thermal flux at R-C boundary from convective region (RC12 Eq. 13)
#Fup_c = sigma*T0**4*math.exp(D*taurc)*((math.exp(-D*tau0)) + (1/((D*tau0)**(4*beta/n)))*(mpmath.gammainc(1+(4*beta/n),D*tau0) - mpmath.gammainc(1+(4*beta/n),D*taurc)))
#Fup_c = (sigma*T0**4)*math.exp(D*taurc)*(math.exp((-D*tau0) + 1/(D*tau0)**(4*beta/n)*mpmath.gammainc(1+4*beta/n)*(mpmath.gammainc(1+4*beta/n, D*tau0) - mpmath.gammainc(1+4*beta/n,D*taurc))))
a1 = D*taurc
a2 = -D*tau0
a3 = D*taurc
a4 = D*tau0
b1 = (4*beta)/n
b2 = b1 + 1
Fup_c1 = (sigma*T0**4)*(math.exp(a1))
Fup_c2 = math.exp(a2)
Fup_c3 = 1/(a4**b1)
Fup_c4 = mpmath.gammainc(b2, a3)
Fup_c5 = mpmath.gammainc(b2, a4)
Fup_c = Fup_c1*(Fup_c2 + (Fup_c3*(Fup_c4 - Fup_c5)))
# temperature at R-C boundary from convective region (RC12 Eq. 11)
T_c = T0*(taurc/tau0)**(beta/n)
print 'D*taurc', D*taurc
print 'Fup_r', Fup_r
print 'Fup_c', Fup_c
print 'T_r', T_r
print 'T_c', T_c
# system is solved when the temperatures and upwelling fluxes are continuous
# at the R-C boundary, so the solver attempts to minimize the difference
# between these quantities as computed from the radiative expressions and the
# convective expressions
#result = np.dtype('f8')
result = [T_r - T_c, Fup_r - Fup_c]
print result
# essentially we solve simultaneous equations:
# T_r - T_c = 0 and Fup_r - Fup_c = 0
broyfunc(x)
#y = np.dtype('f8')
y = scipy.optimize.broyden1(broyfunc, x)
ValueError Traceback (most recent call last) in () 183 y = np.dtype('f8') 184 print x --> 185 y = scipy.optimize.broyden1(broyfunc, x, iter=10) 186
//anaconda/lib/python2.7/site-packages/scipy/optimize/nonlin.pyc in broyden1(F, xin, iter, alpha, reduction_method, max_rank, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, **kw)
//anaconda/lib/python2.7/site-packages/scipy/optimize/nonlin.pyc in nonlin_solve(F, x0, jacobian, iter, verbose, maxiter, f_tol, f_rtol, x_tol, x_rtol, tol_norm, line_search, callback, full_output, raise_exception) 279 dx = np.inf 280 Fx = func(x) --> 281 Fx_norm = norm(Fx) 282 283 jacobian = asjacobian(jacobian)
//anaconda/lib/python2.7/site-packages/scipy/linalg/misc.pyc in norm(a, ord) 115 # Differs from numpy only in non-finite handling and the use of 116 # blas --> 117 a = np.asarray_chkfinite(a) 118 if ord in (None, 2) and (a.ndim == 1) and (a.dtype.char in 'fdFD'): 119 # use blas for fast and stable euclidean norm
//anaconda/lib/python2.7/site-packages/numpy/lib/function_base.pyc in asarray_chkfinite(a, dtype, order) 588 if a.dtype.char in typecodes['AllFloat'] and not np.isfinite(a).all(): 589 raise ValueError( --> 590 "array must not contain infs or NaNs") 591 return a 592
ValueError: array must not contain infs or NaNs
The most immediate problem is that you're not returning anything from broyfunc:
#result = np.dtype('f8')
result = [T_r - T_c, Fup_r - Fup_c]
print result
and so all the broyden integrator is seeing is effectively
>>> scipy.optimize.broyden1(lambda x: None, 0)
Traceback (most recent call last):
[...]
ValueError: array must not contain infs or NaNs
(Terrible error message, admittedly.)
After adding return result I get something like
[...]
k1C 114.773687493
k2C 0.0688650670293
D*taurc 23.3390821136
Fup_r 9.7427665389
Fup_c 3.96039873101838
T_r 113.899157512
T_c 90.592872808
[23.306284704464417, mpf('5.782367807885139')]
[1000.0, 1000.0]
k1C 114.8261008
k2C -0.00260762458087
Traceback (most recent call last):
[...]
OverflowError: math range error
but that's a separate issue and would require actually thinking about what your code is supposed to be doing. ;-)
There's another issue it's worth warning about: you write
n = 4/3 # scaling parameter that relates tau ~ p^n, where tau is gray thermal optical depth (see Eq.6 of RC12)
but you're using Python 2, and so you have integer division:
>>> n = 4/3
>>> n
1
Either upgrade to Python 3 (preferable), ensure that everything's a float by adding a decimal point (4./3), or -- and this is what I'd recommend if you're stuck with Python 2 and are working on numerical code -- add from __future__ import division to the start of your program so that you get
>>> from __future__ import division
>>> 4/3
1.3333333333333333