[cig-commits] [commit] : reverted the modifications I made 5 minutes ago because I introduced a problem (hidden dependency) in the algorithm (14d6ad0)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Thu Nov 14 20:16:28 PST 2013


Repository : ssh://geoshell/specfem3d

On branch  : 
Link       : https://github.com/geodynamics/specfem2d/compare/1e201257d91c794056b990a43329e05d04f77454...0000000000000000000000000000000000000000

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

commit 14d6ad01170882435ce4fe3d7e1bc1ebd5faec9c
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date:   Tue Jan 5 18:43:28 2010 +0000

    reverted the modifications I made 5 minutes ago because I introduced a problem (hidden dependency)
    in the algorithm


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

14d6ad01170882435ce4fe3d7e1bc1ebd5faec9c
 specfem3D.f90 | 286 +++++++++++++++++++++++++++++++---------------------------
 1 file changed, 155 insertions(+), 131 deletions(-)

diff --git a/specfem3D.f90 b/specfem3D.f90
index e2d7b5e..e399598 100644
--- a/specfem3D.f90
+++ b/specfem3D.f90
@@ -1560,14 +1560,20 @@
     ispec2D_moho_bot = 0
   endif
 
+  do ispec = 1,NSPEC_AB
+
+    if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
+      if (is_moho_top(ispec)) then
+        ispec2D_moho_top = ispec2D_moho_top + 1
+      else if (is_moho_bot(ispec)) then
+        ispec2D_moho_bot = ispec2D_moho_bot + 1
+      endif
+    endif
+
 !---------------------------------------------------------------------------------------------------
 ! beginning of nested loops on i,j,k to perform the forward calculations in a given element (ispec)
 !---------------------------------------------------------------------------------------------------
 
-!!!!!!!! DK DK beginning of compute_forces_elastic_direct
-
-  do ispec = 1,NSPEC_AB !  for the direct calculation
-
     do k=1,NGLLZ
       do j=1,NGLLY
         do i=1,NGLLX
@@ -1810,139 +1816,20 @@
       tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
       tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
 
-          tempx1l = 0.
-          tempy1l = 0.
-          tempz1l = 0.
-
-          tempx2l = 0.
-          tempy2l = 0.
-          tempz2l = 0.
-
-          tempx3l = 0.
-          tempy3l = 0.
-          tempz3l = 0.
-
-          do l=1,NGLLX
-            fac1 = hprimewgll_xx(l,i)
-            tempx1l = tempx1l + tempx1(l,j,k)*fac1
-            tempy1l = tempy1l + tempy1(l,j,k)*fac1
-            tempz1l = tempz1l + tempz1(l,j,k)*fac1
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ          enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ          do l=1,NGLLY
-            fac2 = hprimewgll_yy(l,j)
-            tempx2l = tempx2l + tempx2(i,l,k)*fac2
-            tempy2l = tempy2l + tempy2(i,l,k)*fac2
-            tempz2l = tempz2l + tempz2(i,l,k)*fac2
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ          enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ          do l=1,NGLLZ
-            fac3 = hprimewgll_zz(l,k)
-            tempx3l = tempx3l + tempx3(i,j,l)*fac3
-            tempy3l = tempy3l + tempy3(i,j,l)*fac3
-            tempz3l = tempz3l + tempz3(i,j,l)*fac3
           enddo
-
-          fac1 = wgllwgll_yz(j,k)
-          fac2 = wgllwgll_xz(i,k)
-          fac3 = wgllwgll_xy(i,j)
-
-! sum contributions from each element to the global mesh
-  iglob = ibool(i,j,k,ispec)
-  accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
-  accel(2,iglob) = accel(2,iglob) - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
-  accel(3,iglob) = accel(3,iglob) - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
-
-! update memory variables based upon the Runge-Kutta scheme
-
-  if(ATTENUATION_VAL .and. not_fully_in_bedrock(ispec)) then
-
-! use Runge-Kutta scheme to march in time
-  do i_sls = 1,N_SLS
-
-! get coefficients for that standard linear solid
-
-  factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
-  alphaval_loc = alphaval(iselected,i_sls)
-  betaval_loc = betaval(iselected,i_sls)
-  gammaval_loc = gammaval(iselected,i_sls)
-
-! term in xx
-  Sn   = factor_loc * epsilondev_xx(i,j,k,ispec)
-  Snp1   = factor_loc * epsilondev_xx_loc(i,j,k)
-  R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in yy
-  Sn   = factor_loc * epsilondev_yy(i,j,k,ispec)
-  Snp1   = factor_loc * epsilondev_yy_loc(i,j,k)
-  R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in zz not computed since zero trace
-! This is because we only implement Q_\mu attenuation and not Q_\kappa.
-! Note that this does *NOT* imply that there is no attenuation for P waves
-! because for Q_\kappa = infinity one gets (see for instance Dahlen and Tromp (1998)
-! equation (9.59) page 350): Q_\alpha = Q_\mu * 3 * (V_p/V_s)^2 / 4
-! therefore Q_\alpha is not zero; for instance for V_p / V_s = sqrt(3)
-! we get Q_\alpha = (9 / 4) * Q_\mu = 2.25 * Q_\mu
-
-! term in xy
-  Sn   = factor_loc * epsilondev_xy(i,j,k,ispec)
-  Snp1   = factor_loc * epsilondev_xy_loc(i,j,k)
-  R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in xz
-  Sn   = factor_loc * epsilondev_xz(i,j,k,ispec)
-  Snp1   = factor_loc * epsilondev_xz_loc(i,j,k)
-  R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in yz
-  Sn   = factor_loc * epsilondev_yz(i,j,k,ispec)
-  Snp1   = factor_loc * epsilondev_yz_loc(i,j,k)
-  R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-    enddo   ! end of loop on memory variables
-
-  endif  !  end attenuation
-
         enddo
       enddo
-    enddo
 
 !---------------------------------------------------------------------------------------------
 ! end of nested loops on i,j,k to perform the forward calculations in a given element (ispec)
 !---------------------------------------------------------------------------------------------
 
-! save deviatoric strain for Runge-Kutta scheme
-  if(ATTENUATION_VAL .and. not_fully_in_bedrock(ispec)) then
-    epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
-    epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
-    epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
-    epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
-    epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
-  endif
-
-  enddo   ! of the spectral element loop for the direct calculation
-
-!!!!!!!! DK DK end of compute_forces_elastic_direct
-
-!!!!!!!! DK DK beginning of compute_forces_elastic_backward
-
 !----------------------------------------------------------------------------------------------------
 ! beginning of nested loops on i,j,k to perform the backward calculations in a given element (ispec)
 !----------------------------------------------------------------------------------------------------
 
     if (SIMULATION_TYPE == 3) then
 
-  do ispec = 1,NSPEC_AB ! for the backward calculation
-
-    if (SAVE_MOHO_MESH) then
-      if (is_moho_top(ispec)) then
-        ispec2D_moho_top = ispec2D_moho_top + 1
-      else if (is_moho_bot(ispec)) then
-        ispec2D_moho_bot = ispec2D_moho_bot + 1
-      endif
-    endif
-
     do k=1,NGLLZ
       do j=1,NGLLY
         do i=1,NGLLX
@@ -2272,6 +2159,136 @@
         b_tempy3(i,j,k) = jacobianl * (b_sigma_xy*gammaxl + b_sigma_yy*gammayl + b_sigma_yz*gammazl)
         b_tempz3(i,j,k) = jacobianl * (b_sigma_xz*gammaxl + b_sigma_yz*gammayl + b_sigma_zz*gammazl)
 
+          enddo
+        enddo
+      enddo
+
+      endif ! of test if SIMULATION_TYPE == 3
+
+!----------------------------------------------------------------------------------------------
+! end of nested loops on i,j,k to perform the backward calculations in a given element (ispec)
+!----------------------------------------------------------------------------------------------
+
+!---------------------------------------------------------------------------------------------------
+! beginning of nested loops on i,j,k to perform the forward calculations in a given element (ispec)
+!---------------------------------------------------------------------------------------------------
+
+    do k=1,NGLLZ
+      do j=1,NGLLY
+        do i=1,NGLLX
+
+          tempx1l = 0.
+          tempy1l = 0.
+          tempz1l = 0.
+
+          tempx2l = 0.
+          tempy2l = 0.
+          tempz2l = 0.
+
+          tempx3l = 0.
+          tempy3l = 0.
+          tempz3l = 0.
+
+          do l=1,NGLLX
+            fac1 = hprimewgll_xx(l,i)
+            tempx1l = tempx1l + tempx1(l,j,k)*fac1
+            tempy1l = tempy1l + tempy1(l,j,k)*fac1
+            tempz1l = tempz1l + tempz1(l,j,k)*fac1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ          enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ          do l=1,NGLLY
+            fac2 = hprimewgll_yy(l,j)
+            tempx2l = tempx2l + tempx2(i,l,k)*fac2
+            tempy2l = tempy2l + tempy2(i,l,k)*fac2
+            tempz2l = tempz2l + tempz2(i,l,k)*fac2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ          enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ          do l=1,NGLLZ
+            fac3 = hprimewgll_zz(l,k)
+            tempx3l = tempx3l + tempx3(i,j,l)*fac3
+            tempy3l = tempy3l + tempy3(i,j,l)*fac3
+            tempz3l = tempz3l + tempz3(i,j,l)*fac3
+          enddo
+
+          fac1 = wgllwgll_yz(j,k)
+          fac2 = wgllwgll_xz(i,k)
+          fac3 = wgllwgll_xy(i,j)
+
+! sum contributions from each element to the global mesh
+  iglob = ibool(i,j,k,ispec)
+  accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
+  accel(2,iglob) = accel(2,iglob) - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
+  accel(3,iglob) = accel(3,iglob) - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
+
+! update memory variables based upon the Runge-Kutta scheme
+
+  if(ATTENUATION_VAL .and. not_fully_in_bedrock(ispec)) then
+
+! use Runge-Kutta scheme to march in time
+  do i_sls = 1,N_SLS
+
+! get coefficients for that standard linear solid
+
+  factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
+  alphaval_loc = alphaval(iselected,i_sls)
+  betaval_loc = betaval(iselected,i_sls)
+  gammaval_loc = gammaval(iselected,i_sls)
+
+! term in xx
+  Sn   = factor_loc * epsilondev_xx(i,j,k,ispec)
+  Snp1   = factor_loc * epsilondev_xx_loc(i,j,k)
+  R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in yy
+  Sn   = factor_loc * epsilondev_yy(i,j,k,ispec)
+  Snp1   = factor_loc * epsilondev_yy_loc(i,j,k)
+  R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in zz not computed since zero trace
+! This is because we only implement Q_\mu attenuation and not Q_\kappa.
+! Note that this does *NOT* imply that there is no attenuation for P waves
+! because for Q_\kappa = infinity one gets (see for instance Dahlen and Tromp (1998)
+! equation (9.59) page 350): Q_\alpha = Q_\mu * 3 * (V_p/V_s)^2 / 4
+! therefore Q_\alpha is not zero; for instance for V_p / V_s = sqrt(3)
+! we get Q_\alpha = (9 / 4) * Q_\mu = 2.25 * Q_\mu
+
+! term in xy
+  Sn   = factor_loc * epsilondev_xy(i,j,k,ispec)
+  Snp1   = factor_loc * epsilondev_xy_loc(i,j,k)
+  R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in xz
+  Sn   = factor_loc * epsilondev_xz(i,j,k,ispec)
+  Snp1   = factor_loc * epsilondev_xz_loc(i,j,k)
+  R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in yz
+  Sn   = factor_loc * epsilondev_yz(i,j,k,ispec)
+  Snp1   = factor_loc * epsilondev_yz_loc(i,j,k)
+  R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+    enddo   ! end of loop on memory variables
+
+  endif  !  end attenuation
+
+        enddo
+      enddo
+    enddo
+
+!---------------------------------------------------------------------------------------------
+! end of nested loops on i,j,k to perform the forward calculations in a given element (ispec)
+!---------------------------------------------------------------------------------------------
+
+!----------------------------------------------------------------------------------------------------
+! beginning of nested loops on i,j,k to perform the backward calculations in a given element (ispec)
+!----------------------------------------------------------------------------------------------------
+
+    if (SIMULATION_TYPE == 3) then
+
+    do k=1,NGLLZ
+      do j=1,NGLLY
+        do i=1,NGLLX
+
           b_tempx1l = 0.
           b_tempy1l = 0.
           b_tempz1l = 0.
@@ -2363,22 +2380,29 @@
       enddo
     enddo
 
+    endif ! of test if SIMULATION_TYPE == 3
+
 !----------------------------------------------------------------------------------------------
 ! end of nested loops on i,j,k to perform the backward calculations in a given element (ispec)
 !----------------------------------------------------------------------------------------------
 
 ! save deviatoric strain for Runge-Kutta scheme
   if(ATTENUATION_VAL .and. not_fully_in_bedrock(ispec)) then
-    b_epsilondev_xx(:,:,:,ispec) = b_epsilondev_xx_loc(:,:,:)
-    b_epsilondev_yy(:,:,:,ispec) = b_epsilondev_yy_loc(:,:,:)
-    b_epsilondev_xy(:,:,:,ispec) = b_epsilondev_xy_loc(:,:,:)
-    b_epsilondev_xz(:,:,:,ispec) = b_epsilondev_xz_loc(:,:,:)
-    b_epsilondev_yz(:,:,:,ispec) = b_epsilondev_yz_loc(:,:,:)
+    epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
+    epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
+    epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
+    epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
+    epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
+    if (SIMULATION_TYPE == 3) then
+      b_epsilondev_xx(:,:,:,ispec) = b_epsilondev_xx_loc(:,:,:)
+      b_epsilondev_yy(:,:,:,ispec) = b_epsilondev_yy_loc(:,:,:)
+      b_epsilondev_xy(:,:,:,ispec) = b_epsilondev_xy_loc(:,:,:)
+      b_epsilondev_xz(:,:,:,ispec) = b_epsilondev_xz_loc(:,:,:)
+      b_epsilondev_yz(:,:,:,ispec) = b_epsilondev_yz_loc(:,:,:)
+    endif
   endif
 
-  enddo   ! of the spectral element loop for the backward calculation
-
-    endif ! of test if SIMULATION_TYPE == 3
+  enddo   ! of the spectral element loop
 
 ! add Stacey conditions
 



More information about the CIG-COMMITS mailing list