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

[xie.zhinan at geodynamics.org]

Author: xie.zhinan
Date: 2012-07-11 04:01:53 -0700 (Wed, 11 Jul 2012)
New Revision: 20517

Modified:
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
correct the mass matrix formulation when using PML with mpi


Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90	2012-07-11 10:35:15 UTC (rev 20516)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/invert_mass_matrix.F90	2012-07-11 11:01:53 UTC (rev 20517)
@@ -60,7 +60,7 @@
                                 ,coord &
 #endif
                                 ,K_x_store,K_z_store,npoin_PML,ibool_PML,is_PML,&
-                                d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS)
+                                d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem)
 
 !  builds the global mass matrix
 
@@ -128,6 +128,7 @@
   integer, dimension(NGLLX,NGLLZ,nspec) :: ibool_PML
   logical, dimension(nspec) :: is_PML
   logical :: PML_BOUNDARY_CONDITIONS
+  logical, dimension(4,nspec) :: which_PML_elem
 !! DK DK added this for Guenneau, March 2012
 #ifdef USE_GUENNEAU
   double precision, dimension(NDIM,nglob_elastic), intent(in) :: coord
@@ -182,15 +183,25 @@
 
         if (is_PML(ispec) .and. PML_BOUNDARY_CONDITIONS) then
           iPML=ibool_PML(i,j,ispec)
-!          rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob)  &
-!                  + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(iPML) * K_z_store(iPML)&
-!                  + (d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML)) * deltat / 2.d0)
-
+         if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+             .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
           rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob)  &
-                  + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (1.0d0 &
+                  + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(iPML)&
+                  + d_x_store(iPML) * deltat / 2.d0)
+          rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+         elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+                  (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+          rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob)  &
+                  + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_x_store(iPML) * K_z_store(iPML)&
                   + (d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML)) * deltat / 2.d0)
+          rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+         else
+          rmass_inverse_elastic_one(iglob) = rmass_inverse_elastic_one(iglob)  &
+                  + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * (K_z_store(iPML)&
+                  + d_z_store(iPML)* deltat / 2.d0)
+          rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
+         endif         
 
-          rmass_inverse_elastic_three(iglob) = rmass_inverse_elastic_one(iglob)
         else
 
 !! DK DK added this for Guenneau, March 2012
@@ -220,13 +231,27 @@
 
         if (PML_BOUNDARY_CONDITIONS .and. is_PML(ispec)) then
           iPML=ibool_PML(i,j,ispec)
+
+         if ((which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+             .not. (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
           rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob)  &
+                  + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(iPML)&
+                  + d_x_store(iPML) * deltat / 2.d0)
+         elseif ( (which_PML_elem(ILEFT,ispec) .OR. which_PML_elem(IRIGHT,ispec)) .and. &
+                  (which_PML_elem(ITOP,ispec) .OR. which_PML_elem(IBOTTOM,ispec)) ) then
+          rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob)  &
                   + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_x_store(iPML) * K_z_store(iPML)&
                   + (d_x_store(iPML)*k_z_store(iPML)+d_z_store(iPML)*k_x_store(iPML)) * deltat / 2.d0)
          else
+          rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob)  &
+                  + wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * (K_z_store(iPML)&
+                  + d_z_store(iPML)* deltat / 2.d0)
+         endif  
+
+       else
           rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) &
                   + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
-         endif
+       endif
 
         endif
 

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90	2012-07-11 10:35:15 UTC (rev 20516)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90	2012-07-11 11:01:53 UTC (rev 20517)
@@ -2959,7 +2959,7 @@
                                 ,coord &
 #endif
                                 ,K_x_store,K_z_store,npoin_PML,ibool_PML,is_PML,&
-                                d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS)
+                                d_x_store,d_z_store,PML_BOUNDARY_CONDITIONS,which_PML_elem)
 
 #ifdef USE_MPI
   if ( nproc > 1 ) then