KeyError in a dictionary when in the previous iteration of a while loop the key:value pair should have been created

Question

I am using a while loop to integrate various quantities from the surface of a star inwards using appropriate boundary conditions and stellar structure equations.

I am using dictionaries to represent physical variables such as pressure and density, where the plan is for the radii to be keys, and the value to be the pressure or density.

I have a key:value pair for the surface, and then I step inwards iteratively using a while loop updating the dictionaries as below:

import constants
import math
import matplotlib.pyplot as plt


mass=3*constants.solar_mass
radius=1.5*constants.solar_radius
#Variables to guess
core_temperature=1.4109*10**7
core_pressure= 2.6851*10**14
luminosity=pow(3,3.5)*constants.solar_luminosity

            

#Functions we are searching for
temperature={}
#guess
temperature[0]=core_temperature
#From the steffan boltzmann law
temperature[radius]=pow(luminosity/(4*math.pi*pow(radius,2)*constants.stefan_boltzmann_constant),0.25)

        
pressure={}
#guess
pressure[0]=core_pressure
#Pressure surface boundary condition
pressure[radius]=(2*constants.gravitation_constant*mass)/(3*constants.opacity*pow(radius,2))

            
mass_enclosed={}
#boundary conditions
mass_enclosed[0]=0
mass_enclosed[radius]=mass


density={}
#density surface boundary condition
density[radius]=(constants.mean_molecular_weight*pressure[radius])/(constants.gas_constant*temperature[radius])


delta_radius=int(radius/100)

#Polytropic constant
K=(pressure[radius]*constants.mean_molecular_weight)/(constants.gas_constant*pow(density[radius],constants.adiabatic_constant))


def integrate_from_surface():
    i=0
    while radius-i*delta_radius>(0.5*radius):
        #temporary radius just for each loop through
        r=radius-i*delta_radius
        

        #updating pressure
        pressure[r-delta_radius]=pressure[r]+(density[r]*constants.gravitation_constant*mass_enclosed[r]*delta_radius)/pow(r,2)

        #updating density
        density[r-delta_radius]=pow((pressure[r-delta_radius]*constants.mean_molecular_weight)/(constants.gas_constant*K),1.0/constants.adiabatic_constant)

        #updating mass enclosed
        mass_enclosed[r-delta_radius]=mass_enclosed[r]-4*math.pi*pow(r,2)*delta_radius*density[r]

        i=i+1
    
    

integrate_from_surface()

While Loop: Radius and dictionaries are defined above

I am getting a KeyError, as shown below:

Traceback (most recent call last):
  File "main.py", line 63, in <module>
    integrate_from_surface()
  File "main.py", line 51, in integrate_from_surface
    pressure[r-delta_radius]=pressure[r]+(density[r]*constants.gravitation_constant*mass_enclosed[r]*delta_radius)/pow(r,2)
KeyError: 1043966868.09

KeyError message

If I print out the variable r in the while group, the process works perfectly until r is 1043966868.09. I do not understand, surely on the previous iteration I made this a key, so there should be no KeyError.

Constants file below:

solar_mass=1.9891*10**30

solar_radius=6.9598*10**8

solar_luminosity=3.8515*10**26

gas_constant=8.3145*10**3

gravitation_constant=6.6726*10**-11

radiation_constant=7.5646*10**-16

speed_of_light=2.9979*10**8

stefan_boltzmann_constant= radiation_constant*speed_of_light * 0.25



opacity=0.034

adiabatic_constant = 5.0/3

mean_molecular_weight = 8.0/13

Thanks in advance for any help.

Answer 1

As I stated in the comments, this behaviour is probably linked to the difficulty of hashing floats. In the end, using floats as dictionary keys is only as precise as your float precision. There are already extensive articles about the consequences and mechanisms at work here, for instance this one.

As an example (credits to the aforementioned article), you can see that hash(10-9.8) == hash(.2) returns False : using such keys in a dictionnary would create two entries.

A workaround (if you want to stick to float keys) would be to first evaluate all possible keys, then reuse those as many time as needed.

There is another good reason for this approch, as you will have to rewrite your "while loop" and replace it with a "for loop" : while loops are slower than for loops (more hints in here).

In your case, you could evaluate your steps this way :

radiuses = [radius-i*delta_radius for i in range(1,50)]

Afterwards, switching the loop iteration is very straightforward (for r in radiuses:). And more to the point, your code won't raise any exception.

Note that you could also store all your datas in a pandas.DataFrame, combining the way of accessing you steps by order (using .iloc) or by radius (using .loc).

That could be used this way (note that I'm storing all intermediate results in lists rather than dictionaries, this is in relation with the pandas.DataFrame constructor ) :

import math
import pandas as pd

SOLAR_MASS = 1.9891*10**30
SOLAR_RADIUS = 6.9598*10**8
SOLAR_LUMINOSITY = 3.8515*10**26
RADIATION_CONSTANT = 7.5646*10**-16
SPEED_OF_LIGHT = 2.9979*10**8
STEFAN_BOLTZMANN = RADIATION_CONSTANT * SPEED_OF_LIGHT * 0.25
OPACITY = 0.034
GRAVITATION_CONSTANT = 6.6726*10**-11
MEAN_MOLECULAR_WEIGHT = 8.0/13
GAS_CONSTANT = 8.3145*10**3
ADIABATIC_CONSTANT = 5.0/3

mass=3*SOLAR_MASS
radius=1.5*SOLAR_RADIUS
luminosity=(3**3.5)*SOLAR_LUMINOSITY

TEMPERATURE_R = (luminosity/(4*math.pi*(radius**2)*STEFAN_BOLTZMANN))**0.25
PRESSURE_R = (2*GRAVITATION_CONSTANT*mass)/(3*OPACITY*(radius**2))
MASS_R = mass
DENSITY_R = (MEAN_MOLECULAR_WEIGHT*PRESSURE_R)/(GAS_CONSTANT*TEMPERATURE_R)
K=(PRESSURE_R*MEAN_MOLECULAR_WEIGHT)/(GAS_CONSTANT*DENSITY_R**ADIABATIC_CONSTANT)


def integrate_from_surface_df():
    delta_radius=int(radius/100)
    pressure=[PRESSURE_R]
    mass_enclosed=[MASS_R]
    density=[DENSITY_R]
    radiuses = [radius-i*delta_radius for i in range(1,50)]
    for r in radiuses:
        pressure.append(pressure[-1]+(density[-1]*GRAVITATION_CONSTANT*mass_enclosed[-1]*delta_radius)/r**2)
        density.append(((pressure[-2]*MEAN_MOLECULAR_WEIGHT)/(GAS_CONSTANT*K))**(1.0/ADIABATIC_CONSTANT))
        mass_enclosed.append(mass_enclosed[-1]-4*math.pi*r**2*delta_radius*density[-2])

    df = pd.DataFrame({"pressure":pressure, "density":density, "mass_enclosed":mass_enclosed}, index=[radius]+radiuses)
    return df

print(integrate_from_surface_df())

Which would return :

              pressure      density    mass_enclosed
1.043970e+09  7.163524e+03  0.000043   5.967300e+30
1.033530e+09  1.743478e+05  0.000043   5.967300e+30
1.023091e+09  3.449615e+05  0.000292   5.967300e+30
1.012651e+09  1.527176e+06  0.000439   5.967300e+30
1.002211e+09  3.344824e+06  0.001072   5.967300e+30
9.917715e+08  7.876679e+06  0.001716   5.967300e+30
9.813318e+08  1.528565e+07  0.002870   5.967300e+30
9.708921e+08  2.793971e+07  0.004271   5.967299e+30
9.604524e+08  4.718766e+07  0.006134   5.967299e+30
9.500127e+08  7.543907e+07  0.008400   5.967298e+30
9.395730e+08  1.149941e+08  0.011132   5.967297e+30
9.291333e+08  1.685944e+08  0.014335   5.967296e+30
9.186936e+08  2.391984e+08  0.018035   5.967294e+30
9.082539e+08  3.300760e+08  0.022246   5.967292e+30
8.978142e+08  4.447982e+08  0.026988   5.967290e+30
8.873745e+08  5.872670e+08  0.032278   5.967287e+30
8.769348e+08  7.617399e+08  0.038133   5.967284e+30
8.664951e+08  9.728621e+08  0.044575   5.967280e+30
8.560554e+08  1.225701e+09  0.051622   5.967276e+30
8.456157e+08  1.525790e+09  0.059297   5.967271e+30
8.351760e+08  1.879167e+09  0.067624   5.967265e+30
8.247363e+08  2.292433e+09  0.076627   5.967259e+30
8.142966e+08  2.772804e+09  0.086334   5.967253e+30
8.038569e+08  3.328175e+09  0.096773   5.967245e+30
7.934172e+08  3.967188e+09  0.107976   5.967237e+30
7.829775e+08  4.699315e+09  0.119976   5.967229e+30
7.725378e+08  5.534939e+09  0.132808   5.967219e+30
7.620981e+08  6.485454e+09  0.146512   5.967209e+30
7.516584e+08  7.563373e+09  0.161127   5.967198e+30
7.412187e+08  8.782447e+09  0.176699   5.967187e+30
7.307790e+08  1.015780e+10  0.193274   5.967174e+30
7.203393e+08  1.170609e+10  0.210904   5.967161e+30
7.098996e+08  1.344566e+10  0.229642   5.967147e+30
6.994599e+08  1.539674e+10  0.249548   5.967133e+30
6.890202e+08  1.758168e+10  0.270684   5.967117e+30
6.785805e+08  2.002515e+10  0.293117   5.967101e+30
6.681408e+08  2.275445e+10  0.316921   5.967083e+30
6.577011e+08  2.579981e+10  0.342172   5.967066e+30
6.472614e+08  2.919473e+10  0.368956   5.967047e+30
6.368217e+08  3.297638e+10  0.397363   5.967027e+30
6.263820e+08  3.718607e+10  0.427491   5.967007e+30
6.159423e+08  4.186974e+10  0.459445   5.966985e+30
6.055026e+08  4.707856e+10  0.493339   5.966963e+30
5.950629e+08  5.286959e+10  0.529296   5.966940e+30
5.846232e+08  5.930656e+10  0.567451   5.966917e+30
5.741835e+08  6.646075e+10  0.607948   5.966892e+30
5.637438e+08  7.441197e+10  0.650945   5.966867e+30
5.533041e+08  8.324981e+10  0.696611   5.966841e+30
5.428644e+08  9.307487e+10  0.745135   5.966814e+30
5.324247e+08  1.040004e+11  0.796717   5.966786e+30

There is also an other workaround, which you already guessed by yourself : that would be to use integers as keys in your dictionary. I don't think this to be awkward in any way, but I'd rather have keys with a real, direct meaning (that is, your radius) rather than the iteration step... But this is really up to you.