Fast way to calculate source response using Green's function

I am a newbie in the scientific computing in python. The function I am trying to calculate is the response of the wave equation for a given source term (For reference, See: Equation 11.67 in https://webhome.phy.duke.edu/~rgb/Class/phy319/phy319/node75.html).

$$\phi(\pmb{x},t) = \int_V G_+(\pmb{x}, t; \pmb{x'}, t)\rho(\pmb{x'}, t')d^3x'dt'$$

To give context, the term phi in the equation 11.67 is the displacement and the term rho can be thought as a source of disturbance.

Now, in my case, the problem is constructed in spatial dimension of 2 (x-y). Thus, I have to iterate the equation for grid points in x, y and t. This makes the overall calculation extremely time-consuming. Below is my lines of code:

import numpy as np

# Assigning grid points

Nx = 100
Ny = 100
Nt = 100

x = np.linspace(-Nx, Nx, 100)
y = np.linspace(-Ny, Ny, 100)
t = np.linspace(0, Nt-1, 100)

# Equation to solve
# phi(x,y,t) = ∫∫∫ G(x,y,t; x',y',t') . Source(x',y',t') dx' dy' dt'
# G(x,y,t; x',y',t') = Green's Function

# phi = displacement by the wave
phi = np.zeros((Nt,Ny,Nx))

# Define Function to realize Green's Function for Wave Equation
def gw(xx, yy, tt, mm, nn, ss):

numer = (tt-ss) - np.sqrt((xx-mm)**2+(yy-nn)**2)
denom = ((tt-ss)**2-(xx-mm)**2-(yy-nn)**2)

if denom < 0:
denom = 0
else:
denom = np.sqrt(denom)

kk = np.heaviside(numer,1)/(2*np.pi*denom+1)
return (kk)

# Define Function to realize Gaussian Disturbance
def g_source(xx, yy, tt):

spatial_sigma  = 5   # spatial width of the Gaussian
temporal_sigma = 30  # temporal width of the Gaussian
onset_time     = 20  # time when Gaussian Disturbance reaches its peak

kk = np.exp((-(np.sqrt((xx)**2+(yy)**2)/spatial_sigma)**2)))*np.exp(-((tt-onset_time)/temporal_sigma)**2)
return (kk)

for k in range (Nt):
for j in range (Ny):
for i in range (Nx):

# this is the Green's function evaluated in each grid points
green = np.array([gw(x[i],y[j],t[k],i_grid,j_grid,k_grid) for k_grid in t for j_grid in y for i_grid in x])
green = green.reshape(Nt, Ny, Nx)

# this is the source term (rho), here it is a Gaussian Disturbance
gauss = np.array([g_source(i_grid,j_grid,k_grid) for k_grid in t for j_grid in y for i_grid in x])
gauss = gauss.reshape(Nt, Ny, Nx)

# performing integration
phi[k,j,i] = sum(sum(sum(np.multiply(green,gauss))))

Please let me know if there is a faster way to do this integration. Thanks in advance. I deeply appreciate your patience and support.

Performance wise

In general:

You use numpy, but you write it almost like in fortran. Python with numpy is good for scientific programing and computing as long as you don't do many loops

• read this first A beginners guide to using Python for performance computing
• if you really need to do tight loops than use cython.
• But most of the time you can avoid that by expressing the loop as some functional operation on the whole array
• most numpy functions can operator on whole array in one pass, use it.
• always prefer y[:]=np.sqrt(x[:]) and avoid
for i in xrange(N):
y[i]=np.sqrt(x[i])
• operation like dot,convolve,einsum,ufunc.reduce are very potent to express most sciencetific algorithms in more abstract way
• you can express most inner if by boolen mask arrays or integer index array

• e.g. f[ f>15.0 ] = 15.0 clamp all elements in array f which are >15.0 to 15.0; you can also store the boolean mask mask = (f>15.0) and than use it like f[mask]+=g[mask]. This way you can express branching as fully functional programs/expressions with arrays.

• Do not construct new np.array (np.zeros, np.ones etc.) too often (i.e. inside tight loops). Optimal performance is obtained if you prepare all arrays you need at the beggining. Avoid frequent conversion between list and array

To address your code example in particular:

• Green's functions are basically convolutions. I'm pretty sure you can express it using e.g. scipy.ndimage.filters.convolve
• if your convolution kernel is large (i.e. pixels interact more than with few neighbors) than it is often much faster to do it in Fourier space (convolution transforms as multiplication) and using np.fftn with O(nlog(N)) cost.

• cannot def gw(xx, yy, tt, mm, nn, ss): operate on the whole array, rather than individual numbers?

• if denom < 0: can be expressed using boolean mask array like mask = denom >0; denom[mask] = np.sqrt(denom[maks])
• Never do this if you can:
for k in range (Nt):
for j in range (Ny):
for i in range (Nx):

the operations seems to be possible to rewrite that operate on the whole x,y space at once

• This is how you do convolution-like operations in numpy (from here)
# The actual iteration
u[1:-1, 1:-1] = ((u[0:-2, 1:-1] + u[2:, 1:-1])*dy2 +
(u[1:-1,0:-2] + u[1:-1, 2:])*dx2)*dnr_inv
• this is so horrible:
green = np.array([gw(x[i],y[j],t[k],i_grid,j_grid,k_grid) for k_grid in t for j_grid in y for i_grid in x])
• list comprehesion is relatively fast, but still much slower than numpy array operations (which are implemented in C).
• do not create temporary list and convert it to temporary array, you loose lot time doing that.
• why you cannot just prealocate Green[Nx,Ny,:] and Gauss[Nx,Ny,:] and use it as a whole?
• yes, I guess you overdo it a bit. You don't have to work with 6-dimensional arrays right away. Perhaps the best is to keep few outer loops. In this particular case when order of summationd does not matter and the size along each dimension is the same, it does not really matter which loops you keep, in other case you should choose the shorter loops (usually the size of kernel (Green, Gauss) is smaller than size of the function phi). Anyway, as I said, this is convolution, so you can perhaps simply call scipy.ndimage.filters.convolve. – Prokop Hapala Aug 27 at 9:05
• I notice, that the kernel (Green) is translational symmetric, so you can precompute it as 3D arrays depending only of difference `Green[:,:,:]=gw( dx,dy,dz) where dx,dy,dz are 1D arrays representing tt-ss, xx-nn,yy-mm – Prokop Hapala Aug 27 at 9:19
• then you just call phi = scipy.ndimage.filters.convolve( Gauss, Green) – Prokop Hapala Aug 27 at 9:22
• But as I say for such large kernel it is rather inefficient to do it in real space (under the hood the 3D convolve function still does 6-dimensional iteration). Therefore I strongly recommend to do it in Fourier space using FFT if you don't mind periodic boundary conditions. It it can be thousands time faster. Should be something like phi = np.ifftn( np.fftn(Green) * np.fftn(Gauss) ) – Prokop Hapala Aug 27 at 9:28
• when you learn how to use np.fftn and scipy.ndimage.filters.convolve please post the final solution there so other people can use it as a refference. (I don't have time now to test it, and I don't want to post untested code). – Prokop Hapala Aug 27 at 9:32

I'm no mathematician, but here are some suggestions as a long-time programmer:

1. Naming is one of the hardest programming skills to learn, but it is also incredibly important for readability and maintainability. Renaming things can make certain bugs and code smells obvious.

• I would recommend expanding any abbreviations within reason, as long as writing them out makes the code even a smidge more readable, like denominator instead of denom.
• It is not clear what Nx etc. are. It looks like they are half of the extent of the world, since you create x from -Nx to Nx etc. (see below).
• x, y and t are usually used as names for coordinates. In your case they are instead collections of coordinates or axes, so I would name them x_axis or x_coordinates etc.
• xx etc. are coordinates, so I would name them x_coordinate or just x.
• A good rule of thumb is to try to "optimize" comments into names. So rather than

# Define Function to realize Gaussian Disturbance
def g_source(xx, yy, tt):

I would probably write something like

def gaussian_disturbance(x: int, y: int, t: int) -> float:

As far as I can tell this should make at least some of the comments within that function redundant.

2. As far as I can tell np.linspace(-Nx, Nx, 100) creates an array [-100, -98, …, 98, 100]. Is it intentional that points are separated by two, should the range be from 0 to 100, or should num be 200? If the range is meant to be from 0 to 100 that would cut the memory and processing requirements by a factor of eight.

General suggestions:

1. black can automatically format your code to be more idiomatic. In terms of bang for your buck this tool is unbeatable.
2. flake8 with a strict complexity limit will give you more hints to write idiomatic Python:

[flake8]
max-complexity = 4
ignore = W503,E203

It can take some time to understand and correct for all the things it prints out, but doing so will help a lot in the long run.

3. I would then recommend adding type hints and validating them using a strict mypy configuration:

[mypy]
check_untyped_defs = true
disallow_untyped_defs = true
ignore_missing_imports = true
no_implicit_optional = true
warn_redundant_casts = true
warn_return_any = true
warn_unused_ignores = true

This makes it clear what you are passing around, and can highlight issues such as using a variable for more than one thing.

• Thank you so much for the detailed tips and suggestions. I will try to follow these from now on. Appreciate the effort and patience. – Bahauddin Omar Aug 23 at 19:09

Thanks a lot for Prokop Hapala's helpful comments and suggestions. This is the code that has done job.

import numpy as np
from numpy import pi
from scipy import signal

# Assigning grid points

Nx = 100
Ny = 100
Nt = 100

x = np.linspace(-Nx, Nx, 100)
y = np.linspace(-Ny, Ny, 100)
t = np.linspace(0, Nt-1, 100)

# Equation to solve
# phi(x,y,t) = ∫∫∫ G(x,y,t; x',y',t') . Source(x',y',t') dx' dy' dt'
# G(x,y,t; x',y',t') = Green's Function

# phi = displacement by the wave
phi = np.zeros((Nt,Ny,Nx))

# Define Function to realize Green's Function for Wave Equation
def gw(xx, yy, tt):

kk = np.heaviside((tt-np.sqrt(xx**2+yy**2)),1)/(2*np.pi*np.sqrt(np.clip(tt**2-xx**2-yy**2,0,None))+1)
return (kk)

# Define Function to realize Gaussian Disturbance
def g_source(xx, yy, tt):

spatial_sigma  = 5   # spatial width of the Gaussian
temporal_sigma = 30  # temporal width of the Gaussian
onset_time     = 20  # time when Gaussian Disturbance reaches its peak

kk = np.exp((-(np.sqrt((xx)**2+(yy)**2)/spatial_sigma)**2)))*np.exp(-((tt-onset_time)/temporal_sigma)**2)
return (kk)

# Calculate the two function for given grid points
green = gw(x[None,None,:],y[None,:,None],t[:,None,None])
gauss = g_source(x[None,None,:],y[None,:,None],t[:,None,None])

# Calculate Source Response via convolution
phi = signal.convolve(gauss, green, mode='full')

Ideally, one would use scipy.signal.fftconvolve to take the advantage of doing convolution in Fourier space (which brings down the iteration to 3D in this problem, instead of 6D in real space). However, As of v0.19, signal.convolve automatically chooses fft method or the real space method based on an estimation of which is faster.

Note: I tried to employ the scipy.ndimage.filters.convolve for my problem, but it gave back memory error which means it tried to do the 3D convolve function by 6-dimensional iteration under the hood and that was too much for my RAM (32 GB). But the scipy.signal.convolve worked fine even for much denser grid points.

• And it took only 1.01 seconds to solve. Thanks again to Prokop Hapala :D – Bahauddin Omar Aug 28 at 19:10