In [1]:
from sympy import *
from sympy.abc import *
from sympy.galgebra.ga import *
import numpy as np
from numpy import linalg as LA
from __future__ import print_function
init_printing()
We consder differential operation on a vector $u$ at a given point in 3D with a 1D stencil size $k$ (number of points in the stencil) for every order, the subindex $i$ represent the dimension number $1$ for z, $2$ for $x$ and 3 for $y$,
First order : $ \frac{d u}{dx_i} $
Second order : $ \frac{d^2 u}{dx_i^2} $
Second order cross derivative $ \frac{d^2 u}{dx_i dx_j} $
In [2]:
# Arithmetic operations
k = symbols('k')
s = symbols('s')
# 1D stencil
# multiplication addition
AI_dxi = k + k - 1
AI_dxxi = k + 1 + k - 1
AI_dxxij = (k+1)**2 - 2*k + 1 + 1 + (k+1)**2 - 2*k + 1 -1
In [3]:
# I/O operations
# load
IO_dxi = k
IO_dxxi = k
IO_dxxij = (k+1)**2 - 2*k + 1
In [4]:
# Operational intensity in sing;e precision
print(AI_dxi/(4*IO_dxi))
print(AI_dxxi/(4*IO_dxxi))
print(AI_dxxij/(4*IO_dxxij))
OI_dxi = lambdify(k,AI_dxi/(4*IO_dxi))
OI_dxxi = lambdify(k,AI_dxxi/(4*IO_dxxi))
OI_dxxij = lambdify(k,AI_dxxij/(4*IO_dxxij))
We now consider geophysical wave equations to obtain the theoretical expression of the operational intensity. We write directly the expression of a single time step as a function of differential operators. An operation on a wavefield is counted only once as we consider the minimum of arithmetic operations required.
$ u(x,y,z,t+dt) = dt^2 v^2(x,y,z) ( 2 u(x,y,z,t) + u(x,y,z,t-dt) + \nabla^2 u(x,y,z,t) +q ) $
$ p(x,y,z,t+dt) = dt^2 v^2(x,y,z) \left( 2 p(x,y,z,t) + p(x,y,z,t-dt) +(1+2\epsilon)(\frac{d^2 p(x,t)}{dx^2}+ \frac{d^2 p(x,t)}{dy^2}) + \sqrt{(1+2\delta)} \frac{d^2 r(x,t)}{dz^2} + q \right) $ $ r(x,y,z,t+dt) = dt^2 v^2(x,y,z) \left( 2 r(x,y,z,t) + r(x,y,z,t-dt) +\sqrt{(1+2\delta)}(\frac{d^2 p(x,t)}{dx^2}+ \frac{d^2 p(x,t)}{dy^2}) + \frac{d^2 r(x,t)}{dz^2} + q \right) $
$ p(x,y,z,t+dt) = dt^2 v^2(x,y,z) \left( 2 p(x,y,z,t) + p(x,y,z,t-dt) + (1+2\epsilon) (G_{\bar{x}\bar{x}} + G_{\bar{y}\bar{y}}) p(x,y,z,t) + \sqrt{(1+2\delta)} G_{\bar{z}\bar{z}} r(x,y,z,t) + q \right) $ $ r(x,y,z,t+dt) = dt^2 v^2(x,y,z) \left( 2 r(x,y,z,t) + r(x,y,z,t-dt) + \sqrt{(1+2\delta)}(G_{\bar{x}\bar{x}} + G_{\bar{y}\bar{y}}) p(x,y,z,t) + G_{\bar{z}\bar{z}} r(x,y,z) +q \right) $
where $ \begin{cases} G_{\bar{x}\bar{x}} & = cos(\phi)^2 cos(\theta)^2 \frac{d^2}{dx^2} +sin(\phi)^2 cos(\theta)^2 \frac{d^2}{dy^2}+ sin(\theta)^2 \frac{d^2}{dz^2} + sin(2\phi) cos(\theta)^2 \frac{d^2}{dx dy} - sin(\phi) sin(2\theta) \frac{d^2}{dy dz} -cos(\phi) sin(2\theta) \frac{d^2}{dx dz} \\ G_{\bar{y}\bar{y}} & = sin(\phi)^2 \frac{d^2}{dx^2} +cos(\phi)^2 \frac{d^2}{dy^2} - sin(2\phi)^2 \frac{d^2}{dx dy}\\ G_{\bar{z}\bar{z}} & = cos(\phi)^2 sin(\theta)^2 \frac{d^2}{dx^2} +sin(\phi)^2 sin(\theta)^2 \frac{d^2}{dy^2}+ cos(\theta)^2 \frac{d^2}{dz^2} + sin(2\phi) sin(\theta)^2 \frac{d^2}{dx dy} + sin(\phi) sin(2\theta) \frac{d^2}{dy dz} +cos(\phi) sin(2\theta) \frac{d^2}{dx dz} \\ \end{cases} $
In [5]:
# Arithmetic
# dxi dxxi dxxij multiplications additions duplicates
AI_acou = 0*AI_dxi + 3*AI_dxxi + 0*AI_dxxij + 3 + 5 - 2 * 2
AI_vti = 2 * ( 0*AI_dxi + 3*AI_dxxi + 0*AI_dxxij + 5 + 5 - 2 )
AI_tti = 2 * ( 0*AI_dxi + 3*AI_dxxi + 3*AI_dxxij + 44 + 17 - 0 )
In [6]:
# I/O operations (we load a point once only)
# dxi dxxi dxxij duplicate other load/write
IO_acou = 0*IO_dxi + 3*IO_dxxi + 0*IO_dxxij - 2 + 3
IO_acoum = 0*IO_dxi + 3*IO_dxxi + 0*IO_dxxij - 2 + 3 + 1 / s
IO_vti = 2 * ( 0*IO_dxi + 3*IO_dxxi + 0*IO_dxxij - 2 + 2 ) + 3
IO_tti = 2 * ( 0*IO_dxi + 3*IO_dxxi + 3*IO_dxxij - 3*k +2 + 4 ) + 7
In [7]:
print(simplify(AI_acou/(4*IO_acou)))
print(simplify(AI_vti/(4*IO_vti)))
print(simplify(AI_tti/(4*IO_tti)))
print(simplify(AI_acou/(4*IO_acoum)))
OI_acou = lambdify(k,AI_acou/(4*IO_acou))
OI_acoum = lambdify((k,s),AI_acou/(4*IO_acoum))
OI_vti = lambdify(k,AI_vti/(4*IO_vti))
OI_tti = lambdify(k,AI_tti/(4*IO_tti))
In [8]:
print(limit(OI_acou(k),k,oo))
print(limit(OI_vti(k),k,oo))
print(limit(OI_tti(k),k,oo))
In [ ]:
k=[2,4,6,8,10,12,14,16,18,20,22,24,26,28,30]
OI_wave=np.zeros((15,3))
OI=np.zeros((15,3))
for i in range(0,15):
OI_wave[i,0]=OI_acou(k[i])
OI_wave[i,1]=OI_vti(k[i])
OI_wave[i,2]=OI_tti(k[i])
OI[i,0]=OI_dxi(k[i])
OI[i,1]=OI_dxxi(k[i])
OI[i,2]=OI_dxxij(k[i])
In [ ]:
k=[2,4,6,8,10,12,14,16,18,20,22,24,26,28,30]
s=[2,4,8,16,32,64,128,256]
OI_wavem=np.zeros((15,8))
for i in range(0,15):
for j in range(0,8):
OI_wavem[i,j]=OI_acoum(k[i],s[j])
In [ ]:
import matplotlib.pyplot as plt
In [ ]:
fig = plt.figure()
plt.hold("off")
intensity = plt.plot(OI) # this is how you'd plot a single line...
In [ ]:
fig = plt.figure()
plt.hold("off")
acou = plt.plot(OI_wave[:,0],label='acou') # this is how you'd plot a single line...
acou = plt.plot(OI_wavem[:,5],label='acoum') # this is how you'd plot a single line...
vti = plt.plot(OI_wave[:,1],label='vti') # this is how you'd plot a single line...
tti = plt.plot(OI_wave[:,2],label='tti') # this is how you'd plot a single line...
In [ ]:
fig = plt.figure()
plt.hold("off")
intensity = plt.imshow(OI_wavem) # this is how you'd plot a single line...
plt.show()
In [ ]: