[cig-commits] [commit] master: py: Split GLL functions into a separate file. (4573796)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jan 7 16:50:46 PST 2015
Repository : https://github.com/geodynamics/specfem1d
On branch : master
Link : https://github.com/geodynamics/specfem1d/compare/c63d455e0ff0075e43cb487147ab0bc101d4019c...81ba23d88e3baad884dd78c47c7581ffe856f7b2
>---------------------------------------------------------------
commit 4573796a3b519b4f375b09dda02dbf522e459fa2
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date: Wed Jan 7 04:06:32 2015 -0500
py: Split GLL functions into a separate file.
>---------------------------------------------------------------
4573796a3b519b4f375b09dda02dbf522e459fa2
Python_version/defines.py | 62 +++------------
Python_version/functions.py | 91 ----------------------
Python_version/{functions.py => gll.py} | 132 ++++++++++++--------------------
3 files changed, 60 insertions(+), 225 deletions(-)
diff --git a/Python_version/defines.py b/Python_version/defines.py
index 76b5312..1b6f916 100644
--- a/Python_version/defines.py
+++ b/Python_version/defines.py
@@ -22,54 +22,10 @@ except ImportError:
import numpy as np
import matplotlib.pyplot as plt
+import gll
import functions
-# Gauss Lobatto Legendre points and integration weights
-ksiGLL = {
- 4: np.array([-1, -0.6546536707, 0, 0.6546536707, 1]),
- 5: np.array([-1, -0.7650553239, -0.2852315164, 0.2852315164, 0.7650553239,
- 1]),
- 6: np.array([-1, -0.8302238962, -0.4688487934, 0, 0.4688487934,
- 0.8302238962, 1]),
- 7: np.array([-1, -0.8717401485, -0.5917001814, -0.2092992179, 0.2092992179,
- 0.5917001814, 0.8717401485, 1]),
-}
-
-wGLL = {
- 4: np.array([0.1, 0.5444444444, 0.7111111111, 0.5444444444, 0.1]),
- 5: np.array([0.0666666666, 0.3784749562, 0.5548583770, 0.5548583770,
- 0.3784749562, 0.0666666666]),
- 6: np.array([0.0476190476, 0.2768260473, 0.4317453812, 0.4876190476,
- 0.4317453812, 0.2768260473, 0.0476190476]),
- 7: np.array([0.0357142857, 0.2107042271, 0.3411226924, 0.4124587946,
- 0.4124587946, 0.3411226924, 0.2107042271, 0.0357142857]),
-}
-
-# Gauss Lobatto Jacobi points and integration weights
-
-ksiGLJ = {
- 4: np.array([-1.0, -0.5077876295, 0.1323008207, 0.7088201421, 1.0]),
- 5: np.array([-1.0, -0.6507788566, -0.1563704318, 0.3734893787,
- 0.7972962733, 1.0]),
- 6: np.array([-1.0, -0.7401236486, -0.3538526341, 0.09890279315,
- 0.5288423045, 0.8508465697, 1.0]),
- 7: np.array([-1.0, -0.7993818545, -0.4919057913, -0.1117339354,
- 0.2835397724, 0.6337933270, 0.8856884817, 1.0]),
-}
-
-wGLJ = {
- 4: np.array([0.01333333333, 0.2896566946, 0.7360043695, 0.794338936,
- 0.1666666667]),
- 5: np.array([0.006349206349, 0.1503293754, 0.452292685, 0.6858215721,
- 0.5909214468, 0.1142857143]),
- 6: np.array([0.003401360544, 0.08473655296, 0.2803032119, 0.5016469619,
- 0.5945754451, 0.4520031342, 0.08333333333]),
- 7: np.array([0.001984126984, 0.0510689152, 0.1792187805, 0.3533996177,
- 0.4909749105, 0.5047839706, 0.3550776151, 0.06349206349]),
-}
-
-
class FakeGlobalSectionHead(object):
def __init__(self, fp):
self.fp = fp
@@ -160,22 +116,22 @@ class Parameter(object):
# Gauss Lobatto Legendre points and integration weights :
try:
# Position of the GLL points in [-1,1]
- self.ksiGLL = ksiGLL[self.N]
+ self.ksiGLL = gll.GLL_POINTS[self.N]
# Integration weights
- self.wGLL = wGLL[self.N]
+ self.wGLL = gll.GLL_WEIGHTS[self.N]
except KeyError:
raise ValueError('N = %d is invalid!' % (self.N, ))
try:
# Position of the GLJ points in [-1,1]
- self.ksiGLJ = ksiGLJ[self.NGLJ]
+ self.ksiGLJ = gll.GLJ_POINTS[self.NGLJ]
# Integration weights
- self.wGLJ = wGLJ[self.NGLJ]
+ self.wGLJ = gll.GLJ_WEIGHTS[self.NGLJ]
except KeyError:
raise ValueError('NGLJ = %d is invalid!' % (self.NGLJ, ))
# Derivatives of the Lagrange polynomials at the GLL points
- self.deriv = functions.lagrangeDeriv(self.ksiGLL)
- self.derivGLJ = functions.GLJderiv(self.ksiGLJ)
+ self.deriv = gll.lagrange_derivative(self.ksiGLL)
+ self.derivGLJ = gll.glj_derivative(self.ksiGLJ)
class OneDgrid:
@@ -215,8 +171,8 @@ class OneDgrid:
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)
+ self.dXdKsi = gll.jacobian(self.ticks, param)
+ self.dKsiDx = gll.jacobian_inverse(self.ticks, param)
def plotGrid(self,fig=0):
"""Plot the grid
diff --git a/Python_version/functions.py b/Python_version/functions.py
index 270b1d0..a130cbc 100644
--- a/Python_version/functions.py
+++ b/Python_version/functions.py
@@ -9,7 +9,6 @@ elements simulations.
"""
import numpy as np
-from scipy.interpolate import interp1d
def globalArray(nSpec,nGLL):
@@ -37,96 +36,6 @@ def invProjection(ksi, elt_number, ticks):
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 = len(ksiGLL) - 1
- 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')
diff --git a/Python_version/functions.py b/Python_version/gll.py
similarity index 50%
copy from Python_version/functions.py
copy to Python_version/gll.py
index 270b1d0..c6bebfa 100644
--- a/Python_version/functions.py
+++ b/Python_version/gll.py
@@ -1,43 +1,58 @@
# -*- 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)
+Functions for working with GLL points.
"""
import numpy as np
from scipy.interpolate import interp1d
-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
+# Gauss Lobatto Legendre points and integration weights
+GLL_POINTS = {
+ 4: np.array([-1, -0.6546536707, 0, 0.6546536707, 1]),
+ 5: np.array([-1, -0.7650553239, -0.2852315164, 0.2852315164, 0.7650553239,
+ 1]),
+ 6: np.array([-1, -0.8302238962, -0.4688487934, 0, 0.4688487934,
+ 0.8302238962, 1]),
+ 7: np.array([-1, -0.8717401485, -0.5917001814, -0.2092992179, 0.2092992179,
+ 0.5917001814, 0.8717401485, 1]),
+}
+
+GLL_WEIGHTS = {
+ 4: np.array([0.1, 0.5444444444, 0.7111111111, 0.5444444444, 0.1]),
+ 5: np.array([0.0666666666, 0.3784749562, 0.5548583770, 0.5548583770,
+ 0.3784749562, 0.0666666666]),
+ 6: np.array([0.0476190476, 0.2768260473, 0.4317453812, 0.4876190476,
+ 0.4317453812, 0.2768260473, 0.0476190476]),
+ 7: np.array([0.0357142857, 0.2107042271, 0.3411226924, 0.4124587946,
+ 0.4124587946, 0.3411226924, 0.2107042271, 0.0357142857]),
+}
+
+# Gauss Lobatto Jacobi points and integration weights
-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
+GLJ_POINTS = {
+ 4: np.array([-1.0, -0.5077876295, 0.1323008207, 0.7088201421, 1.0]),
+ 5: np.array([-1.0, -0.6507788566, -0.1563704318, 0.3734893787,
+ 0.7972962733, 1.0]),
+ 6: np.array([-1.0, -0.7401236486, -0.3538526341, 0.09890279315,
+ 0.5288423045, 0.8508465697, 1.0]),
+ 7: np.array([-1.0, -0.7993818545, -0.4919057913, -0.1117339354,
+ 0.2835397724, 0.6337933270, 0.8856884817, 1.0]),
+}
-def lagrangeDeriv(ksiGLL):
+GLJ_WEIGHTS = {
+ 4: np.array([0.01333333333, 0.2896566946, 0.7360043695, 0.794338936,
+ 0.1666666667]),
+ 5: np.array([0.006349206349, 0.1503293754, 0.452292685, 0.6858215721,
+ 0.5909214468, 0.1142857143]),
+ 6: np.array([0.003401360544, 0.08473655296, 0.2803032119, 0.5016469619,
+ 0.5945754451, 0.4520031342, 0.08333333333]),
+ 7: np.array([0.001984126984, 0.0510689152, 0.1792187805, 0.3533996177,
+ 0.4909749105, 0.5047839706, 0.3550776151, 0.06349206349]),
+}
+
+
+def lagrange_derivative(ksiGLL):
"""Calculates the values of the derivative of the Lagrange polynomials
at the GLL points"""
N = len(ksiGLL) - 1
@@ -65,7 +80,8 @@ def lagrangeDeriv(ksiGLL):
deriv[i,j]=prod1*sum1
return deriv
-def GLJderiv(ksiGLJ):
+
+def glj_derivative(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
@@ -103,7 +119,8 @@ def GLJderiv(ksiGLJ):
deriv[i,j] = (N*(N+2)-1.)/4.
return deriv
-def jacob(ticks,param):
+
+def jacobian(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
@@ -115,7 +132,8 @@ def jacob(ticks,param):
dXdKsi[e,k] = dE[e]/2 # Here it does not depend on ksi
return dXdKsi
-def jacobInv(ticks,param):
+
+def jacobian_inverse(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
@@ -126,51 +144,3 @@ def jacobInv(ticks,param):
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
More information about the CIG-COMMITS
mailing list