numpy function giving incorrect results - checked by hand and excel

Question

I'm writing some functions in numpy for rock physics modelling and have noticed that one of my functions gives erroneous results. The function is my implimentation of Hertz-Mindlin sphere modelling:

Summary of the Hertz-Mindlin model

Here is my function currently:

# Hertz-Mindlin sphere pack model: 

import numpy as np 

def hertzmindlin(K0, G0, PHIC, P, f=1.0):
'''
Hertz-Mindlin sphere-pack model, adapted from:
'Dvorkin, J. and Nur, A., 1996. Elasticity of high-porosity sandstones: 
Theory for two North Sea data sets. Geophysics, 61(5), pp.1363-1370."

Arguments:
K0 = Bulk modulus of mineral in GPa
G0 = Shear modulus of mineral in GPa
PHIC = Critical porosity for mineral-fluid mixture. Calculate using Dvorkin-Nuir (1995) or use literature
P = Confining pressure in GPa
f = Shear modulus correction factor. Default = 1

Results:
V0 = Theoretical poissons ratio of mineral
n = Coordination number of sphere-pack, calculated from Murphy's (1982) empirical relation
K_HM = Hertz-Mindlin effective dry Bulk modulus at pressure, P, in GPa
G_HM = Hertz-Mindlin effective dry Shear modulus at pressure, P, in GPa

'''
V0 = (3*K0-2*G0)/(6*K0+2*G0) # Calculated theoretical poissons ratio of bulk rock
n = 20-(34*PHIC)+(14*(PHIC**2)) # Coordination number at critical porosity (Murphy 1982)
K_HM = (P*(n**2*(1-PHIC)**2*G0**2) / (18*np.pi**2*(1-V0)**2))**(1/3)
G_HM = ((2+3*f-V0*(1+3*f))/(5*(2-V0))) * ((P*(3*n**2*(1-PHIC)**2*G0**2)/(2*np.pi**2*(1-V0)**2)))**(1/3)
return K_HM, G_HM

The problem is that when I run this function for inputs of:

K, G, = 36, 45

PHIC = 0.4

P = 0.001

I get a result of K_HM = 1.0, G_HM = 0.49009009009009

The hand calculated and excel calculated values show this is incorrect, I should be outputting K_HM = 0.763265313, G_HM = 1.081083984

I am fairly certain something is going wrong in the function based on the fact that for the inputs K, G, the output G should be larger than K (it is currently smaller)

Any help would be appreciated! I can do this in excel, but ideally want everything running in python.

unutbu · Accepted Answer · 2017-12-05T19:06:52.577

In Python2, division of integers (using /) returns an integer. For example, 1/3 = 0. In Python3, division of integers (using /) may return a float.

It appears you are using Python2. To get floating-point division (in both Python2 and Python3), ensure each division operation involves at least one float: for example, change 1/3 to 1.0/3 or 1/3.0 or (acceptable but perhaps less readable, 1/3.):

import numpy as np
def hertzmindlin(K0, G0, PHIC, P, f=1.0):
    K0, G0 = map(float, (K0, G0))
    V0 = (3*K0-2*G0)/(6*K0+2*G0) # Calculated theoretical poissons ratio of bulk rock
    n = 20-(34*PHIC)+(14*(PHIC**2)) # Coordination number at critical porosity (Murphy 1982)
    K_HM = (P*(n**2*(1-PHIC)**2*G0**2) / (18*np.pi**2*(1-V0)**2))**(1/3.0)
    G_HM = ((2+3*f-V0*(1+3*f))/(5*(2-V0))) * ((P*(3*n**2*(1-PHIC)**2*G0**2)/(2*np.pi**2*(1-V0)**2)))**(1/3.0)
    return K_HM, G_HM

K, G, = 36, 45
PHIC = 0.4
P = 0.001

print(hertzmindlin(K, G, PHIC, P))

Alternatively, in later versions of Python2 (such as Python2.7) you could place

from __future__ import division

at the top of your script (before all other import statements) to activate Python3-style floating-point division.

Awesome! Figured it would be something relating to number format. I'm unfortunately unable to change my python version or environment at my institution, but will consider this in the future! — 8556732, Dec 05 '17 at 19:05

numpy function giving incorrect results - checked by hand and excel

1 Answers1