[cig-commits] [commit] master: py: Simplify stiffness matrix creation. (a433f54)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/specfem1d

On branch  : master
Link       : https://github.com/geodynamics/specfem1d/compare/2ccfba13481e26296eebf034518cd31d2c727906...33b512de01491f2c70de206939855d95b22febf4

>---------------------------------------------------------------

commit a433f542c8f89b62757b42ebc2374602775f91f2
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date:   Sun Jan 11 03:04:59 2015 -0500

    py: Simplify stiffness matrix creation.
    
    Pull out some loop invariants and use some temporary variables to reduce
    conditions inside loops. Also use some builtin product+sum function to
    get a little speed boost. The optimization here doesn't make a huge
    difference since it's called only once. However, I think it reads a bit
    better than a whole bunch of nested loops.


>---------------------------------------------------------------

a433f542c8f89b62757b42ebc2374602775f91f2
 Python_version/functions.py | 47 ++++++++++++++++++---------------------------
 1 file changed, 19 insertions(+), 28 deletions(-)

diff --git a/Python_version/functions.py b/Python_version/functions.py
index 84f0022..049ec09 100644
--- a/Python_version/functions.py
+++ b/Python_version/functions.py
@@ -30,38 +30,29 @@ def project_inverse(ksi, elt_number, ticks):
 
 def make_stiffness_matrix(grid, param):
     """Computation of stiffness matrices"""
-    Ke=np.zeros((param.nSpec,param.nGLL,param.nGLL),dtype='d')
-    for e in np.arange(param.nSpec):
+    Ke = np.zeros((param.nSpec, param.nGLL, param.nGLL))
+    for e in range(param.nSpec):
         if param.axisym and e == 0:
-            Ke[e,:,:] = _make_stiffness_first_elem(grid, param)
+            gllj = param.wGLJ
+            ngllj = param.nGLJ
+            deriv = param.derivGLJ
         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]*project_inverse(param.ksiGLL[k],e,grid.ticks)/project_inverse(param.ksiGLL[i],e,grid.ticks)
-                    Ke[e,i,j]= sumk
-    return Ke
+            gllj = param.wGLL
+            ngllj = param.nGLL
+            deriv = param.deriv
+
+        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
+            if param.axisym and e != 0:
+                tmpi /= project_inverse(param.ksiGLL[i], e, grid.ticks)
+
+            Ke[e,i,:] = np.inner(tmpi, deriv)
 
-def _make_stiffness_first_elem(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
+    return Ke
 
 
 def make_mass_matrix(grid,param):