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

[xie.zhinan at geodynamics.org]

Author: xie.zhinan
Date: 2013-05-27 04:54:05 -0700 (Mon, 27 May 2013)
New Revision: 22138

Modified:
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
Log:
set potential and its derivatives to be zero on outer boundary of PML to remove some potential instability problem


Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90	2013-05-26 22:10:00 UTC (rev 22137)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90	2013-05-27 11:54:05 UTC (rev 22138)
@@ -941,5 +941,70 @@
     enddo
   endif  ! end of absorbing boundaries
 
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!
+!--- absorbing boundaries
+!
+ if( PML_BOUNDARY_CONDITIONS .and. anyabs) then
+
+! we have to put Dirichlet on the boundary of the PML
+   do ispecabs = 1,nelemabs
+
+     ispec = numabs(ispecabs)
+
+     if (is_PML(ispec)) then
+      ispec_PML=spec_to_PML(ispec)  
+!--- left absorbing boundary
+      if(codeabs(IEDGE4,ispecabs)) then
+        i = 1
+        do j = 1,NGLLZ
+          iglob = ibool(i,j,ispec)
+          potential_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+       enddo
+      endif
+!--- right absorbing boundary
+      if(codeabs(IEDGE2,ispecabs)) then
+        i = NGLLX
+        do j = 1,NGLLZ
+          iglob = ibool(i,j,ispec)
+          potential_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+        enddo
+
+      endif
+!--- bottom absorbing boundary
+      if(codeabs(IEDGE1,ispecabs)) then
+        j = 1
+! exclude corners to make sure there is no contradiction on the normal
+        ibegin = 1
+        iend = NGLLX
+        do i = ibegin,iend
+          iglob = ibool(i,j,ispec)
+          potential_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+        enddo
+      endif
+!--- top absorbing boundary
+      if(codeabs(IEDGE3,ispecabs)) then
+        j = NGLLZ
+! exclude corners to make sure there is no contradiction on the normal
+        ibegin = 1
+        iend = NGLLX
+        do i = ibegin,iend
+          iglob = ibool(i,j,ispec)
+          potential_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_acoustic(iglob) = 0._CUSTOM_REAL
+          potential_dot_dot_acoustic(iglob) = 0._CUSTOM_REAL
+        enddo
+      endif  !  end of top absorbing boundary
+      endif ! end of is_PML
+    enddo ! end specabs loop
+  endif  ! end of absorbing boundaries PML_BOUNDARY_CONDITIONS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
   end subroutine compute_forces_acoustic