[cig-commits] [commit] master: PEP8 all Python files. (31278bf)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Jan 15 23:35:59 PST 2015
Repository : https://github.com/geodynamics/specfem1d
On branch : master
Link : https://github.com/geodynamics/specfem1d/compare/33b512de01491f2c70de206939855d95b22febf4...58c6e68e5fc4e335b2fa135ce471179cc1f417fc
>---------------------------------------------------------------
commit 31278bf92dafaec38280b24ed5b1d9de3fdcd80b
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date: Fri Jan 16 00:32:17 2015 -0500
PEP8 all Python files.
>---------------------------------------------------------------
31278bf92dafaec38280b24ed5b1d9de3fdcd80b
Fortran_version/plot_script_using_python.py | 4 +--
Python_version/config.py | 52 ++++++++++++++++-------------
Python_version/functions.py | 25 ++++++++------
Python_version/gll.py | 1 -
Python_version/grid.py | 14 ++++----
Python_version/specfem1d.py | 31 ++++++++---------
6 files changed, 66 insertions(+), 61 deletions(-)
diff --git a/Fortran_version/plot_script_using_python.py b/Fortran_version/plot_script_using_python.py
index 1cd8b3c..a90811b 100755
--- a/Fortran_version/plot_script_using_python.py
+++ b/Fortran_version/plot_script_using_python.py
@@ -34,8 +34,8 @@ for seismo in args.files:
if not args.hold:
plt.figure()
- plt.plot(data[:,0], data[:,1])
- plt.xlim(data[0,0], data[-1,0])
+ plt.plot(data[:, 0], data[:, 1])
+ plt.xlim(data[0, 0], data[-1, 0])
plt.grid(args.grid)
plt.hold(args.hold)
if not args.hold:
diff --git a/Python_version/config.py b/Python_version/config.py
index e40ad52..8f7a893 100644
--- a/Python_version/config.py
+++ b/Python_version/config.py
@@ -23,13 +23,13 @@ except ImportError:
import numpy as np
import gll
-import functions
class FakeGlobalSectionHead(object):
def __init__(self, fp):
self.fp = fp
self.sechead = '[global]\n'
+
def readline(self):
if self.sechead:
try:
@@ -117,13 +117,17 @@ class Parameter(object):
args = parser.parse_args()
self.plot = self.plot and not args.no_plot
- self.nGLL = self.N + 1 # Number of GLL points per elements
- self.nGLJ = self.NGLJ + 1 # Number of GLJ in the first element
- self.nGlob = (self.nSpec - 1) * self.N + self.NGLJ + 1 # Number of points in the array
- self.ibool = self.make_global_index() # TODO: add GLJ
- self.dt = 0 # Time step (will be updated)
-
- # Gauss Lobatto Legendre points and integration weights :
+ # Number of GLL points per elements
+ self.nGLL = self.N + 1
+ # Number of GLJ in the first element
+ self.nGLJ = self.NGLJ + 1
+ # Number of points in the array
+ self.nGlob = (self.nSpec - 1) * self.N + self.NGLJ + 1
+ self.ibool = self.make_global_index()
+ # Time step (will be updated)
+ self.dt = 0
+
+ # Gauss Lobatto Legendre points and integration weights:
try:
# Position of the GLL points in [-1,1]
self.ksiGLL = gll.GLL_POINTS[self.N]
@@ -145,41 +149,43 @@ class Parameter(object):
def make_global_index(self):
"""Returns a matrix A. A[element_number,GLL_considered] -> index in the
- global array, if we work in axisym and that the element number is 0 the points
- are GLJ points"""
+ global array, if we work in axisym and that the element number is 0 the
+ points are GLJ points"""
+ # TODO: add GLJ
ibool = np.zeros((self.nSpec, self.nGLL), dtype=np.intp)
for e in np.arange(self.nSpec):
for i in np.arange(self.nGLL):
- ibool[e,i] = (self.nGLL - 1) * e + i
+ ibool[e, i] = (self.nGLL - 1) * e + i
return ibool
class Source(object):
"""Contains the source properties"""
- def __init__(self,param):
+ def __init__(self, param):
"""Init"""
- self.typeOfSource=param.sourceType
- self.ampl=param.maxAmpl
+ self.typeOfSource = param.sourceType
+ self.ampl = param.maxAmpl
if self.typeOfSource == 'ricker':
- self.hdur = param.tSource*param.dt # Duration of the source (s)
- self.decayRate = param.decayRate #2.628
- self.alpha = self.decayRate/self.hdur
+ self.hdur = param.tSource * param.dt # Duration of the source (s)
+ self.decayRate = param.decayRate
+ self.alpha = self.decayRate / self.hdur
else:
print "Unknown source's type"
raise
- def __getitem__(self,t):
+ def __getitem__(self, t):
"""What happens when we do source[t]"""
- t-=self.hdur
- return self.ampl*-2.*(self.alpha**3)*t*np.exp(-self.alpha*self.alpha*t*t)/np.sqrt(np.pi)
+ t -= self.hdur
+ return -2 * self.ampl * self.alpha**3 / np.sqrt(np.pi) * (
+ t * np.exp(-(self.alpha * t)**2))
- def plotSource(self,fig=1):
+ def plotSource(self):
"""Plot the source"""
import matplotlib.pyplot as plt
t = np.linspace(0, self.hdur, 1000)
- plt.figure(fig)
- plt.plot(t,self[t],'b')
+ plt.figure()
+ plt.plot(t, self[t], 'b')
plt.title('Source(t)')
plt.grid(True)
plt.show()
diff --git a/Python_version/functions.py b/Python_version/functions.py
index f49b7de..5c849fb 100644
--- a/Python_version/functions.py
+++ b/Python_version/functions.py
@@ -12,10 +12,12 @@ import numpy as np
def estimate_timestep(grid, param):
- dzMin=(grid.z[1:]-grid.z[:len(grid.z)-1]).min()
- dh=param.length/(len(grid.z)-1)
- vMax = np.sqrt(param.meanMu/param.meanRho).max()
- dt = param.cfl*dh/vMax # param.CFL*dzMin/vMax # TODO : See why...?!
+ dh = param.length / (len(grid.z) - 1)
+ vMax = np.sqrt(param.meanMu / param.meanRho)
+ dt = param.cfl * dh / vMax
+ # TODO : See why...?!
+ # dzMin = np.diff(grid.z).min()
+ # dt = param.CFL * dzMin / vMax
return dt
@@ -24,7 +26,8 @@ def project_inverse(ksi, elt_number, ticks):
the number of the element to which it belongs (elt_number) and
the abscissa of the borders of the elements (ticks) returns
its abscissa onto the interval [zmin, zmax]"""
- z0=1./2.*((ksi+1.)*ticks[elt_number+1]+(1.-ksi)*ticks[elt_number])
+ z0 = 1. / 2. * ((ksi + 1) * ticks[elt_number + 1] +
+ (1 - ksi) * ticks[elt_number])
return z0
@@ -41,16 +44,16 @@ def make_stiffness_matrix(grid, param):
ngllj = param.nGLL
deriv = param.deriv
- tmpe = gllj * grid.mu[e,:] * grid.dKsiDx[e,:]**2 * grid.dXdKsi[e,:]
+ tmpe = gllj * grid.mu[e, :] * grid.dKsiDx[e, :]**2 * grid.dXdKsi[e, :]
if param.axisym and e != 0:
tmpe *= project_inverse(param.ksiGLL, e, grid.ticks)
for i in range(ngllj):
- tmpi = deriv[i,:] * tmpe
+ tmpi = deriv[i, :] * tmpe
if param.axisym and e != 0:
tmpi /= project_inverse(param.ksiGLL[i], e, grid.ticks)
- Ke[e,i,:] = np.inner(tmpi, deriv)
+ Ke[e, i, :] = np.inner(tmpi, deriv)
return Ke
@@ -63,10 +66,10 @@ def make_mass_matrix(grid, param):
# tried to slice it in.
for e in range(param.nSpec):
if param.axisym and e == 0:
- Me = param.wGLJ * grid.rho[e,:] * grid.dXdKsi[e,:]
+ Me = param.wGLJ * grid.rho[e, :] * grid.dXdKsi[e, :]
else:
- Me = param.wGLL * grid.rho[e,:] * grid.dXdKsi[e,:]
+ Me = param.wGLL * grid.rho[e, :] * grid.dXdKsi[e, :]
# We can still slice along the second index because the individual GLL
# points are assumed to not repeat a point within an spectral element.
- M[param.ibool[e,:]] += Me
+ M[param.ibool[e, :]] += Me
return M
diff --git a/Python_version/gll.py b/Python_version/gll.py
index 7f755df..9187971 100644
--- a/Python_version/gll.py
+++ b/Python_version/gll.py
@@ -4,7 +4,6 @@ Functions for working with GLL points.
"""
import numpy as np
-from scipy.interpolate import interp1d
# Gauss Lobatto Legendre points and integration weights
diff --git a/Python_version/grid.py b/Python_version/grid.py
index bee1d63..b2692f8 100644
--- a/Python_version/grid.py
+++ b/Python_version/grid.py
@@ -12,13 +12,13 @@ import gll
class OneDimensionalGrid(object):
"""Contains the grid properties"""
- def __init__(self,param):
+ def __init__(self, param):
"""Init"""
- self.param=param
- self.z=np.zeros(param.nGlob)
- self.rho=np.zeros((param.nSpec,param.nGLL))
- self.mu=np.zeros((param.nSpec,param.nGLL))
- self.ticks=np.zeros(param.nSpec+1)
+ self.param = param
+ self.z = np.zeros(param.nGlob)
+ self.rho = np.zeros((param.nSpec, param.nGLL))
+ self.mu = np.zeros((param.nSpec, param.nGLL))
+ self.ticks = np.zeros(param.nSpec + 1)
if param.gridType == 'homogeneous':
self.ticks = np.linspace(0, param.length, param.nSpec + 1)
@@ -47,7 +47,7 @@ class OneDimensionalGrid(object):
elif param.gridType == 'file':
self.z, self.rho, self.mu = np.loadtxt(param.gridFile, unpack=True)
self.ticks = np.loadtxt(param.ticksFile)
- else :
+ else:
print "Unknown grid's type"
raise
# Jacobians at the GLL (and GLJ for the first element in axisym)
diff --git a/Python_version/specfem1d.py b/Python_version/specfem1d.py
index 3585f5f..309e269 100755
--- a/Python_version/specfem1d.py
+++ b/Python_version/specfem1d.py
@@ -19,18 +19,18 @@ import numpy as np
from config import Parameter
from config import Source
-import functions # Contains fundamental functions
+import functions
from grid import OneDimensionalGrid
-### --- MAIN BODY --- ###
+# --- MAIN BODY --- #
# Initialization
-param=Parameter() # Initialize all the parameters of the run. Store them
+param = Parameter()
grid = OneDimensionalGrid(param)
if param.plot:
grid.plot()
param.dt = functions.estimate_timestep(grid, param)
-source=Source(param) # Initialize the source
+source = Source(param)
if param.plot:
source.plotSource()
@@ -41,8 +41,9 @@ Ke = functions.make_stiffness_matrix(grid, param)
M = functions.make_mass_matrix(grid, param)
# Time integration
-u=np.zeros(param.nGlob,dtype='d')
-vel,acc=np.zeros_like(u),np.zeros_like(u)
+u = np.zeros(param.nGlob)
+vel = np.zeros_like(u)
+acc = np.zeros_like(u)
if param.plot:
import matplotlib.pyplot as plt
@@ -55,19 +56,16 @@ if param.axisym and (param.plot or param.snapshot):
# Main time loop :
for it in range(param.nts):
- print 'it = ',it,' (t = ',it*param.dt,'s)'
+ print 'it = ', it, ' (t = ', it * param.dt, 's)'
if it > 0:
u += param.dt * vel + acc * param.dt**2 / 2
vel += param.dt / 2 * acc
acc.fill(0)
for e in range(param.nSpec):
- acc[param.ibool[e,:]] -= np.dot(Ke[e,:,:], u[param.ibool[e,:]])
+ acc[param.ibool[e, :]] -= np.dot(Ke[e, :, :], u[param.ibool[e, :]])
acc[param.iSource] += source[it*param.dt]
-# # Boundary conditions :
-# acc[0] = 0.
-# acc[len(acc)-1] = 0.
acc /= M
vel += param.dt / 2 * acc
@@ -82,13 +80,12 @@ for it in range(param.nts):
if param.plot and it % param.dplot == 0:
if param.axisym:
c = np.concatenate((u[:0:-1], u))
- plt.plot(cz,c)
- plt.xlim([-max(grid.z),max(grid.z)])
- plt.ylim([-10,10]) #plt.ylim([0,0.01])
+ plt.plot(cz, c)
+ plt.xlim([-max(grid.z), max(grid.z)])
+ plt.ylim([-10, 10])
plt.grid(True)
else:
- plt.plot(grid.z,u)
- plt.ylim([-0.10,0.10]) #plt.ylim([0,0.01])
+ plt.plot(grid.z, u)
+ plt.ylim([-0.10, 0.10])
plt.title("it : {}".format(it))
plt.draw()
-
More information about the CIG-COMMITS
mailing list