1
\$\begingroup\$

Please see the following code. I am using it to calculate the double integration. Please help me to improve the code.

import numpy as np
from scipy import integrate
from scipy.special import j1
from scipy.special import j0
from numpy import sin
from numpy import cos
from numpy import pi
import time
import numba as nb
from numba import jit
from numba import cfunc
from numba.types import intc, CPointer, float64
from scipy import LowLevelCallable
from matplotlib import pyplot as plt

q = np.linspace(0.03, 1.0, 1000)
start = time.time()
length = 20000
pd = 0.1
radius = 40
start = time.time()

def jit_integrand_function(integrand_function):
    jitted_function = jit(integrand_function, nopython=True)
    @cfunc(float64(intc, CPointer(float64)),error_model="numpy",fastmath=True)
    def wrapped(n, xx):
        ar = nb.carray(xx, n)
        return jitted_function(ar[0], ar[1], ar[2])
    return LowLevelCallable(wrapped.ctypes)

@jit_integrand_function
def f(theta,z,q):
    def ff(theta, z, q):
        l = length / 2
        qz = q * sin(theta)
        qr = q * cos(theta)
        return (4 * pi * sin(qz * l) * z * j1(qr * z) / (qr * qz)) ** 2

    def g_distribution(z, radius, pd):
        c = radius
        s = pd * c
        return (1 / s)) * np.exp(-0.5 * (((z - radius) / s) ** 2))

    return ff(theta,z,q)*g_distribution(z,radius,pd)* sin(theta)



def lower_inner(z):
    return 0.

def upper_inner(z):
    return 120.

y = np.empty(len(q))
for n in range(len(q)):
    y[n] = integrate.dblquad(f, 0, 1*pi/180, lower_inner, upper_inner,args=(q[n],))[0]

end = time.time()
print(end - start)
plt.loglog(q,y)
plt.show()

Also, I was wondering if I can use GPU for the above code? How much gain I would have for the trouble(I have zero knowledge of programming for GPU use)?

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.