[cig-commits] [commit] devel: added the 1D Python version developed by Alexis Bottero (1579d3c)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Dec 3 07:59:10 PST 2014
Repository : https://github.com/geodynamics/specfem1d
On branch : devel
Link : https://github.com/geodynamics/specfem1d/compare/e399d273d20b66094ece6211e0fe1708ab3c3120...3348030bcde1c1c39ea71b701b545cc19ec491d3
>---------------------------------------------------------------
commit 1579d3c0f8db55b4e101b966e2805c5cc6ce6b5b
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Wed Dec 3 16:57:30 2014 +0100
added the 1D Python version developed by Alexis Bottero
>---------------------------------------------------------------
1579d3c0f8db55b4e101b966e2805c5cc6ce6b5b
README_for_Python_version.txt | 5 +
defines.py | 236 ++++++++++++++++++++++++++++++++++++++++++
functions.py | 185 +++++++++++++++++++++++++++++++++
main_for_Python_version.py | 88 ++++++++++++++++
4 files changed, 514 insertions(+)
diff --git a/README_for_Python_version.txt b/README_for_Python_version.txt
new file mode 100644
index 0000000..b69059a
--- /dev/null
+++ b/README_for_Python_version.txt
@@ -0,0 +1,5 @@
+
+Version created by Alexis Bottero, CNRS Marseille, France, 2014.
+
+This python project is still in progress, there are not a lot of comments so far.
+Please contact alexis DOT bottero at gmail dot com if any question.
diff --git a/defines.py b/defines.py
new file mode 100644
index 0000000..395d7d4
--- /dev/null
+++ b/defines.py
@@ -0,0 +1,236 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Nov 13 10:51:04 2013
+
+This file gathers all the constants, classes and parameters used for 1D
+spectral elements simulations.
+This file is the only one that should be edited. Feel free to change the
+constants defined at the beginning.
+For the moment just one type of source (ricker) and one type of medium
+is implemented.
+
+ at author: Alexis Bottero (alexis.bottero at gmail.com)
+"""
+
+### THIS PART CAN BE MODIFIED -> ###
+AXISYM=True
+# Grid description
+GRID_TYPE='homogeneous' # Grid's type
+LENGTH=3000 # "Physical" length of the domain (in meters)
+DENSITY=2500 #kg/m^3
+RIGIDITY=30000000000 #Pa
+GRID_FILE='grid_homogeneous.txt'
+TICKS_FILE='ticks_homogeneous.txt'
+
+# Source description
+TSOURCE=100 #500. # Duration of the source in dt
+ISOURCE=0 #501 #501 #501 #501 #501 # GLL point number on which the source is situated
+MAX_AMPL=1e7 # Maximum amplitude
+SOURCE_TYPE='ricker' # Source's type
+DECAY_RATE=2.628 #10 #2.628
+
+# Other constants
+NSPEC=250 # Number of elements
+N=4 # Degree of the basis functions
+NGLJ=N # Degree of basis functions in the first element
+# For the moment NGLJ need to be = N
+NTS=20000 # Number of time steps
+CFL=0.45 # Courant CFL number
+
+# Plots
+DPLOT=10 # One image is displayed each DPLOT time step
+### <- THIS PART CAN BE MODIFIED ###
+
+### THIS PART IS NOT SUPPOSED TO BE MODIFIED -> ###
+# --- MODULES AND PACKAGES ---
+import numpy as np # NumPy (multidimensional arrays, linear algebra, ...)
+import scipy as sp # SciPy (signal and image processing library)
+import matplotlib as mpl # Matplotlib (2D/3D plotting library)
+import matplotlib.pyplot as plt # Matplotlib's pyplot: MATLAB-like syntax
+from pylab import * # Matplotlib's pylab interface
+
+# --- FUNCTIONS --- #
+import functions # Contains fundamental functions
+
+# Gauss Lobatto Legendre points and integration weights
+ksiGLL4=np.array([-1,-0.6546536707,0,0.6546536707,1])
+wGLL4=np.array([0.1,0.5444444444,0.7111111111,0.5444444444,0.1])
+ksiGLL5=np.array([-1,-0.7650553239,-0.2852315164,0.2852315164,0.7650553239,1])
+wGLL5=np.array([0.0666666666,0.3784749562,0.5548583770,0.5548583770,0.3784749562,0.0666666666])
+ksiGLL6=np.array([-1,-0.8302238962,-0.4688487934,0,0.4688487934,0.8302238962,1])
+wGLL6=np.array([0.0476190476,0.2768260473,0.4317453812,0.4876190476,0.4317453812,0.2768260473,0.0476190476])
+ksiGLL7=np.array([-1,-0.8717401485,-0.5917001814,-0.2092992179,0.2092992179,0.5917001814,0.8717401485,1])
+wGLL7=np.array([0.0357142857,0.2107042271,0.3411226924,0.4124587946,0.4124587946,0.3411226924,0.2107042271,0.0357142857])
+
+if N == 4:
+ ksiGLL=ksiGLL4
+ wGLL=wGLL4
+elif N == 5:
+ ksiGLL=ksiGLL5
+ wGLL=wGLL5
+elif N == 6:
+ ksiGLL=ksiGLL6
+ wGLL=wGLL6
+else:
+ ksiGLL=ksiGLL7
+ wGLL=wGLL7
+
+# Gauss Lobatto Jacobi points and integration weights
+
+ksiGLJ4=np.array([-1.0,-0.5077876295,0.1323008207,0.7088201421,1.0])
+wGLJ4=np.array([0.01333333333,0.2896566946,0.7360043695,0.794338936,0.1666666667])
+ksiGLJ5=np.array([-1.0,-0.6507788566,-0.1563704318,0.3734893787,0.7972962733,1.0])
+wGLJ5=np.array([0.006349206349,0.1503293754,0.452292685,0.6858215721,0.5909214468,0.1142857143])
+ksiGLJ6=np.array([-1.0,-0.7401236486,-0.3538526341,0.09890279315,0.5288423045,0.8508465697,1.0])
+wGLJ6=np.array([0.003401360544,0.08473655296,0.2803032119,0.5016469619,0.5945754451,0.4520031342,0.08333333333])
+ksiGLJ7=np.array([-1.0,-0.7993818545,-0.4919057913,-0.1117339354,0.2835397724,0.6337933270,0.8856884817,1.0])
+wGLJ7=np.array([0.001984126984,0.0510689152,0.1792187805,0.3533996177,0.4909749105,0.5047839706,0.3550776151,0.06349206349])
+
+if NGLJ == 4:
+ ksiGLJ=ksiGLJ4
+ wGLJ=wGLJ4
+elif NGLJ == 5:
+ ksiGLJ=ksiGLJ5
+ wGLJ=wGLJ5
+elif NGLJ == 6:
+ ksiGLJ=ksiGLJ6
+ wGLJ=wGLJ6
+else:
+ ksiGLJ=ksiGLJ7
+ wGLJ=wGLJ7
+
+class Parameter:
+ """Contains all the constants and parameters necessary for 1D spectral
+ element simulation"""
+
+ def __init__(self):
+ """Init"""
+ self.axisym=AXISYM # True if axisymmetry
+ self.length=LENGTH # "Physical" length of the domain (in meters)
+ self.nSpec=NSPEC # Number of elements
+ self.N=N # Degree of the basis functions
+ self.NGLJ=NGLJ # Degree of the basis functions in the first element
+ 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=functions.globalArray(self.nSpec,self.nGLL) # Global array TODO add GLJ
+ self.nts=NTS # Number of time steps
+ self.cfl=CFL # Courant CFL number
+ self.dt=0 # Time step (will be updated)
+ #Grid description
+ self.gridType=GRID_TYPE # Grid's type
+ self.meanRho=DENSITY
+ self.meanMu=RIGIDITY
+ # Source description :
+ self.tSource=TSOURCE # Duration of the source in dt
+ self.iSource=ISOURCE # GLL point number on which the source is situated
+ self.maxAmpl=MAX_AMPL # Maximum amplitude
+ self.sourceType=SOURCE_TYPE # Source's type
+ self.decayRate=DECAY_RATE # Decay rate for the ricker
+ # Gauss Lobatto Legendre points and integration weights :
+ self.ksiGLL=ksiGLL # Position of the GLL points in [-1,1]
+ self.wGLL=wGLL # Integration weights
+ self.ksiGLJ=ksiGLJ # Position of the GLJ points in [-1,1]
+ self.wGLJ=wGLJ # Integration weights
+ # Derivatives of the Lagrange polynomials at the GLL points
+ self.deriv=functions.lagrangeDeriv(self.ksiGLL)
+ self.derivGLJ=functions.GLJderiv(self.ksiGLJ)
+
+class OneDgrid:
+ """Contains the grid properties"""
+
+ def __init__(self,param):
+ """Init"""
+ self.param=param
+ self.gridFile=GRID_FILE
+ self.ticksFile=TICKS_FILE
+ 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=linspace(0,param.length,param.nSpec+1)
+ for e in arange(param.nSpec):
+ for i in np.arange(param.nGLL):
+ self.rho[e,i] = param.meanRho
+ self.mu[e,i] = param.meanMu
+ for i in np.arange(param.nGlob-1)+1:
+ if i < param.nGLJ:
+ self.z[i]=functions.invProjection(param.ksiGLJ[i],0,self.ticks)
+ if i >= param.nGLL and i < param.nGlob-1:
+ self.z[i]=functions.invProjection(param.ksiGLL[i%param.N],i//param.N,self.ticks)
+ else:
+ self.z[param.nGlob-1]=self.ticks[len(self.ticks)-1]
+ elif param.gridType == 'gradient':
+ print "typeOfGrid == 'gradient' Has not been implemented yet"
+ raise
+ elif param.gridType == 'miscellaneous':
+ print "typeOfGrid == 'miscellaneous' Has not been implemented yet"
+ raise
+ elif param.gridType == 'file':
+ [self.z,self.rho,self.mu]=np.loadtxt(defines.GRID_FILE).transpose()
+ self.ticks=np.loadtxt(defines.TICKS_FILE)
+ else :
+ print "Unknown grid's type"
+ raise
+ # Jacobians at the GLL (and GLJ for the first element in axisym)
+ # points (arrays nSpec*(N+1) elements)
+ self.dXdKsi=functions.jacob(self.ticks,param)
+ self.dKsiDx=functions.jacobInv(self.ticks,param)
+
+ def plotGrid(self,fig=0):
+ """Plot the grid
+ my_ticks gives the abscissa of the borders
+ TODO I should test : _the types of the parameters
+ _their sizes"""
+ plt.figure(fig)
+ sub1=plt.subplot(211)
+ plt.hold(True)
+ for e in arange(self.param.nSpec):
+ for i in np.arange(self.param.nGLL):
+ plt.plot(self.z[self.param.ibool[e,i]],self.rho[e,i],'b+')
+ sub1.set_title(r'$\rho(z)$')
+ plt.xticks(self.ticks)
+ plt.grid(True)
+ sub2=plt.subplot(212)
+ for e in arange(self.param.nSpec):
+ for i in np.arange(self.param.nGLL):
+ plt.plot(self.z[self.param.ibool[e,i]],self.mu[e,i],'r+')
+ sub2.set_title(r'$\mu(z)$')
+ plt.suptitle(" Grid ")
+ plt.xticks(self.ticks)
+ plt.grid(True)
+ plt.show()
+ plt.hold(False)
+
+class Source:
+ """Contains the source properties"""
+
+ def __init__(self,param):
+ """Init"""
+ 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
+ else:
+ print "Unknown source's type"
+ raise
+
+ 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)
+
+ def plotSource(self,fig=1):
+ """Plot the source"""
+ t=linspace(0,self.hdur,1000)
+ plt.figure(fig)
+ plt.plot(t,self[t],'b')
+ plt.title('Source(t)')
+ plt.grid(True)
+ plt.show()
+
+### <- THIS PART IS NOT SUPPOSED TO BE MODIFIED ###
+
diff --git a/functions.py b/functions.py
new file mode 100644
index 0000000..5645484
--- /dev/null
+++ b/functions.py
@@ -0,0 +1,185 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Nov 13 10:36:57 2013
+
+This file gathers all the fundamental functions used for 1D spectral
+elements simulations.
+
+ at author: Alexis Bottero (alexis.bottero at gmail.com)
+"""
+
+# --- MODULES AND PACKAGES ---
+import numpy as np # NumPy (multidimensional arrays, linear algebra, ...)
+from numpy import diff
+import scipy as sp # SciPy (signal and image processing library)
+from scipy import misc
+from scipy.interpolate import interp1d
+import matplotlib as mpl # Matplotlib (2D/3D plotting library)
+import matplotlib.pyplot as plt # Matplotlib's pyplot: MATLAB-like syntax
+from pylab import * # Matplotlib's pylab interface
+
+# --- FUNCTIONS ---
+import defines # Contains all the constants and parameters
+
+def globalArray(nSpec,nGLL):
+ """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"""
+ ibool = np.zeros((nSpec,nGLL),dtype='d')
+ for e in np.arange(nSpec):
+ for i in np.arange(nGLL):
+ ibool[e,i] = (nGLL-1)*e+i
+ return ibool
+
+def estimateDt(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 # defines.CFL*dzMin/vMax # TODO : See why...?!
+ return dt
+# z0=((ksi+1)/(1-ksi)*ticks[elt_number+1]+ticks[elt_number])*1/(1+(ksi+1)/(1-ksi))
+def invProjection(ksi, elt_number, ticks):
+ """From the abscissa of a point (ksi) onto the reference interval,
+ 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])
+ return z0
+
+def lagrangeDeriv(ksiGLL):
+ """Calculates the values of the derivative of the Lagrange polynomials
+ at the GLL points"""
+ N=defines.N
+ deriv=np.zeros((N+1,N+1),dtype='d')
+ for i in range(N+1):
+ for j in range(N+1):
+ if i == 0 and j == 0:
+ deriv[i,j] = -N*(N+1.)/4.
+ elif i == N and j == N:
+ deriv[i,j] = N*(N+1.)/4.
+ elif i == j:
+ deriv[i,j] = 0.
+ else:
+ prod1=1.
+ for m in range(N+1):
+ if m!=i:
+ prod1 *= 1/(ksiGLL[i]-ksiGLL[m])
+ sum1=0.
+ for m in range(N+1):
+ prod2=1.
+ for k in range(N+1):
+ if k!=m and k!=i:
+ prod2 *= (ksiGLL[j]-ksiGLL[k])
+ sum1 += prod2
+ deriv[i,j]=prod1*sum1
+ return deriv
+
+def GLJderiv(ksiGLJ):
+ """Calculates the values of the derivative of the polynomials associated
+ to the Gauss-Lobatto-Jacobi quadrature at the GLJ points"""
+ N=len(ksiGLJ)-1
+ deriv=np.zeros((N+1,N+1),dtype='d')
+ lenX=100000
+ xi = np.linspace(-1,1,lenX)
+ xi2=(xi[1:len(xi)]+xi[0:len(xi)-1])/2
+ dxi=xi[1]-xi[0]
+ P = np.polynomial.legendre.legval(xi, np.identity(8))
+ Pb=np.zeros_like(P[N])
+ Pb[1:]=(P[N][1:]+P[N+1][1:])/(1.+xi[1:])
+ Pb[0]=round(Pb[1]) # =+-(N+1)
+ PbInterp = interp1d(xi, Pb)
+ for i in np.arange(N+1):
+ for j in np.arange(N+1):
+ if i == 0 and j == 0:
+ deriv[i,j] = -N*(N+2.)/6.
+ elif i == 0 and 0 < j < N:
+ deriv[i,j] = 2*(-1)**N*PbInterp(ksiGLJ[j])/((1+ksiGLJ[j])*(N+1))
+ elif i == 0 and j == N:
+ deriv[i,j] = (-1)**N/(N+1)
+ elif 0 < i < N and j == 0:
+ deriv[i,j] = (-1)**(N+1)*(N+1)/(2*PbInterp(ksiGLJ[i])*(1+ksiGLJ[i]))
+ elif 0 < i < N and 0 < j < N and i != j:
+ deriv[i,j] = 1/(ksiGLJ[j]-ksiGLJ[i])*PbInterp(ksiGLJ[j])/PbInterp(ksiGLJ[i])
+ elif 0 < i < N and i == j:
+ deriv[i,j] = -1/(2*(1+ksiGLJ[i]))
+ elif 0 < i < N and j == N:
+ deriv[i,j] = 1/(1-ksiGLJ[i])*1/PbInterp(ksiGLJ[i])
+ elif i == N and j == 0:
+ deriv[i,j] = (-1)**(N+1)*(N+1)/4
+ elif i == N and 0 < j < N:
+ deriv[i,j] = -1/(1-ksiGLJ[j])*PbInterp(ksiGLJ[j])
+ elif i == N and j == N:
+ deriv[i,j] = (N*(N+2)-1.)/4.
+ return deriv
+
+def jacob(ticks,param):
+ """Calculates the jacobian dx/dksi of the substitution on the GLL points
+ (and GLJ points for the first element in axisymmetric).
+ Returns a matrix nSpec*(N+1) containing its value for each element and
+ each points"""
+ dE=ticks[1:]-ticks[:len(ticks)-1] # dE[i] : length of element number i
+ dXdKsi=np.zeros((len(dE),len(param.ksiGLL)),dtype='d')
+ for e in np.arange(len(dE)):
+ for k in np.arange(len(param.ksiGLL)):
+ dXdKsi[e,k] = dE[e]/2 # Here it does not depend on ksi
+ return dXdKsi
+
+def jacobInv(ticks,param):
+ """Calculate the jacobian dksi/dx of the substitution on the GLL points
+ (and GLJ points for the first element in axisymmetric).
+ Returns a matrix nSpec*(N+1) containing its value for each element and
+ each points"""
+ dE=ticks[1:]-ticks[:len(ticks)-1] # dE[i] : length of element number i
+ dKsiDx=np.zeros((len(dE),len(param.ksiGLL)),dtype='d')
+ for e in np.arange(len(dE)):
+ for k in np.arange(len(param.ksiGLL)):
+ dKsiDx[e,k] = 2/dE[e] # Here it does not depend on ksi
+ return dKsiDx
+
+def makeStiffness(grid,param):
+ """Computation of stiffness matrices"""
+ Ke=np.zeros((param.nSpec,param.nGLL,param.nGLL),dtype='d')
+ for e in np.arange(param.nSpec):
+ if param.axisym and e == 0 : # The first element
+ Ke[e,:,:]=makeStiffness0(grid,param) # Build the stiffness matrix for first element
+ else:
+ for i in np.arange(param.nGLL):
+ for j in np.arange(param.nGLL):
+ sumk=0
+ for k in np.arange(param.nGLL):
+ if param.axisym is not True:
+ sumk+=param.wGLL[k]*grid.mu[e,k]*param.deriv[i,k]* \
+ param.deriv[j,k]*(grid.dKsiDx[e,k]**2)* \
+ grid.dXdKsi[e,k]
+ else:
+ sumk+=param.wGLL[k]*grid.mu[e,k]*param.deriv[i,k]* \
+ param.deriv[j,k]*(grid.dKsiDx[e,k]**2)* \
+ grid.dXdKsi[e,k]*invProjection(param.ksiGLL[k],e,grid.ticks)/invProjection(param.ksiGLL[i],e,grid.ticks)
+ Ke[e,i,j]= sumk
+ return Ke
+
+def makeStiffness0(grid,param):
+ """Computation of stiffness matrix for the first element"""
+ K0=np.zeros((param.nGLJ,param.nGLJ),dtype='d')
+ for i in np.arange(param.nGLJ):
+ for j in np.arange(param.nGLJ):
+ sumk=0
+ for k in np.arange(param.nGLJ):
+ sumk+=param.wGLJ[k]*grid.mu[0,k]*param.derivGLJ[i,k]* \
+ param.derivGLJ[j,k]*(grid.dKsiDx[0,k]**2)*grid.dXdKsi[0,k]
+ K0[i,j]= sumk
+ return K0
+
+def makeMass(grid,param):
+ """Computation of global mass matrix"""
+ M=np.zeros(param.nGlob,dtype='d')
+ for e in np.arange(param.nSpec):
+ for i in np.arange(param.nGLL):
+ if param.axisym and e != 0:
+ Me=param.wGLL[i]*grid.rho[e,i]*grid.dXdKsi[e,i]#*invProjection(param.ksiGLL[i],e,grid.ticks) # Axisym
+ if param.axisym and e == 0:
+ Me=param.wGLJ[i]*grid.rho[e,i]*grid.dXdKsi[e,i]
+ else:
+ Me=param.wGLL[i]*grid.rho[e,i]*grid.dXdKsi[e,i]#*invProjection(param.ksiGLL[i],e,grid.ticks)
+ M[param.ibool[e,i]]+=Me
+ return M
diff --git a/main_for_Python_version.py b/main_for_Python_version.py
new file mode 100644
index 0000000..84bc869
--- /dev/null
+++ b/main_for_Python_version.py
@@ -0,0 +1,88 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Nov 13 10:30:04 2013
+
+Main script for 1D spectral elements.
+ u(z,t) = sum(ui(t)li(z)) for i=1...nSpec*(N+1)
+ with nSpec : number of elements:
+ N+1 : number of basis functions per element
+ li(z): Lagrange polynomial of degree N
+ -> ui(t): Time dependent expansion coeffs (what we are looking for)
+
+ The parameters of the run can be edited in the file defines.py
+
+ at author: Alexis Bottero, CNRS Marseille, France (alexis.bottero at gmail.com)
+"""
+
+### --- MODULES AND PACKAGES --- ###
+import numpy as np # NumPy (multidimensional arrays, linear algebra, ...)
+import scipy as sp # SciPy (signal and image processing library)
+import matplotlib as mpl # Matplotlib (2D/3D plotting library)
+import matplotlib.pyplot as plt # Matplotlib's pyplot: MATLAB-like syntax
+from pylab import * # Matplotlib's pylab interface
+
+### --- FUNCTIONS --- ###
+import defines # Contains all the constants and parameters
+from defines import Parameter
+from defines import OneDgrid
+from defines import Source
+import functions # Contains fundamental functions
+
+### --- MAIN BODY --- ###
+# Initialization
+param=Parameter() # Initialize all the parameters of the run. Store them
+grid=OneDgrid(param) # Initialize the grid from these parameters
+#grid.plotGrid() # Plot the grid
+
+param.dt=functions.estimateDt(grid,param) # estimate the time step in seconds
+source=Source(param) # Initialize the source
+#source.plotSource() # Plot the source
+
+# Computation of stiffness matrices
+Ke=functions.makeStiffness(grid,param)
+
+# Computation of global mass matrices
+M=functions.makeMass(grid,param)
+
+# Time integration
+u=np.zeros(param.nGlob,dtype='d')
+vel,acc=np.zeros_like(u),np.zeros_like(u)
+
+plt.ion()
+fig=plt.figure()
+plt.hold(False)
+bz=-np.array([i for i in reversed(grid.z)])
+cz=append(bz,grid.z)
+
+# Main time loop :
+for it in np.arange(param.nts):
+ 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[:] = 0.
+
+ for e in np.arange(param.nSpec):
+ acc[list(param.ibool[e,:])] -= np.dot(Ke[e,:,:],u[list(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[:]
+
+ if it%defines.DPLOT == 0:
+ if param.axisym:
+ b=np.array([i for i in reversed(u)])
+ c=append(b,u)
+ plt.plot(cz,c)
+ plt.xlim([-max(grid.z),max(grid.z)])
+ plt.ylim([-10,10]) #plt.ylim([0,0.01])
+ plt.grid(True)
+ else:
+ plt.plot(grid.z,u)
+ plt.ylim([-0.10,0.10]) #plt.ylim([0,0.01])
+ plt.title("it : {}".format(it))
+ plt.draw()
+
More information about the CIG-COMMITS
mailing list