[cig-commits] r21565 - in seismo/3D/SPECFEM3D/trunk: examples/CPML/HOMO8_FREE_SURFACE/DATA examples/CPML/HOMO8_NOFREESURFACE/DATA src/specfem3D
joseph.charles at geodynamics.org
joseph.charles at geodynamics.org
Tue Mar 19 06:32:42 PDT 2013
Author: joseph.charles
Date: 2013-03-19 06:32:42 -0700 (Tue, 19 Mar 2013)
New Revision: 21565
Modified:
seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_FREE_SURFACE/DATA/Par_file
seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_NOFREESURFACE/DATA/Par_file
seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
Log:
fixes one small bug in pml_compute_memory_variables.f90 related to acoustic C-PML elements
Modified: seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_FREE_SURFACE/DATA/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_FREE_SURFACE/DATA/Par_file 2013-03-19 01:13:26 UTC (rev 21564)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_FREE_SURFACE/DATA/Par_file 2013-03-19 13:32:42 UTC (rev 21565)
@@ -23,7 +23,7 @@
# If you use our internal mesher, the only option is 8-node bricks (27-node elements are not supported)
# CUBIT does not support HEX27 elements either (it can generate them, but they are flat, i.e. identical to HEX8).
# To generate HEX27 elements with curvature properly taken into account, you can use Gmsh http://geuz.org/gmsh/
-NGNOD = 8
+NGNOD = 8
# models:
# available options are:
@@ -35,7 +35,7 @@
MODEL = default
# parameters describing the model
-APPROXIMATE_OCEAN_LOAD = .false.
+APPROXIMATE_OCEAN_LOAD = .false.
TOPOGRAPHY = .false.
ATTENUATION = .false.
ANISOTROPY = .false.
@@ -55,7 +55,7 @@
ABSORB_INSTEAD_OF_FREE_SURFACE = .false.
# C-PML boundary conditions for a regional simulation
-PML_CONDITIONS = .true.
+PML_CONDITIONS = .true.
# C-PML top surface
PML_INSTEAD_OF_FREE_SURFACE = .false.
Modified: seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_NOFREESURFACE/DATA/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_NOFREESURFACE/DATA/Par_file 2013-03-19 01:13:26 UTC (rev 21564)
+++ seismo/3D/SPECFEM3D/trunk/examples/CPML/HOMO8_NOFREESURFACE/DATA/Par_file 2013-03-19 13:32:42 UTC (rev 21565)
@@ -23,7 +23,7 @@
# If you use our internal mesher, the only option is 8-node bricks (27-node elements are not supported)
# CUBIT does not support HEX27 elements either (it can generate them, but they are flat, i.e. identical to HEX8).
# To generate HEX27 elements with curvature properly taken into account, you can use Gmsh http://geuz.org/gmsh/
-NGNOD = 8
+NGNOD = 8
# models:
# available options are:
@@ -35,7 +35,7 @@
MODEL = default
# parameters describing the model
-APPROXIMATE_OCEAN_LOAD = .false.
+APPROXIMATE_OCEAN_LOAD = .false.
TOPOGRAPHY = .false.
ATTENUATION = .false.
ANISOTROPY = .false.
@@ -55,7 +55,7 @@
ABSORB_INSTEAD_OF_FREE_SURFACE = .false.
# C-PML boundary conditions for a regional simulation
-PML_CONDITIONS = .true.
+PML_CONDITIONS = .true.
# C-PML top surface
PML_INSTEAD_OF_FREE_SURFACE = .true.
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-03-19 01:13:26 UTC (rev 21564)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/pml_compute_memory_variables.f90 2013-03-19 13:32:42 UTC (rev 21565)
@@ -68,11 +68,10 @@
real(kind=CUSTOM_REAL) :: A18,A19,A20 ! for convolution of strain(simple)
if( CPML_regions(ispec_CPML) == 1 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
!------------------------------------------------------------------------------
!---------------------------- X-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
@@ -317,16 +316,41 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
+
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
endif
-
+ enddo
enddo
enddo
- enddo
else if( CPML_regions(ispec_CPML) == 2 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
!------------------------------------------------------------------------------
!---------------------------- Y-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
@@ -371,12 +395,12 @@
else if( ispec_is_acoustic(ispec) ) then
- rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) &
- + PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_1
- rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = 0.d0
+ rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) = coef0_1 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) &
+ + PML_dpotential_dxl_new(i,j,k,ispec_CPML) * coef1_1 + PML_dpotential_dxl(i,j,k,ispec_CPML) * coef2_1
+ rmemory_dpotential_dxl(i,j,k,ispec_CPML,2) = 0.d0
- dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
- + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
+ dpotentialdxl = A6 * PML_dpotential_dxl(i,j,k,ispec_CPML) &
+ + A7 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,1) + A8 * rmemory_dpotential_dxl(i,j,k,ispec_CPML,2)
endif
!---------------------- A9, A10 and A11 --------------------------
@@ -570,17 +594,41 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
+
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
endif
-
+ enddo
enddo
enddo
- enddo
- else if( CPML_regions(ispec_CPML) == 3 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
+ else if( CPML_regions(ispec_CPML) == 3 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
!------------------------------------------------------------------------------
!---------------------------- Z-surface C-PML ---------------------------------
!------------------------------------------------------------------------------
@@ -823,17 +871,41 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
+
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
endif
-
+ enddo
enddo
enddo
- enddo
else if( CPML_regions(ispec_CPML) == 4 ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
!------------------------------------------------------------------------------
!---------------------------- XY-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
@@ -1115,17 +1187,41 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
+
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
endif
-
+ enddo
enddo
enddo
- enddo
- endif
- if(CPML_regions(ispec_CPML)==5) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ else if( CPML_regions(ispec_CPML) == 5 ) then
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
!------------------------------------------------------------------------------
!---------------------------- XZ-edge C-PML -----------------------------------
@@ -1408,17 +1504,41 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
+
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
endif
-
+ enddo
enddo
enddo
- enddo
else if( CPML_regions(ispec_CPML) == 6 ) then
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
!------------------------------------------------------------------------------
!---------------------------- YZ-edge C-PML -----------------------------------
!------------------------------------------------------------------------------
@@ -1699,17 +1819,41 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
+
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
endif
-
+ enddo
enddo
enddo
- enddo
else if( CPML_regions(ispec_CPML) == 7 ) then
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
!------------------------------------------------------------------------------
!---------------------------- XYZ-corner C-PML --------------------------------
!------------------------------------------------------------------------------
@@ -2144,51 +2288,38 @@
+ A20 * rmemory_dux_dyl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dyl_y(i,j,k,ispec_CPML,2)
duxdxl_y = A19 * PML_dux_dxl(i,j,k,ispec_CPML) &
+ A20 * rmemory_dux_dxl_y(i,j,k,ispec_CPML,1) + rmemory_dux_dxl_y(i,j,k,ispec_CPML,2)
- endif
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
+ sigma_yx = mul*duxdyl_x + mul*duydxl_x
+ sigma_zx = mul*duzdxl_x + mul*duxdzl_x
- enddo
- enddo
- enddo
- endif
+ sigma_xy = mul*duxdyl_y + mul*duydxl_y
+ sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
+ sigma_zy = mul*duzdyl_y + mul*duydzl_y
+ sigma_xz = mul*duzdxl_z + mul*duxdzl_z
+ sigma_yz = mul*duzdyl_z + mul*duydzl_z
+ sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
- if( ispec_is_elastic(ispec) ) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl_x + lambdal*duydyl_x + lambdal*duzdzl_x
- sigma_yx = mul*duxdyl_x + mul*duydxl_x
- sigma_zx = mul*duzdxl_x + mul*duxdzl_x
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
- sigma_xy = mul*duxdyl_y + mul*duydxl_y
- sigma_yy = lambdal*duxdxl_y + lambdalplus2mul*duydyl_y + lambdal*duzdzl_y
- sigma_zy = mul*duzdyl_y + mul*duydzl_y
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
- sigma_xz = mul*duzdxl_z + mul*duxdzl_z
- sigma_yz = mul*duzdyl_z + mul*duydzl_z
- sigma_zz = lambdal*duxdxl_z + lambdal*duydyl_z + lambdalplus2mul*duzdzl_z
-
- ! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl) ! this goes to accel_x
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl) ! this goes to accel_y
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl) ! this goes to accel_z
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl) ! this goes to accel_x
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl) ! this goes to accel_y
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl) ! this goes to accel_z
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl) ! this goes to accel_x
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl) ! this goes to accel_y
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl) ! this goes to accel_z
+ endif
+ enddo
enddo
enddo
- enddo
-
else
-
- stop 'wrong PML flag in PML memory variable calculation routine'
-
+ stop 'wrong PML flag in PML memory variable calculation routine'
endif
end subroutine pml_compute_memory_variables
More information about the CIG-COMMITS
mailing list