[cig-commits] r20593 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
cmorency at geodynamics.org
cmorency at geodynamics.org
Sun Aug 19 12:35:20 PDT 2012
Author: cmorency
Date: 2012-08-19 12:35:19 -0700 (Sun, 19 Aug 2012)
New Revision: 20593
Added:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_poroelastic.f90
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
Log:
Poroelastic implementation: debugged poroelastic kernels and added absorbing boundary.
To do: need debbug of backpropagation with absorbing boundary.
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2012-08-19 18:11:05 UTC (rev 20592)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/Makefile.in 2012-08-19 19:35:19 UTC (rev 20593)
@@ -181,6 +181,7 @@
$O/compute_coupling_poroelastic_el.o \
$O/compute_stacey_acoustic.o \
$O/compute_stacey_elastic.o \
+ $O/compute_stacey_poroelastic.o \
$O/compute_gradient.o \
$O/compute_interpolated_dva.o \
$O/initialize_simulation.o \
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-08-19 18:11:05 UTC (rev 20592)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2012-08-19 19:35:19 UTC (rev 20593)
@@ -241,9 +241,11 @@
! adjoint simulations
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
- stop 'adjoint poroelastic simulation not implemented yet'
-! add adjoint source following elastic block by block consideration
+ ! adds adjoint source in this partitions
+ if( nadj_rec_local > 0 ) then
+
+ ! add adjoint source following elastic block by block consideration
! read in adjoint sources block by block (for memory consideration)
! e.g., in exploration experiments, both the number of receivers (nrec) and
! the number of time steps (NSTEP) are huge,
@@ -301,6 +303,9 @@
if (myrank == islice_selected_rec(irec)) then
irec_local = irec_local + 1
+ ispec = ispec_selected_rec(irec)
+ if( ispec_is_poroelastic(ispec) ) then
+
! checks if element is in phase_is_inner run
if (ispec_is_inner(ispec_selected_rec(irec)) .eqv. phase_is_inner) then
@@ -309,7 +314,6 @@
do j = 1,NGLLY
do i = 1,NGLLX
iglob = ibool(i,j,k,ispec_selected_rec(irec))
- iglob = ibool(i,j,k,ispec)
! get poroelastic parameters of current local GLL
phil = phistore(i,j,k,ispec_selected_rec(irec))
rhol_s = rhoarraystore(1,i,j,k,ispec_selected_rec(irec))
@@ -333,11 +337,12 @@
enddo
endif ! phase_is_inner
+ endif ! ispec_is_poroelastic
endif
enddo ! nrec
endif ! it
-
+ endif ! nadj_rec_local
endif !adjoint : if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
! note: b_displ() is read in after Newmark time scheme, thus
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-08-19 18:11:05 UTC (rev 20592)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-08-19 19:35:19 UTC (rev 20593)
@@ -49,8 +49,6 @@
phase_is_inner = .true.
endif
-!Note: Contrary to the elastic & acoustic case, the absorbing implementation is within compute_forces
-!due to the number of material properties which were needed
if( .NOT. GPU_MODE ) then
! solid phase
@@ -87,10 +85,8 @@
num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
phase_ispec_inner_poroelastic )
-
! adjoint simulations: backward/reconstructed wavefield
if( SIMULATION_TYPE == 3 ) then
-stop 'adjoint poroelastic simulation not fully implemented yet'
! solid phase
call compute_forces_solid( iphase, &
@@ -131,6 +127,20 @@
stop 'GPU for poroelastic simulation not implemented'
endif ! GPU_MODE
+! adds poroelastic absorbing boundary terms to accelerations (type Stacey conditions)
+ if(ABSORBING_CONDITIONS) &
+ call compute_stacey_poroelastic(NSPEC_AB,NGLOB_AB,accels_poroelastic,accelw_poroelastic, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ abs_boundary_normal,abs_boundary_jacobian2Dw, &
+ abs_boundary_ijk,abs_boundary_ispec, &
+ num_abs_boundary_faces, &
+ velocs_poroelastic,velocw_poroelastic,rho_vpI,rho_vpII,rho_vsI, &
+ rhoarraystore,phistore,tortstore, &
+ ispec_is_poroelastic,SIMULATION_TYPE,SAVE_FORWARD, &
+ NSTEP,it,NGLOB_ADJOINT,b_accels_poroelastic,b_accelw_poroelastic, &
+ b_num_abs_boundary_faces,b_reclen_field_poro,b_absorb_fields, &
+ b_absorb_fieldw)
+
! acoustic coupling
if( ACOUSTIC_SIMULATION ) then
call compute_coupling_poroelastic_ac(NSPEC_AB,NGLOB_AB, &
@@ -251,14 +261,13 @@
b_request_send_vector_ext_meshw,b_request_recv_vector_ext_meshw)
endif !adjoint
- endif
-
!! DK DK May 2009: removed this because now each slice of a CUBIT + SCOTCH mesh
!! DK DK May 2009: has a different number of spectral elements and therefore
!! DK DK May 2009: only the general non-blocking MPI routines assemble_MPI_vector_ext_mesh_s
!! DK DK May 2009: and assemble_MPI_vector_ext_mesh_w above can be used.
!! DK DK May 2009: For adjoint runs below (SIMULATION_TYPE == 3) they should be used as well.
+ endif
enddo
! solid phase
@@ -318,8 +327,6 @@
if (SIMULATION_TYPE == 3) b_velocw_poroelastic(:,:) = b_velocw_poroelastic(:,:) + &
b_deltatover2*b_accelw_poroelastic(:,:)
-
-
! elastic coupling
if( ELASTIC_SIMULATION ) &
call compute_continuity_disp_po_el(NSPEC_AB,NGLOB_AB,ibool,&
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_poroelastic.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_stacey_poroelastic.f90 2012-08-19 19:35:19 UTC (rev 20593)
@@ -0,0 +1,206 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 2 . 1
+! ---------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA and CNRS / INRIA / University of Pau
+! (c) Princeton University / California Institute of Technology and CNRS / INRIA / University of Pau
+! July 2012
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! for poroelastic solver
+
+! absorbing boundary terms for poroelastic media (type Stacey conditions)
+
+ subroutine compute_stacey_poroelastic(NSPEC_AB,NGLOB_AB,accels,accelw, &
+ ibool,ispec_is_inner,phase_is_inner, &
+ abs_boundary_normal,abs_boundary_jacobian2Dw, &
+ abs_boundary_ijk,abs_boundary_ispec, &
+ num_abs_boundary_faces, &
+ velocs,velocw,rho_vpI,rho_vpII,rho_vsI, &
+ rhoarraystore,phistore,tortstore, &
+ ispec_is_poroelastic,SIMULATION_TYPE,SAVE_FORWARD, &
+ NSTEP,it,NGLOB_ADJOINT,b_accels,b_accelw, &
+ b_num_abs_boundary_faces,b_reclen_field_poro,b_absorb_fields, &
+ b_absorb_fieldw)
+
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accels,accelw
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! type Stacey conditions
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: velocs,velocw
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vpI,rho_vpII,rho_vsI
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: phistore,tortstore
+
+ logical, dimension(NSPEC_AB) :: ispec_is_poroelastic
+
+! absorbing boundary surface
+ integer :: num_abs_boundary_faces
+ real(kind=CUSTOM_REAL) :: abs_boundary_normal(NDIM,NGLLSQUARE,num_abs_boundary_faces)
+ real(kind=CUSTOM_REAL) :: abs_boundary_jacobian2Dw(NGLLSQUARE,num_abs_boundary_faces)
+ integer :: abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces)
+ integer :: abs_boundary_ispec(num_abs_boundary_faces)
+
+! adjoint simulations
+ integer:: SIMULATION_TYPE
+ integer:: NSTEP,it,NGLOB_ADJOINT
+ integer:: b_num_abs_boundary_faces,b_reclen_field_poro
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLLSQUARE,b_num_abs_boundary_faces):: b_absorb_fields, &
+ b_absorb_fieldw
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_ADJOINT):: b_accels,b_accelw
+ logical:: SAVE_FORWARD
+
+! local parameters
+ real(kind=CUSTOM_REAL) vx,vy,vz,vfx,vfy,vfz,nx,ny,nz,tx,ty,tz,tfx,tfy,tfz,vn,vfn,jacobianw
+ integer :: ispec,iglob,i,j,k,iface,igll
+ real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,rhol_bar,phil,tortl
+ !integer:: reclen1,reclen2
+
+ ! checks if anything to do
+ if( num_abs_boundary_faces == 0 ) return
+
+! adjoint simulations:
+ if (SIMULATION_TYPE == 3 .and. num_abs_boundary_faces > 0) then
+ ! reads in absorbing boundary array when first phase is running
+ if( phase_is_inner .eqv. .false. ) then
+ ! note: the index NSTEP-it+1 is valid if b_displ is read in after the Newmark scheme
+ ! uses fortran routine
+ !read(IOABS,rec=NSTEP-it+1) reclen1,b_absorb_field,reclen2
+ !if (reclen1 /= b_reclen_field .or. reclen1 /= reclen2) &
+ ! call exit_mpi(0,'Error reading absorbing contribution b_absorb_field')
+ ! uses c routine for faster reading
+ call read_abs(0,b_absorb_fields,b_reclen_field_poro,NSTEP-it+1)
+ call read_abs(0,b_absorb_fieldw,b_reclen_field_poro,NSTEP-it+1)
+ endif
+ endif !adjoint
+
+
+ ! absorbs absorbing-boundary surface using Stacey condition (Clayton & Enquist)
+ do iface=1,num_abs_boundary_faces
+
+ ispec = abs_boundary_ispec(iface)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ if( ispec_is_poroelastic(ispec) ) then
+
+ ! reference gll points on boundary face
+ do igll = 1,NGLLSQUARE
+
+ ! gets local indices for GLL point
+ i = abs_boundary_ijk(1,igll,iface)
+ j = abs_boundary_ijk(2,igll,iface)
+ k = abs_boundary_ijk(3,igll,iface)
+
+ ! gets velocity
+ iglob=ibool(i,j,k,ispec)
+ vx=velocs(1,iglob)
+ vy=velocs(2,iglob)
+ vz=velocs(3,iglob)
+ !
+ vfx=velocw(1,iglob)
+ vfy=velocw(2,iglob)
+ vfz=velocw(3,iglob)
+
+ ! gets properties
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+
+ ! gets associated normal
+ nx = abs_boundary_normal(1,igll,iface)
+ ny = abs_boundary_normal(2,igll,iface)
+ nz = abs_boundary_normal(3,igll,iface)
+
+ ! velocity component in normal direction (normal points out of element)
+ vn = vx*nx + vy*ny + vz*nz
+ vfn = vfx*nx + vfy*ny + vfz*nz
+
+ ! stacey term: velocity vector component * vp * rho in normal direction + vs * rho component tangential to it
+ tx = rho_vpI(i,j,k,ispec)*vn*nx + rho_vsI(i,j,k,ispec)*(vx-vn*nx)
+ ty = rho_vpI(i,j,k,ispec)*vn*ny + rho_vsI(i,j,k,ispec)*(vy-vn*ny)
+ tz = rho_vpI(i,j,k,ispec)*vn*nz + rho_vsI(i,j,k,ispec)*(vz-vn*nz)
+ !
+ tfx = rho_vpII(i,j,k,ispec)/(phil*rhol_bar)*tortl*rhol_f*vfn*nx - &
+ rho_vsI(i,j,k,ispec)/rhol_bar*rhol_f*(vx-vn*nx)
+ tfy = rho_vpII(i,j,k,ispec)/(phil*rhol_bar)*tortl*rhol_f*vfn*ny - &
+ rho_vsI(i,j,k,ispec)/rhol_bar*rhol_f*(vy-vn*ny)
+ tfz = rho_vpII(i,j,k,ispec)/(phil*rhol_bar)*tortl*rhol_f*vfn*nz - &
+ rho_vsI(i,j,k,ispec)/rhol_bar*rhol_f*(vz-vn*nz)
+
+ ! gets associated, weighted jacobian
+ jacobianw = abs_boundary_jacobian2Dw(igll,iface)
+
+ ! adds stacey term (weak form)
+ accels(1,iglob) = accels(1,iglob) - tx*jacobianw
+ accels(2,iglob) = accels(2,iglob) - ty*jacobianw
+ accels(3,iglob) = accels(3,iglob) - tz*jacobianw
+ !
+ accelw(1,iglob) = accelw(1,iglob) - tfx*jacobianw
+ accelw(2,iglob) = accelw(2,iglob) - tfy*jacobianw
+ accelw(3,iglob) = accelw(3,iglob) - tfz*jacobianw
+
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ b_accels(:,iglob) = b_accels(:,iglob) - b_absorb_fields(:,igll,iface)
+ !
+ b_accelw(:,iglob) = b_accelw(:,iglob) - b_absorb_fieldw(:,igll,iface)
+ else if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then
+ b_absorb_fields(1,igll,iface) = tx*jacobianw
+ b_absorb_fields(2,igll,iface) = ty*jacobianw
+ b_absorb_fields(3,igll,iface) = tz*jacobianw
+ !
+ b_absorb_fieldw(1,igll,iface) = tfx*jacobianw
+ b_absorb_fieldw(2,igll,iface) = tfy*jacobianw
+ b_absorb_fieldw(3,igll,iface) = tfz*jacobianw
+ endif !adjoint
+
+ enddo
+ endif ! ispec_is_poroelastic
+ endif ! ispec_is_inner
+ enddo
+
+ ! adjoint simulations: stores absorbed wavefield part
+ if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. num_abs_boundary_faces > 0 ) then
+ ! writes out absorbing boundary value only when second phase is running
+ if( phase_is_inner .eqv. .true. ) then
+ ! uses fortran routine
+ !write(IOABS,rec=it) b_reclen_field,b_absorb_field,b_reclen_field
+ ! uses c routine
+ call write_abs(0,b_absorb_fields,b_reclen_field_poro,it)
+ call write_abs(0,b_absorb_fieldw,b_reclen_field_poro,it)
+ endif
+ endif
+
+ end subroutine compute_stacey_poroelastic
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2012-08-19 18:11:05 UTC (rev 20592)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/iterate_time.f90 2012-08-19 19:35:19 UTC (rev 20593)
@@ -233,6 +233,11 @@
b_Usolidnorm = maxval(abs(b_potential_dot_dot_acoustic(:)))
endif
endif
+ else
+ if( POROELASTIC_SIMULATION ) then
+ b_Usolidnorm = maxval(sqrt(b_displs_poroelastic(1,:)**2 + b_displs_poroelastic(2,:)**2 + &
+ b_displs_poroelastic(3,:)**2))
+ endif
endif
! check stability of the code, exit if unstable
! negative values can occur with some compilers when the unstable value is greater
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-08-19 18:11:05 UTC (rev 20592)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-08-19 19:35:19 UTC (rev 20593)
@@ -701,7 +701,6 @@
if(FIX_UNDERFLOW_PROBLEM) b_potential_acoustic = VERYSMALLVAL
endif
- endif
! poroelastic domain
if( POROELASTIC_SIMULATION ) then
@@ -728,6 +727,7 @@
if(FIX_UNDERFLOW_PROBLEM) b_displw_poroelastic = VERYSMALLVAL
endif
+ endif
! initialize Moho boundary index
if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
@@ -850,6 +850,66 @@
endif
endif
+
+ ! poroelastic domains
+ if( POROELASTIC_SIMULATION) then
+ ! allocates wavefields for solid and fluid phases
+ allocate(b_absorb_fields(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ allocate(b_absorb_fieldw(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_absorb_fields and b_absorb_fieldw'
+
+ ! size of single record
+ b_reclen_field_poro = CUSTOM_REAL * NDIM * NGLLSQUARE * num_abs_boundary_faces
+
+ ! check integer size limit: size of b_reclen_field must fit onto an
+ ! 4-byte integer
+ if( num_abs_boundary_faces > 2147483647 / (CUSTOM_REAL * NDIM * NGLLSQUARE) ) then
+ print *,'reclen needed exceeds integer 4-byte limit: ',b_reclen_field_poro
+ print *,' ',CUSTOM_REAL, NDIM, NGLLSQUARE, num_abs_boundary_faces
+ print*,'bit size fortran: ',bit_size(b_reclen_field_poro)
+ call exit_MPI(myrank,"error b_reclen_field_poro integer limit")
+ endif
+
+ ! total file size
+ filesize = b_reclen_field_poro
+ filesize = filesize*NSTEP
+
+ if (SIMULATION_TYPE == 3) then
+ ! opens existing files
+
+ ! uses fortran routines for reading
+ !open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='old',&
+ ! action='read',form='unformatted',access='direct', &
+ ! recl=b_reclen_field+2*4,iostat=ier )
+ !if( ier /= 0 ) call exit_mpi(myrank,'error opening
+ !proc***_absorb_field.bin file')
+ ! uses c routines for faster reading
+ call open_file_abs_r(0,trim(prname)//'absorb_fields.bin', &
+ len_trim(trim(prname)//'absorb_fields.bin'), &
+ filesize)
+ call open_file_abs_r(0,trim(prname)//'absorb_fieldw.bin', &
+ len_trim(trim(prname)//'absorb_fieldw.bin'), &
+ filesize)
+
+ else
+ ! opens new file
+ ! uses fortran routines for writing
+ !open(unit=IOABS,file=trim(prname)//'absorb_field.bin',status='unknown',&
+ ! form='unformatted',access='direct',&
+ ! recl=b_reclen_field+2*4,iostat=ier )
+ !if( ier /= 0 ) call exit_mpi(myrank,'error opening
+ !proc***_absorb_field.bin file')
+ ! uses c routines for faster writing (file index 0 for acoutic domain
+ ! file)
+ call open_file_abs_w(0,trim(prname)//'absorb_fields.bin', &
+ len_trim(trim(prname)//'absorb_fields.bin'), &
+ filesize)
+ call open_file_abs_w(0,trim(prname)//'absorb_fieldw.bin', &
+ len_trim(trim(prname)//'absorb_fieldw.bin'), &
+ filesize)
+
+ endif
+ endif
else
! needs dummy array
b_num_abs_boundary_faces = 1
@@ -862,6 +922,12 @@
allocate(b_absorb_potential(NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array b_absorb_potential'
endif
+
+ if( POROELASTIC_SIMULATION ) then
+ allocate(b_absorb_fields(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ allocate(b_absorb_fieldw(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_absorb_fields and b_absorb_fieldw'
+ endif
endif
else ! ABSORBING_CONDITIONS
! needs dummy array
@@ -875,6 +941,12 @@
allocate(b_absorb_potential(NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
if( ier /= 0 ) stop 'error allocating array b_absorb_potential'
endif
+
+ if( POROELASTIC_SIMULATION ) then
+ allocate(b_absorb_fields(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ allocate(b_absorb_fieldw(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_absorb_fields and b_absorb_fieldw'
+ endif
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90 2012-08-19 18:11:05 UTC (rev 20592)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90 2012-08-19 19:35:19 UTC (rev 20593)
@@ -38,7 +38,7 @@
real(kind=CUSTOM_REAL) :: rhol,mul,kappal
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: weights_kernel
real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,rhol_bar,phil,tortl
- real(kind=CUSTOM_REAL) :: kappal_s ! mul_s
+ real(kind=CUSTOM_REAL) :: kappal_s
real(kind=CUSTOM_REAL) :: kappal_f,etal_f
real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr
real(kind=CUSTOM_REAL) :: permlxx,permlxy,permlxz,permlyz,permlyy,permlzz
@@ -147,6 +147,7 @@
D_biot = kappal_s*(1._CUSTOM_REAL + phil*(kappal_s/kappal_f - 1._CUSTOM_REAL))
H_biot = (kappal_s - kappal_fr)*(kappal_s - kappal_fr)/(D_biot - kappal_fr) + &
kappal_fr + 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
+ B_biot = H_biot - 4._CUSTOM_REAL*mul_fr/3._CUSTOM_REAL
C_biot = kappal_s*(kappal_s - kappal_fr)/(D_biot - kappal_fr)
M_biot = kappal_s*kappal_s/(D_biot - kappal_fr)
! Approximated velocities (no viscous dissipation)
@@ -175,9 +176,9 @@
!at the moment suitable for constant permeability
eta_kl(i,j,k,ispec) = - etal_f/permlxx * eta_kl(i,j,k,ispec)
mufr_kl(i,j,k,ispec) = - 2._CUSTOM_REAL * mul_fr * mufr_kl(i,j,k,ispec)
- B_kl(i,j,k,ispec) = - B_kl(i,j,k,ispec)
- C_kl(i,j,k,ispec) = - C_kl(i,j,k,ispec)
- M_kl(i,j,k,ispec) = - M_kl(i,j,k,ispec)
+ B_kl(i,j,k,ispec) = - B_biot * B_kl(i,j,k,ispec)
+ C_kl(i,j,k,ispec) = - C_biot * C_kl(i,j,k,ispec)
+ M_kl(i,j,k,ispec) = - M_biot * M_kl(i,j,k,ispec)
! density kernels
rhob_kl(i,j,k,ispec) = rhot_kl(i,j,k,ispec) + B_kl(i,j,k,ispec) + mufr_kl(i,j,k,ispec)
@@ -384,6 +385,101 @@
endif
+ ! save kernels to binary files
+ if( POROELASTIC_SIMULATION ) then
+ ! primary kernels
+ open(unit=27,file=prname(1:len_trim(prname))//'rhot_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhot_primeporo_kernel.bin'
+ write(27) rhot_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'rhof_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhof_primeporo_kernel.bin'
+ write(27) rhof_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'sm_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file sm_primeporo_kernel.bin'
+ write(27) sm_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'eta_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file eta_primeporo_kernel.bin'
+ write(27) eta_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'mufr_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file mufr_primeporo_kernel.bin'
+ write(27) mufr_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'B_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file B_primeporo_kernel.bin'
+ write(27) B_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'C_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file C_primeporo_kernel.bin'
+ write(27) C_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'M_primeporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file M_primeporo_kernel.bin'
+ write(27) M_kl
+ close(27)
+ ! density kernels
+ open(unit=27,file=prname(1:len_trim(prname))//'rhob_densityporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhob_densityporo_kernel.bin'
+ write(27) rhob_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'rhofb_densityporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhofb_densityporo_kernel.bin'
+ write(27) rhofb_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'phi_densityporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file phi_densityporo_kernel.bin'
+ write(27) phi_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'mufrb_densityporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file mufrb_densityporo_kernel.bin'
+ write(27) mufrb_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'Bb_densityporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file Bb_densityporo_kernel.bin'
+ write(27) Bb_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'Cb_densityporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file Cb_densityporo_kernel.bin'
+ write(27) Cb_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'Mb_densityporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file Mb_densityporo_kernel.bin'
+ write(27) Mb_kl
+ close(27)
+ ! wavespeed kernels
+ open(unit=27,file=prname(1:len_trim(prname))//'rhobb_waveporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhobb_waveporo_kernel.bin'
+ write(27) rhobb_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'rhofbb_waveporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhofbb_waveporo_kernel.bin'
+ write(27) rhofbb_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'phib_waveporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file phib_waveporo_kernel.bin'
+ write(27) phib_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'cs_waveporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file cs_waveporo_kernel.bin'
+ write(27) cs_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'cpI_waveporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file cpI_waveporo_kernel.bin'
+ write(27) cpI_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'cpII_waveporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file cpII_waveporo_kernel.bin'
+ write(27) cpII_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'ratio_waveporo_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file ratio_waveporo_kernel.bin'
+ write(27) ratio_kl
+ close(27)
+ endif
+
! save weights for volume integration, in order to benchmark the kernels with analytical expressions
if( SAVE_WEIGHTS ) then
allocate(weights_kernel(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-08-19 18:11:05 UTC (rev 20592)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-08-19 19:35:19 UTC (rev 20593)
@@ -496,8 +496,8 @@
rhobb_kl, rhofbb_kl, phib_kl, cpI_kl, cpII_kl, cs_kl, ratio_kl
! absorbing stacey wavefield parts
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: b_absorb_fields_poroelastic,b_absorb_fieldw_poroelastic
- integer :: b_reclen_fields,b_reclen_fieldw
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: b_absorb_fields,b_absorb_fieldw
+ integer :: b_reclen_field_poro
! for assembling backward field
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: b_buffer_send_vector_ext_meshs
More information about the CIG-COMMITS
mailing list