[cig-commits] r20525 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
cmorency at geodynamics.org
cmorency at geodynamics.org
Sun Jul 15 08:31:06 PDT 2012
Author: cmorency
Date: 2012-07-15 08:31:06 -0700 (Sun, 15 Jul 2012)
New Revision: 20525
Modified:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_fluid.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_solid.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_kernels.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
Log:
Adjoint poroelastic implementation : defined poroelastic kernels (primary, density & wavespeed), saved kernels.
To do: check treatment of absorbing boundary.
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_fluid.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_fluid.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_fluid.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -34,6 +34,8 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
phistore,tortstore,jacobian,ibool,&
+ epsilonwdev_xx,epsilonwdev_yy,epsilonwdev_xy,&
+ epsilonwdev_xz,epsilonwdev_yz,epsilonw_trace_over_3, &
SIMULATION_TYPE,NSPEC_ADJOINT, &
num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
phase_ispec_inner_poroelastic )
@@ -53,6 +55,10 @@
velocw_poroelastic
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: &
+ epsilonwdev_xx,epsilonwdev_yy,epsilonwdev_xy,epsilonwdev_xz,epsilonwdev_yz
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: epsilonw_trace_over_3
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -330,7 +336,12 @@
sigmap = C_biot*duxdxl_plus_duydyl_plus_duzdzl + M_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
if(SIMULATION_TYPE == 3) then ! kernels calculation
- l = NSPEC_ADJOINT ! to avoid compilation warning
+ epsilonw_trace_over_3(i,j,k,ispec) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilonwdev_xx(i,j,k,ispec) = duxdxl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilonwdev_yy(i,j,k,ispec) = duydyl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilonwdev_xy(i,j,k,ispec) = 0.5 * duxdyl_plus_duydxl
+ epsilonwdev_xz(i,j,k,ispec) = 0.5 * duzdxl_plus_duxdzl
+ epsilonwdev_yz(i,j,k,ispec) = 0.5 * duzdyl_plus_duydzl
endif
! endif !if(VISCOATTENUATION)
@@ -457,8 +468,6 @@
accelw_poroelastic(3,iglob) = accelw_poroelastic(3,iglob) + ( fac1*(rhol_f/rhol_bar*tempz1ls - tempz1lw) &
+ fac2*(rhol_f/rhol_bar*tempz2ls - tempz2lw) + fac3*(rhol_f/rhol_bar*tempz3ls - tempz3lw) )
- if(SIMULATION_TYPE == 3) then ! kernels calculation
- endif
!
!---- viscous damping
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -50,8 +50,9 @@
endif
!Note: Contrary to the elastic & acoustic case, the absorbing implementation is within compute_forces
-!due to the number of properties which were needed
+!due to the number of material properties which were needed
+ if( .NOT. GPU_MODE ) then
! solid phase
call compute_forces_solid( iphase, &
@@ -63,6 +64,8 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
phistore,tortstore,jacobian,ibool,&
+ epsilonsdev_xx,epsilonsdev_yy,epsilonsdev_xy,&
+ epsilonsdev_xz,epsilonsdev_yz,epsilons_trace_over_3, &
SIMULATION_TYPE,NSPEC_ADJOINT, &
num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
phase_ispec_inner_poroelastic )
@@ -78,6 +81,8 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
phistore,tortstore,jacobian,ibool,&
+ epsilonwdev_xx,epsilonwdev_yy,epsilonwdev_xy,&
+ epsilonwdev_xz,epsilonwdev_yz,epsilonw_trace_over_3, &
SIMULATION_TYPE,NSPEC_ADJOINT, &
num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
phase_ispec_inner_poroelastic )
@@ -97,6 +102,8 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
phistore,tortstore,jacobian,ibool,&
+ b_epsilonsdev_xx,b_epsilonsdev_yy,b_epsilonsdev_xy,&
+ b_epsilonsdev_xz,b_epsilonsdev_yz,b_epsilons_trace_over_3, &
SIMULATION_TYPE,NSPEC_ADJOINT, &
num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
phase_ispec_inner_poroelastic )
@@ -112,11 +119,17 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
phistore,tortstore,jacobian,ibool,&
+ b_epsilonwdev_xx,b_epsilonwdev_yy,b_epsilonwdev_xy,&
+ b_epsilonwdev_xz,b_epsilonwdev_yz,b_epsilonw_trace_over_3, &
SIMULATION_TYPE,NSPEC_ADJOINT, &
num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
phase_ispec_inner_poroelastic )
endif
+ else
+ ! on GPU
+stop 'GPU for poroelastic simulation not implemented'
+ endif ! GPU_MODE
! acoustic coupling
if( ACOUSTIC_SIMULATION ) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_solid.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_solid.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_solid.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -34,6 +34,8 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
phistore,tortstore,jacobian,ibool,&
+ epsilonsdev_xx,epsilonsdev_yy,epsilonsdev_xy,&
+ epsilonsdev_xz,epsilonsdev_yz,epsilons_trace_over_3, &
SIMULATION_TYPE,NSPEC_ADJOINT, &
num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
phase_ispec_inner_poroelastic )
@@ -52,6 +54,10 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic,accels_poroelastic
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displw_poroelastic,velocw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: &
+ epsilonsdev_xx,epsilonsdev_yy,epsilonsdev_xy,epsilonsdev_xz,epsilonsdev_yz
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: epsilons_trace_over_3
+
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
@@ -332,31 +338,12 @@
! endif !if(ATTENUATION)
if(SIMULATION_TYPE == 3) then ! kernels calculation
- l = NSPEC_ADJOINT ! to avoid compilation warnings
-! kernels calculation
-! if(isolver == 2) then
-! iglob = ibool(i,j,ispec)
-! dsxx = dux_dxl
-! dsxz = HALF * (duz_dxl + dux_dzl)
-! dszz = duz_dzl
-!
-! dwxx = dwx_dxl
-! dwxz = HALF * (dwz_dxl + dwx_dzl)
-! dwzz = dwz_dzl
-!
-! b_dsxx = b_dux_dxl
-! b_dsxz = HALF * (b_duz_dxl + b_dux_dzl)
-! b_dszz = b_duz_dzl
-!
-! b_dwxx = b_dwx_dxl
-! b_dwxz = HALF * (b_dwz_dxl + b_dwx_dzl)
-! b_dwzz = b_dwz_dzl
-!
-! B_k(iglob) = (dux_dxl + duz_dzl) * (b_dux_dxl + b_duz_dzl) * (H_biot - FOUR_THIRDS * mul_fr)
-! mufr_k(iglob) = (dsxx * b_dsxx + dszz * b_dszz + &
-! 2._CUSTOM_REAL * dsxz * b_dsxz - &
-! 1._CUSTOM_REAL/3._CUSTOM_REAL * (dux_dxl + duz_dzl) * (b_dux_dxl + b_duz_dzl) ) * mul_fr
-! endif
+ epsilons_trace_over_3(i,j,k,ispec) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilonsdev_xx(i,j,k,ispec) = duxdxl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilonsdev_yy(i,j,k,ispec) = duydyl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilonsdev_xy(i,j,k,ispec) = 0.5 * duxdyl_plus_duydxl
+ epsilonsdev_xz(i,j,k,ispec) = 0.5 * duzdxl_plus_duxdzl
+ epsilonsdev_yz(i,j,k,ispec) = 0.5 * duzdyl_plus_duydzl
endif
@@ -391,8 +378,6 @@
tempy3p(i,j,k) = jacobianl * sigmap*gammayl
tempz3p(i,j,k) = jacobianl * sigmap*gammazl
- if(SIMULATION_TYPE == 3) then ! kernels calculation
- endif
enddo
enddo
@@ -480,9 +465,6 @@
accels_poroelastic(3,iglob) = accels_poroelastic(3,iglob) - ( fac1*(tempz1ls - phil/tortl*tempz1lw) &
+ fac2*(tempz2ls - phil/tortl*tempz2lw) + fac3*(tempz3ls - phil/tortl*tempz3lw) )
- if(SIMULATION_TYPE == 3) then ! kernels calculation
- endif
-
!
!---- viscous damping
!
@@ -523,7 +505,7 @@
bl_relaxed(5) = etal_f*invpermlyz
bl_relaxed(6) = etal_f*invpermlzz
-! ifVISCOATTENUATION) then
+! if(VISCOATTENUATION) then
! bl_unrelaxed(1) = etal_f*invpermlxx*theta_e/theta_s
! bl_unrelaxed(2) = etal_f*invpermlxz*theta_e/theta_s
! bl_unrelaxed(3) = etal_f*invpermlzz*theta_e/theta_s
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_kernels.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_kernels.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -27,11 +27,13 @@
subroutine compute_kernels()
! kernel calculations
-! see e.g. Tromp et al. (2005)
+! see e.g. Tromp et al. (2005) for elastic calculation
+! and Morency et al. (2009) for poroelastic calculation
use specfem_par
use specfem_par_elastic
use specfem_par_acoustic
+ use specfem_par_poroelastic
implicit none
! elastic simulations
@@ -44,6 +46,11 @@
call compute_kernels_ac()
endif
+ ! poroelastic simulations
+ if( POROELASTIC_SIMULATION ) then
+ call compute_kernels_po()
+ endif
+
! computes an approximative hessian for preconditioning kernels
if ( APPROXIMATE_HESS_KL ) then
call compute_kernels_hessian()
@@ -235,6 +242,112 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine compute_kernels_po()
+
+! kernel calculations
+! see e.g. Morency et al. (2009)
+
+ use specfem_par
+ use specfem_par_poroelastic
+
+ implicit none
+ ! local parameters
+ integer :: i,j,k,ispec,iglob
+ real(kind=CUSTOM_REAL), dimension(5) :: epsilonsdev_loc,b_epsilonsdev_loc
+
+ if( .not. GPU_MODE ) then
+ ! updates kernels on CPU
+ do ispec = 1, NSPEC_AB
+
+ ! poroelastic domains
+ if( ispec_is_poroelastic(ispec) ) then
+
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+ iglob = ibool(i,j,k,ispec)
+
+ ! 8 isotropic kernels
+ ! note: takes displacement from backward/reconstructed
+ ! (forward) field b_displ
+ ! and acceleration from adjoint field accel
+ ! (containing adjoint sources)
+ !
+ ! note: : time integral summation uses deltat
+ ! 3 density kernels, see Morency et al. (2009),
+ ! equations (39)-(41)
+ rhot_kl(i,j,k,ispec) = rhot_kl(i,j,k,ispec) &
+ + deltat * dot_product(accels_poroelastic(:,iglob), b_displs_poroelastic(:,iglob))
+
+ rhof_kl(i,j,k,ispec) = rhof_kl(i,j,k,ispec) &
+ + deltat * ( dot_product(accelw_poroelastic(:,iglob), b_displs_poroelastic(:,iglob)) &
+ + dot_product(accels_poroelastic(:,iglob), b_displw_poroelastic(:,iglob)) )
+
+ sm_kl(i,j,k,ispec) = sm_kl(i,j,k,ispec) &
+ + deltat * dot_product(accelw_poroelastic(:,iglob), b_displw_poroelastic(:,iglob))
+
+ ! kernel for viscous damping, see e.g. Morency et al. (2009),
+ ! equation (42)
+ eta_kl(i,j,k,ispec) = eta_kl(i,j,k,ispec) &
+ + deltat * dot_product(velocw_poroelastic(:,iglob), b_displw_poroelastic(:,iglob))
+
+ ! kernel for frame shear modulus, see e.g. Morency et al. (2009),
+ ! equation (46)
+ ! note: multiplication with 2*mufr(x) will be done after the
+ ! time loop
+ epsilonsdev_loc(1) = epsilonsdev_xx(i,j,k,ispec)
+ epsilonsdev_loc(2) = epsilonsdev_yy(i,j,k,ispec)
+ epsilonsdev_loc(3) = epsilonsdev_xy(i,j,k,ispec)
+ epsilonsdev_loc(4) = epsilonsdev_xz(i,j,k,ispec)
+ epsilonsdev_loc(5) = epsilonsdev_yz(i,j,k,ispec)
+
+ b_epsilonsdev_loc(1) = b_epsilonsdev_xx(i,j,k,ispec)
+ b_epsilonsdev_loc(2) = b_epsilonsdev_yy(i,j,k,ispec)
+ b_epsilonsdev_loc(3) = b_epsilonsdev_xy(i,j,k,ispec)
+ b_epsilonsdev_loc(4) = b_epsilonsdev_xz(i,j,k,ispec)
+ b_epsilonsdev_loc(5) = b_epsilonsdev_yz(i,j,k,ispec)
+
+ mufr_kl(i,j,k,ispec) = mufr_kl(i,j,k,ispec) &
+ + deltat * (epsilonsdev_loc(1)*b_epsilonsdev_loc(1) + epsilonsdev_loc(2)*b_epsilonsdev_loc(2) &
+ + (epsilonsdev_loc(1)+epsilonsdev_loc(2)) * (b_epsilonsdev_loc(1)+b_epsilonsdev_loc(2)) &
+ + 2 * (epsilonsdev_loc(3)*b_epsilonsdev_loc(3) + epsilonsdev_loc(4)*b_epsilonsdev_loc(4) + &
+ epsilonsdev_loc(5)*b_epsilonsdev_loc(5)) )
+
+ ! 3 kernels for the bulk moduli, see e.g. Morency et al. (2009),
+ ! equations (43)-(45)
+ ! note: multiplication by respective moduli will be done after the
+ ! time loop
+ B_kl(i,j,k,ispec) = B_kl(i,j,k,ispec) &
+ + deltat * (9 * epsilons_trace_over_3(i,j,k,ispec) &
+ * b_epsilons_trace_over_3(i,j,k,ispec))
+
+ C_kl(i,j,k,ispec) = C_kl(i,j,k,ispec) &
+ + deltat * (9 * epsilons_trace_over_3(i,j,k,ispec) &
+ * b_epsilonw_trace_over_3(i,j,k,ispec) &
+ + 9 * epsilonw_trace_over_3(i,j,k,ispec) &
+ * b_epsilons_trace_over_3(i,j,k,ispec) )
+
+ M_kl(i,j,k,ispec) = M_kl(i,j,k,ispec) &
+ + deltat * (9 * epsilonw_trace_over_3(i,j,k,ispec) &
+ * b_epsilonw_trace_over_3(i,j,k,ispec))
+ enddo
+ enddo
+ enddo
+ endif !ispec_is_poroelastic
+
+ enddo
+
+ else
+ ! TO DO
+ ! updates kernels on GPU
+ !call compute_kernels_poroelastic_cuda(Mesh_pointer,deltat)
+ endif
+
+ end subroutine compute_kernels_po
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine compute_kernels_hessian()
use specfem_par
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/finalize_simulation.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -31,6 +31,7 @@
use specfem_par
use specfem_par_elastic
use specfem_par_acoustic
+ use specfem_par_poroelastic
implicit none
@@ -67,6 +68,15 @@
endif
endif
+ if( POROELASTIC_SIMULATION ) then
+ write(27) displs_poroelastic
+ write(27) velocs_poroelastic
+ write(27) accels_poroelastic
+ write(27) displw_poroelastic
+ write(27) velocw_poroelastic
+ write(27) accelw_poroelastic
+ endif
+
close(27)
! adjoint simulations
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.F90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -703,6 +703,32 @@
endif
endif
+ ! poroelastic domain
+ if( POROELASTIC_SIMULATION ) then
+ rhot_kl(:,:,:,:) = 0._CUSTOM_REAL
+ rhof_kl(:,:,:,:) = 0._CUSTOM_REAL
+ sm_kl(:,:,:,:) = 0._CUSTOM_REAL
+ eta_kl(:,:,:,:) = 0._CUSTOM_REAL
+ mufr_kl(:,:,:,:) = 0._CUSTOM_REAL
+ B_kl(:,:,:,:) = 0._CUSTOM_REAL
+ C_kl(:,:,:,:) = 0._CUSTOM_REAL
+ M_kl(:,:,:,:) = 0._CUSTOM_REAL
+
+ !if ( APPROXIMATE_HESS_KL ) &
+ ! hess_kl(:,:,:,:) = 0._CUSTOM_REAL
+
+ ! reconstructed/backward elastic wavefields
+ b_displs_poroelastic = 0._CUSTOM_REAL
+ b_velocs_poroelastic = 0._CUSTOM_REAL
+ b_accels_poroelastic = 0._CUSTOM_REAL
+ b_displw_poroelastic = 0._CUSTOM_REAL
+ b_velocw_poroelastic = 0._CUSTOM_REAL
+ b_accelw_poroelastic = 0._CUSTOM_REAL
+ if(FIX_UNDERFLOW_PROBLEM) b_displs_poroelastic = VERYSMALLVAL
+ if(FIX_UNDERFLOW_PROBLEM) b_displw_poroelastic = VERYSMALLVAL
+
+ endif
+
! initialize Moho boundary index
if (SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3) then
ispec2D_moho_top = 0
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -261,6 +261,24 @@
rho_vsI(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array poroelastic properties'
+ ! needed for kernel computations
+ allocate(epsilonsdev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonsdev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonsdev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonsdev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonsdev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonwdev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonwdev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonwdev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonwdev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonwdev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array epsilonsdev_xx etc.'
+
+ allocate(epsilons_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ epsilonw_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array epsilons_trace_over_3 etc.'
+
+
read(27) rmass_solid_poroelastic
read(27) rmass_fluid_poroelastic
read(27) rhoarraystore
@@ -802,4 +820,147 @@
if( ier /= 0 ) stop 'error allocating array dsdx_top etc.'
endif
+ ! allocates adjoint arrays for poroelastic simulations
+ if( POROELASTIC_SIMULATION .and. SIMULATION_TYPE == 3 ) then
+ ! backward displacement,velocity,acceleration for the solid (s) & fluid (w) phases
+ allocate(b_displs_poroelastic(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_displs_poroelastic'
+ allocate(b_velocs_poroelastic(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_velocs_poroelastic'
+ allocate(b_accels_poroelastic(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_accels_poroelastic'
+ allocate(b_displw_poroelastic(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_displw_poroelastic'
+ allocate(b_velocw_poroelastic(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_velocw_poroelastic'
+ allocate(b_accelw_poroelastic(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_accelw_poroelastic'
+
+ ! adjoint kernels
+
+ ! primary, isotropic kernels
+ allocate(rhot_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ rhof_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ sm_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ eta_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rhot_kl etc.'
+ allocate(mufr_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mufr_kl'
+ allocate(B_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ C_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ M_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), stat=ier)
+ if( ier /= 0 ) stop 'error allocating array B_kl etc.'
+
+ ! density, isotropic kernels
+ allocate(rhob_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ rhofb_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ phi_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rhob_kl etc.'
+ allocate(mufrb_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mufrb_kl'
+ allocate(Bb_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ Cb_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ Mb_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), stat=ier)
+ if( ier /= 0 ) stop 'error allocating array Bb_kl etc.'
+
+ ! wavespeed, isotropic kernels
+ allocate(rhobb_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ rhofbb_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ phib_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ ratio_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rhobb_kl etc.'
+ allocate(cs_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array cs_kl'
+ allocate(cpI_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ cpII_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), stat=ier)
+ if( ier /= 0 ) stop 'error allocating array cpI_kl etc.'
+
+ ! MPI handling
+ allocate(b_request_send_vector_ext_meshs(num_interfaces_ext_mesh), &
+ b_request_recv_vector_ext_meshs(num_interfaces_ext_mesh), &
+ b_buffer_send_vector_ext_meshs(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ b_buffer_recv_vector_ext_meshs(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_request_send_vector_ext_meshs etc.'
+
+ allocate(b_request_send_vector_ext_meshw(num_interfaces_ext_mesh), &
+ b_request_recv_vector_ext_meshw(num_interfaces_ext_mesh), &
+ b_buffer_send_vector_ext_meshw(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ b_buffer_recv_vector_ext_meshw(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_request_send_vector_ext_meshw etc.'
+
+ ! arrays needed for kernel computations
+ allocate(b_epsilonsdev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonsdev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonsdev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonsdev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonsdev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonwdev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonwdev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonwdev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonwdev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonwdev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_epsilonsdev_xx etc.'
+
+ allocate(b_epsilons_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ b_epsilonw_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_epsilons_trace_over_3 etc.'
+
+ else ! dummy arrays
+
+ ! backward displacement,velocity,acceleration for the solid (s) & fluid (w)
+ ! phases
+ allocate(b_displs_poroelastic(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_displs_poroelastic'
+ allocate(b_velocs_poroelastic(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_velocs_poroelastic'
+ allocate(b_accels_poroelastic(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_accels_poroelastic'
+ allocate(b_displw_poroelastic(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_displw_poroelastic'
+ allocate(b_velocw_poroelastic(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_velocw_poroelastic'
+ allocate(b_accelw_poroelastic(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_accelw_poroelastic'
+
+ ! adjoint kernels
+
+ ! primary, isotropic kernels
+ allocate(rhot_kl(1,1,1,1), &
+ rhof_kl(1,1,1,1), &
+ sm_kl(1,1,1,1), &
+ eta_kl(1,1,1,1), stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array rhot_kl etc.'
+ allocate(mufr_kl(1,1,1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array mufr_kl'
+ allocate(B_kl(1,1,1,1), &
+ C_kl(1,1,1,1), &
+ M_kl(1,1,1,1), stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array B_kl etc.'
+
+ ! density, isotropic kernels
+ allocate(rhob_kl(1,1,1,1), &
+ rhofb_kl(1,1,1,1), &
+ phi_kl(1,1,1,1), stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array rhob_kl etc.'
+ allocate(mufrb_kl(1,1,1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array mufrb_kl'
+ allocate(Bb_kl(1,1,1,1), &
+ Cb_kl(1,1,1,1), &
+ Mb_kl(1,1,1,1), stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array Bb_kl etc.'
+
+ ! wavespeed, isotropic kernels
+ allocate(rhobb_kl(1,1,1,1), &
+ rhofbb_kl(1,1,1,1), &
+ phib_kl(1,1,1,1), &
+ ratio_kl(1,1,1,1), stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array rhobb_kl etc.'
+ allocate(cs_kl(1,1,1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array cs_kl'
+ allocate(cpI_kl(1,1,1,1), &
+ cpII_kl(1,1,1,1), stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array cpI_kl etc.'
+
+ endif
+
end subroutine read_mesh_databases_adjoint
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -31,11 +31,20 @@
use specfem_par
use specfem_par_elastic
use specfem_par_acoustic
+ use specfem_par_poroelastic
implicit none
integer:: ispec,i,j,k,iglob,ier
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) :: mul_s,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
+ real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,B_biot,cpIsquare,cpIIsquare,cssquare
+ real(kind=CUSTOM_REAL) :: rholb,ratio,dd1,gamma1,gamma2,gamma3,gamma4
+ real(kind=CUSTOM_REAL) :: afactor,bfactor,cfactor
! flag to save GLL weights
logical,parameter :: SAVE_WEIGHTS = .false.
@@ -108,7 +117,213 @@
endif ! acoustic
+ ! poroelastic simulations
+ if( ispec_is_poroelastic(ispec) ) then
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ ! isotropic adjoint kernels (see e.g. Morency et al. 2009)
+
+! get poroelastic parameters of current local GLL
+ 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
+ kappal_s = kappaarraystore(1,i,j,k,ispec)
+ kappal_f = kappaarraystore(2,i,j,k,ispec)
+ kappal_fr = kappaarraystore(3,i,j,k,ispec)
+ mul_fr = mustore(i,j,k,ispec)
+ etal_f = etastore(i,j,k,ispec)
+ permlxx = permstore(1,i,j,k,ispec)
+ permlxy = permstore(2,i,j,k,ispec)
+ permlxz = permstore(3,i,j,k,ispec)
+ permlyy = permstore(4,i,j,k,ispec)
+ permlyz = permstore(5,i,j,k,ispec)
+ permlzz = permstore(6,i,j,k,ispec)
+! Biot coef
+ 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
+ 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)
+ afactor = rhol_bar - phil/tortl*rhol_f
+ bfactor = H_biot + phil*rhol_bar/(tortl*rhol_f)*M_biot - 2._CUSTOM_REAL*phil/tortl*C_biot
+ cfactor = phil/(tortl*rhol_f)*(H_biot*M_biot - C_biot*C_biot)
+ cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(2._CUSTOM_REAL*afactor)
+ cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4._CUSTOM_REAL*afactor*cfactor))/(2._CUSTOM_REAL*afactor)
+ cssquare = mul_fr/afactor
+! extras needed
+! Approximated ratio r = amplitude "w" field/amplitude "s" field (no viscous
+! dissipation)
+ gamma1 = H_biot - phil/tortl*C_biot
+ gamma2 = C_biot - phil/tortl*M_biot
+ gamma3 = phil/tortl*( M_biot*(afactor/rhol_f + phil/tortl) - C_biot)
+ gamma4 = phil/tortl*( C_biot*(afactor/rhol_f + phil/tortl) - H_biot)
+ ratio = 1._CUSTOM_REAL/2._CUSTOM_REAL*(gamma1 - gamma3)/gamma4 + &
+ 1._CUSTOM_REAL/2._CUSTOM_REAL*sqrt((gamma1-gamma3)**2/gamma4**2 + 4._CUSTOM_REAL * gamma2/gamma4)
+ rholb = rhol_bar - phil*rhol_f/tortl
+ dd1 = (1._CUSTOM_REAL+rholb/rhol_f)*ratio**2 + 2._CUSTOM_REAL*ratio + tortl/phil
+
+ ! primary kernels
+ rhot_kl(i,j,k,ispec) = - rhol_bar * rhot_kl(i,j,k,ispec)
+ rhof_kl(i,j,k,ispec) = - rhol_f * rhof_kl(i,j,k,ispec)
+ sm_kl(i,j,k,ispec) = - rhol_f*tortl/phil * sm_kl(i,j,k,ispec)
+ !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)
+
+ ! 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)
+ rhofb_kl(i,j,k,ispec) = rhof_kl(i,j,k,ispec) + C_kl(i,j,k,ispec) + M_kl(i,j,k,ispec) + sm_kl(i,j,k,ispec)
+ Bb_kl(i,j,k,ispec) = B_kl(i,j,k,ispec)
+ Cb_kl(i,j,k,ispec) = C_kl(i,j,k,ispec)
+ Mb_kl(i,j,k,ispec) = M_kl(i,j,k,ispec)
+ mufrb_kl(i,j,k,ispec) = mufr_kl(i,j,k,ispec)
+ phi_kl(i,j,k,ispec) = - sm_kl(i,j,k,ispec) - M_kl(i,j,k,ispec)
+
+ ! wavespeed kernels
+ rhobb_kl(i,j,k,ispec) = rhob_kl(i,j,k,ispec) - &
+ phil*rhol_f/(tortl*B_biot) * &
+ (cpIIsquare + (cpIsquare - cpIIsquare)*( (phil / &
+ tortl*ratio +1._CUSTOM_REAL)/dd1 + &
+ (rhol_bar**2*ratio**2/rhol_f**2*(phil / &
+ tortl*ratio+1)*(phil/tortl*ratio + &
+ phil/tortl * &
+ (1+rhol_f/rhol_bar)-1))/dd1**2 ) - &
+ 4._CUSTOM_REAL/3._CUSTOM_REAL*cssquare )*Bb_kl(i,j,k,ispec) - &
+ rhol_bar*ratio**2/M_biot * (cpIsquare - cpIIsquare)* &
+ (phil/tortl*ratio + &
+ 1._CUSTOM_REAL)**2/dd1**2*Mb_kl(i,j,k,ispec) + &
+ rhol_bar*ratio/C_biot * (cpIsquare - cpIIsquare)* (&
+ (phil/tortl*ratio+1._CUSTOM_REAL)/dd1 - &
+ phil*ratio/tortl*(phil / &
+ tortl*ratio+1._CUSTOM_REAL)*&
+ (1+rhol_bar*ratio/rhol_f)/dd1**2)*Cb_kl(i,j,k,ispec)+ &
+ phil*rhol_f*cssquare / &
+ (tortl*mul_fr)*mufrb_kl(i,j,k,ispec)
+ rhofbb_kl(i,j,k,ispec) = rhofb_kl(i,j,k,ispec) + &
+ phil*rhol_f/(tortl*B_biot) * &
+ (cpIIsquare + (cpIsquare - cpIIsquare)*( (phil/ &
+ tortl*ratio +1._CUSTOM_REAL)/dd1+&
+ (rhol_bar**2*ratio**2/rhol_f**2*(phil/ &
+ tortl*ratio+1)*(phil/tortl*ratio+ &
+ phil/tortl*&
+ (1+rhol_f/rhol_bar)-1))/dd1**2 )- &
+ 4._CUSTOM_REAL/3._CUSTOM_REAL*cssquare )*Bb_kl(i,j,k,ispec) + &
+ rhol_bar*ratio**2/M_biot * (cpIsquare - cpIIsquare)* &
+ (phil/tortl*ratio + &
+ 1._CUSTOM_REAL)**2/dd1**2*Mb_kl(i,j,k,ispec) - &
+ rhol_bar*ratio/C_biot * (cpIsquare - cpIIsquare)* (&
+ (phil/tortl*ratio+1._CUSTOM_REAL)/dd1 - &
+ phil*ratio/tortl*(phil/ &
+ tortl*ratio+1._CUSTOM_REAL)*&
+ (1+rhol_bar*ratio/rhol_f)/dd1**2)*Cb_kl(i,j,k,ispec)- &
+ phil*rhol_f*cssquare/ &
+ (tortl*mul_fr)*mufrb_kl(i,j,k,ispec)
+ phib_kl(i,j,k,ispec) = phi_kl(i,j,k,ispec) - &
+ phil*rhol_bar/(tortl*B_biot) &
+ * ( cpIsquare - rhol_f/rhol_bar*cpIIsquare- &
+ (cpIsquare-cpIIsquare)*( (2._CUSTOM_REAL*ratio**2*phil/ &
+ tortl + (1._CUSTOM_REAL+&
+ rhol_f/rhol_bar)* &
+ (2._CUSTOM_REAL*ratio*phil/tortl+&
+ 1._CUSTOM_REAL))/dd1 + (phil/tortl*ratio+ &
+ 1._CUSTOM_REAL)*(phil*&
+ ratio/tortl+phil/tortl* &
+ (1._CUSTOM_REAL+rhol_f/&
+ rhol_bar)-1._CUSTOM_REAL)*((1._CUSTOM_REAL+ &
+ rhol_bar/rhol_f-&
+ 2._CUSTOM_REAL*phil/tortl)*ratio**2+2._CUSTOM_REAL*ratio)/dd1**2) - &
+ 4._CUSTOM_REAL/3._CUSTOM_REAL*rhol_f*cssquare/rhol_bar)*Bb_kl(i,j,k,ispec) + &
+ rhol_f/M_biot * (cpIsquare-cpIIsquare)*(&
+ 2._CUSTOM_REAL*ratio*(phil/tortl*ratio+1._CUSTOM_REAL)/dd1 - &
+ (phil/tortl*ratio+1._CUSTOM_REAL)**2*( &
+ (1._CUSTOM_REAL+rhol_bar/&
+ rhol_f-2._CUSTOM_REAL*phil/tortl)*ratio**2+2._CUSTOM_REAL*ratio)/dd1**2 &
+ )*Mb_kl(i,j,k,ispec) + &
+ phil*rhol_f/(tortl*C_biot)* &
+ (cpIsquare-cpIIsquare)*ratio* (&
+ (1._CUSTOM_REAL+rhol_f/rhol_bar*ratio)/dd1 - &
+ (phil/tortl*ratio+1._CUSTOM_REAL)* &
+ (1._CUSTOM_REAL+rhol_bar/&
+ rhol_f*ratio)*((1._CUSTOM_REAL+rhol_bar/rhol_f-2._CUSTOM_REAL*&
+ phil/tortl)*ratio+2._CUSTOM_REAL)/dd1**2&
+ )*Cb_kl(i,j,k,ispec) -&
+ phil*rhol_f*cssquare &
+ /(tortl*mul_fr)*mufrb_kl(i,j,k,ispec)
+ cpI_kl(i,j,k,ispec) = 2._CUSTOM_REAL*cpIsquare/B_biot*rhol_bar*( &
+ 1._CUSTOM_REAL-phil/tortl + &
+ (phil/tortl*ratio+ &
+ 1._CUSTOM_REAL)*(phil/tortl*&
+ ratio+phil/tortl* &
+ (1._CUSTOM_REAL+rhol_f/rhol_bar)-&
+ 1._CUSTOM_REAL)/dd1 &
+ )* Bb_kl(i,j,k,ispec) +&
+ 2._CUSTOM_REAL*cpIsquare*rhol_f*tortl/(phil*M_biot) *&
+ (phil/tortl*ratio+1._CUSTOM_REAL)**2/dd1*Mb_kl(i,j,k,ispec)+&
+ 2._CUSTOM_REAL*cpIsquare*rhol_f/C_biot * &
+ (phil/tortl*ratio+1._CUSTOM_REAL)* &
+ (1._CUSTOM_REAL+rhol_bar/&
+ rhol_f*ratio)/dd1*Cb_kl(i,j,k,ispec)
+ cpII_kl(i,j,k,ispec) = 2._CUSTOM_REAL*cpIIsquare*rhol_bar/B_biot * (&
+ phil*rhol_f/(tortl*rhol_bar) - &
+ (phil/tortl*ratio+ &
+ 1._CUSTOM_REAL)*(phil/tortl*&
+ ratio+phil/tortl* &
+ (1._CUSTOM_REAL+rhol_f/rhol_bar)-&
+ 1._CUSTOM_REAL)/dd1 ) * Bb_kl(i,j,k,ispec) +&
+ 2._CUSTOM_REAL*cpIIsquare*rhol_f*tortl/(phil*M_biot) * (&
+ 1._CUSTOM_REAL - (phil/tortl*ratio+ &
+ 1._CUSTOM_REAL)**2/dd1 )*Mb_kl(i,j,k,ispec) + &
+ 2._CUSTOM_REAL*cpIIsquare*rhol_f/C_biot * (&
+ 1._CUSTOM_REAL - (phil/tortl*ratio+ &
+ 1._CUSTOM_REAL)*(1._CUSTOM_REAL+&
+ rhol_bar/rhol_f*ratio)/dd1)*Cb_kl(i,j,k,ispec)
+ cs_kl(i,j,k,ispec) = - 8._CUSTOM_REAL/3._CUSTOM_REAL*cssquare* &
+ rhol_bar/B_biot*(1._CUSTOM_REAL-&
+ phil*rhol_f/(tortl* &
+ rhol_bar))*Bb_kl(i,j,k,ispec) + &
+ 2._CUSTOM_REAL*(rhol_bar-rhol_f*&
+ phil/tortl)/&
+ mul_fr*cssquare*mufrb_kl(i,j,k,ispec)
+ ratio_kl(i,j,k,ispec) = ratio*rhol_bar*phil/(tortl*B_biot) * &
+ (cpIsquare-cpIIsquare) * ( &
+ phil/tortl*(2._CUSTOM_REAL*ratio+1._CUSTOM_REAL+rhol_f/ &
+ rhol_bar)/dd1 - (phil/tortl*ratio+1._CUSTOM_REAL)*&
+ (phil/tortl*ratio+phil/tortl*(&
+ 1._CUSTOM_REAL+rhol_f/rhol_bar)-1._CUSTOM_REAL)*(2._CUSTOM_REAL*ratio*(&
+ 1._CUSTOM_REAL+rhol_bar/rhol_f-phil/tortl) +&
+ 2._CUSTOM_REAL)/dd1**2 )*Bb_kl(i,j,k,ispec) + &
+ ratio*rhol_f*tortl/(phil*M_biot)*(cpIsquare-cpIIsquare) * &
+ 2._CUSTOM_REAL*phil/tortl * (&
+ (phil/tortl*ratio+1._CUSTOM_REAL)/dd1 - &
+ (phil/tortl*ratio+1._CUSTOM_REAL)**2*( &
+ (1._CUSTOM_REAL+rhol_bar/&
+ rhol_f-phil/tortl)*ratio+ &
+ 1._CUSTOM_REAL)/dd1**2 )*Mb_kl(i,j,k,ispec) +&
+ ratio*rhol_f/C_biot*(cpIsquare-cpIIsquare) * (&
+ (2._CUSTOM_REAL*phil*rhol_bar* &
+ ratio/(tortl*rhol_f)+&
+ phil/tortl+rhol_bar/rhol_f)/dd1 - &
+ 2._CUSTOM_REAL*phil/tortl*(phil/tortl*ratio+&
+ 1._CUSTOM_REAL)*(1._CUSTOM_REAL+rhol_bar/rhol_f*ratio)*((1._CUSTOM_REAL+&
+ rhol_bar/rhol_f- &
+ phil/tortl)*ratio+1._CUSTOM_REAL)/&
+ dd1**2 )*Cb_kl(i,j,k,ispec)
+
+ enddo
+ enddo
+ enddo
+
+ endif ! poroelastic
+
enddo
! save kernels to binary files
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-07-15 13:59:46 UTC (rev 20524)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-07-15 15:31:06 UTC (rev 20525)
@@ -459,6 +459,12 @@
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accels_poroelastic,velocs_poroelastic,displs_poroelastic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accelw_poroelastic,velocw_poroelastic,displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ epsilonsdev_xx,epsilonsdev_yy,epsilonsdev_xy,epsilonsdev_xz,epsilonsdev_yz, &
+ epsilonwdev_xx,epsilonwdev_yy,epsilonwdev_xy,epsilonwdev_xz,epsilonwdev_yz
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ epsilons_trace_over_3,epsilonw_trace_over_3
+
! material properties
! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: mustore
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: etastore,tortstore
@@ -488,8 +494,14 @@
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_accels_poroelastic,b_velocs_poroelastic,b_displs_poroelastic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: b_accelw_poroelastic,b_velocw_poroelastic,b_displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ b_epsilonsdev_xx,b_epsilonsdev_yy,b_epsilonsdev_xy,b_epsilonsdev_xz,b_epsilonsdev_yz, &
+ b_epsilonwdev_xx,b_epsilonwdev_yy,b_epsilonwdev_xy,b_epsilonwdev_xz,b_epsilonwdev_yz
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ b_epsilons_trace_over_3,b_epsilonw_trace_over_3
+
! adjoint kernels [primary kernels, density kernels, wavespeed kernels]
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhot_kl, rhof_kl, sm_kl, eta_kl, mufr_kl, B_kl, &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhot_kl, rhof_kl, sm_kl, eta_kl, mufr_kl, B_kl, &
C_kl, M_kl, rhob_kl, rhofb_kl, phi_kl, Bb_kl, Cb_kl, Mb_kl, mufrb_kl, &
rhobb_kl, rhofbb_kl, phib_kl, cpI_kl, cpII_kl, cs_kl, ratio_kl
More information about the CIG-COMMITS
mailing list