[cig-commits] r19170 - seismo/2D/SPECFEM2D/trunk/src/specfem2D

[dkomati1 at geodynamics.org]

Author: dkomati1
Date: 2011-11-09 14:59:03 -0800 (Wed, 09 Nov 2011)
New Revision: 19170

Modified:
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_color_image.F90
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90
Log:
re-added the corners for initial plane waves (but not for Stacey), but carefully checking that there is no bug this time


Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90	2011-11-09 21:57:45 UTC (rev 19169)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.f90	2011-11-09 22:59:03 UTC (rev 19170)
@@ -662,11 +662,11 @@
 
            j = 1
 
-           ! exclude corners to make sure there is no contradiction on the normal
+!! DK DK not needed           ! exclude corners to make sure there is no contradiction on the normal
            ibegin = 1
            iend = NGLLX
-           if(codeabs(ILEFT,ispecabs)) ibegin = 2
-           if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+!! DK DK not needed           if(codeabs(ILEFT,ispecabs)) ibegin = 2
+!! DK DK not needed           if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
 
            do i = ibegin,iend
 
@@ -725,6 +725,15 @@
                  ty = rho_vs*vy
                  tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
 
+! exclude corners to make sure there is no contradiction on the normal
+! for Stacey absorbing conditions but not for incident plane waves;
+! thus subtract nothing i.e. zero in that case
+                 if((codeabs(ILEFT,ispecabs) .and. i == 1) .or. (codeabs(IRIGHT,ispecabs) .and. i == NGLLX)) then
+                   tx = 0
+                   ty = 0
+                   tz = 0
+                 endif
+
                  accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
                  accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
                  accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
@@ -759,11 +768,11 @@
 
            j = NGLLZ
 
-           ! exclude corners to make sure there is no contradiction on the normal
+!! DK DK not needed           ! exclude corners to make sure there is no contradiction on the normal
            ibegin = 1
            iend = NGLLX
-           if(codeabs(ILEFT,ispecabs)) ibegin = 2
-           if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+!! DK DK not needed           if(codeabs(ILEFT,ispecabs)) ibegin = 2
+!! DK DK not needed           if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
 
            do i = ibegin,iend
 
@@ -814,6 +823,15 @@
                  ty = rho_vs*vy
                  tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
 
+! exclude corners to make sure there is no contradiction on the normal
+! for Stacey absorbing conditions but not for incident plane waves;
+! thus subtract nothing i.e. zero in that case
+                 if((codeabs(ILEFT,ispecabs) .and. i == 1) .or. (codeabs(IRIGHT,ispecabs) .and. i == NGLLX)) then
+                   tx = 0
+                   ty = 0
+                   tz = 0
+                 endif
+
                  accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
                  accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
                  accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_color_image.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_color_image.F90	2011-11-09 21:57:45 UTC (rev 19169)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_color_image.F90	2011-11-09 22:59:03 UTC (rev 19170)
@@ -352,7 +352,7 @@
           endif
         enddo
       enddo
-    endif !if(poroelastic(ispec)) then
+    endif ! of if(poroelastic(ispec)) then
 
 ! now display acoustic layers as constant blue, because they likely correspond to water in the case of ocean acoustics
 ! or in the case of offshore oil industry experiments.

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90	2011-11-09 21:57:45 UTC (rev 19169)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/prepare_initialfield.F90	2011-11-09 22:59:03 UTC (rev 19170)
@@ -402,11 +402,11 @@
     endif
     if(codeabs(IBOTTOM,ispecabs)) then
        j = 1
-       ! exclude corners to make sure there is no contradiction regarding the normal
+!! DK DK not needed       ! exclude corners to make sure there is no contradiction regarding the normal
        ibegin = 1
        iend = NGLLX
-       if(codeabs(ILEFT,ispecabs)) ibegin = 2
-       if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+!! DK DK not needed       if(codeabs(ILEFT,ispecabs)) ibegin = 2
+!! DK DK not needed       if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
        do i = ibegin,iend
           count_bottom=count_bottom+1
           iglob = ibool(i,j,ispec)