[cig-commits] r19024 - seismo/3D/SPECFEM3D/trunk/src/specfem3D
cmorency at geodynamics.org
cmorency at geodynamics.org
Wed Oct 5 19:26:05 PDT 2011
Author: cmorency
Date: 2011-10-05 19:26:04 -0700 (Wed, 05 Oct 2011)
New Revision: 19024
Added:
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90
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
Log:
Merging of 3D Biot into SPECFEM3D with CUBIT support...
Added new files in specfem3D/ to handle poroelastic and acoustic-poro coupling.
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_poroelastic.f90 2011-10-06 02:26:04 UTC (rev 19024)
@@ -0,0 +1,247 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! 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
+
+ subroutine compute_add_sources_poroelastic( NSPEC_AB,NGLOB_AB, &
+ accels,accelw,&
+ rhoarraystore,phistore,tortstore,&
+ ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
+ xi_source,eta_source,gamma_source,nu_source, &
+ hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
+ ispec_is_poroelastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ nrec,islice_selected_rec,ispec_selected_rec, &
+ nadj_rec_local,adj_sourcearrays )
+
+ use specfem_par,only: PRINT_SOURCE_TIME_FUNCTION,stf_used_total
+ implicit none
+
+ include "constants.h"
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accels,accelw
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ phistore,tortstore
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! source
+ integer :: NSOURCES,myrank,it
+ integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
+ double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
+ double precision, dimension(3,3,NSOURCES) :: nu_source
+ double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,tshift_cmt
+ double precision :: dt,t0
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
+
+ double precision, external :: comp_source_time_function,comp_source_time_function_rickr,&
+ comp_source_time_function_gauss
+
+ logical, dimension(NSPEC_AB) :: ispec_is_poroelastic
+
+!adjoint simulations
+ integer:: SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT
+ integer:: nrec
+ integer,dimension(nrec) :: islice_selected_rec,ispec_selected_rec
+ integer:: nadj_rec_local
+ real(kind=CUSTOM_REAL),dimension(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ):: adj_sourcearrays
+
+! local parameters
+ double precision :: f0
+ double precision :: stf
+ real(kind=CUSTOM_REAL) stf_used,stf_used_total_all,time_source
+ integer :: isource,iglob,i,j,k,ispec
+! integer :: irec_local,irec
+ real(kind=CUSTOM_REAL) :: phil,tortl,rhol_s,rhol_f,rhol_bar
+
+! plotting source time function
+ if(PRINT_SOURCE_TIME_FUNCTION .and. .not. phase_is_inner ) then
+ ! initializes total
+ stf_used_total = 0.0_CUSTOM_REAL
+ endif
+
+! forward simulations
+ if (SIMULATION_TYPE == 1) then
+
+ do isource = 1,NSOURCES
+
+ ! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
+
+ ispec = ispec_selected_source(isource)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ if( ispec_is_poroelastic(ispec) ) then
+
+ if(USE_FORCE_POINT_SOURCE) then
+
+ ! note: for use_force_point_source xi/eta/gamma are in the range [1,NGLL*]
+ iglob = ibool(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+
+ ! get poroelastic parameters of current local GLL
+ phil = phistore(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ tortl = tortstore(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ rhol_s = rhoarraystore(1,nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ rhol_f = rhoarraystore(2,nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+ t0 = 1.2d0/f0
+
+ if (it == 1 .and. myrank == 0) then
+ print *,'using a source of dominant frequency ',f0
+!chris
+! print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+! print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ endif
+
+ ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
+ !stf_used = 1.d10 * comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+ stf_used = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_cmt(isource),f0)
+
+ ! we use nu_source(:,3) here because we want a source normal to the surface.
+! solid phase
+ accels(:,iglob) = accels(:,iglob) &
+ + (1._CUSTOM_REAL - phil/tortl) * sngl( nu_source(:,3,isource) ) * stf_used
+! fluid phase
+ accelw(:,iglob) = accelw(:,iglob) &
+ + (1._CUSTOM_REAL - rhol_f/rhol_bar) * sngl( nu_source(:,3,isource) ) * stf_used
+
+ else
+
+ !stf = comp_source_time_function(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+ !t0 = 1.2d0/hdur(isource)
+ !stf = comp_source_time_function_rickr(dble(it-1)*DT-t0-tshift_cmt(isource),hdur(isource))
+ stf = comp_source_time_function_gauss(dble(it-1)*DT-t0-tshift_cmt(isource),hdur_gaussian(isource))
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+
+ ! add source array
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ! 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
+! source in the solid phase only
+! solid phase
+ accels(:,iglob) = accels(:,iglob) &
+ + sourcearrays(isource,:,i,j,k)*stf_used
+! + (1._CUSTOM_REAL - phil/tortl)*sourcearrays(isource,:,i,j,k)*stf_used
+! fluid phase
+ accelw(:,iglob) = accelw(:,iglob) &
+ - rhol_f/rhol_bar*sourcearrays(isource,:,i,j,k)*stf_used
+! + (1._CUSTOM_REAL - rhol_f/rhol_bar)*sourcearrays(isource,:,i,j,k)*stf_used
+ enddo
+ enddo
+ enddo
+
+ endif ! USE_FORCE_POINT_SOURCE
+
+ stf_used_total = stf_used_total + stf_used
+
+ endif ! ispec_is_poroelastic
+ endif ! ispec_is_inner
+ endif ! myrank
+
+ enddo ! NSOURCES
+ endif ! forward
+
+! NOTE: adjoint sources and backward wavefield timing:
+! idea is to start with the backward field b_displ,.. at time (T)
+! and convolve with the adjoint field at time (T-t)
+!
+! backward/reconstructed wavefields:
+! time for b_displ( it ) corresponds to (NSTEP - it - 1 )*DT - t0 ...
+! since we start with saved wavefields b_displ( 0 ) = displ( NSTEP ) which correspond
+! to a time (NSTEP - 1)*DT - t0
+! (see sources for simulation_type 1 and seismograms)
+! now, at the beginning of the time loop, the numerical Newark time scheme updates
+! the wavefields, that is b_displ( it=1) corresponds now to time (NSTEP -1 - 1)*DT - t0
+!
+! let's define the start time t to (1-1)*DT - t0 = -t0, and the end time T to (NSTEP-1)*DT - t0
+! these are the start and end times of all seismograms
+!
+! adjoint wavefields:
+! since the adjoint source traces were derived from the seismograms,
+! it follows that for the adjoint wavefield, the time equivalent to ( T - t ) uses the time-reversed
+! adjoint source traces which start at -t0 and end at time (NSTEP-1)*DT - t0
+! for it=1: (NSTEP -1 - 1)*DT - t0 for backward wavefields corresponds to time T-1
+! and time (T-1) corresponds now to index (NSTEP -1) in the adjoint source array
+
+
+! adjoint simulations
+ if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
+ stop 'adjoint poroelastic simulation not implemented yet'
+ endif !adjoint
+
+! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ stop 'adjoint poroelastic simulation not implemented yet'
+ endif ! adjoint
+
+ ! master prints out source time function to file
+ if(PRINT_SOURCE_TIME_FUNCTION .and. phase_is_inner) then
+ time_source = (it-1)*DT - t0
+ call sum_all_cr(stf_used_total,stf_used_total_all)
+ if( myrank == 0 ) write(IOSTF,*) time_source,stf_used_total_all
+ endif
+
+
+ end subroutine compute_add_sources_poroelastic
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_acoustic_po.f90 2011-10-06 02:26:04 UTC (rev 19024)
@@ -0,0 +1,122 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! 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 acoustic solver
+
+ subroutine compute_coupling_acoustic_po(NSPEC_AB,NGLOB_AB, &
+ ibool,displs_poroelastic,displw_poroelastic, &
+ potential_dot_dot_acoustic, &
+ num_coupling_ac_po_faces, &
+ coupling_ac_po_ispec,coupling_ac_po_ijk, &
+ coupling_ac_po_normal, &
+ coupling_ac_po_jacobian2Dw, &
+ ispec_is_inner,phase_is_inner)
+
+! returns the updated pressure array: potential_dot_dot_acoustic
+
+ implicit none
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacement and pressure
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic,displw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic
+
+! global indexing
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+! acoustic-poroelastic coupling surface
+ integer :: num_coupling_ac_po_faces
+ real(kind=CUSTOM_REAL) :: coupling_ac_po_normal(NDIM,NGLLSQUARE,num_coupling_ac_po_faces)
+ real(kind=CUSTOM_REAL) :: coupling_ac_po_jacobian2Dw(NGLLSQUARE,num_coupling_ac_po_faces)
+ integer :: coupling_ac_po_ijk(3,NGLLSQUARE,num_coupling_ac_po_faces)
+ integer :: coupling_ac_po_ispec(num_coupling_ac_po_faces)
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! local parameters
+ real(kind=CUSTOM_REAL) :: displ_x,displ_y,displ_z,displ_n
+ real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw
+
+ integer :: iface,igll,ispec,iglob
+ integer :: i,j,k
+
+! loops on all coupling faces
+ do iface = 1,num_coupling_ac_po_faces
+
+ ! gets corresponding elements
+ ispec = coupling_ac_po_ispec(iface)
+
+ if( ispec_is_inner(ispec) .eqv. phase_is_inner ) then
+
+ ! loops over common GLL points
+ do igll = 1, NGLLSQUARE
+ i = coupling_ac_po_ijk(1,igll,iface)
+ j = coupling_ac_po_ijk(2,igll,iface)
+ k = coupling_ac_po_ijk(3,igll,iface)
+
+ ! gets global index of this common GLL point
+ ! (note: should be the same as for corresponding i',j',k',ispec_poroelastic or ispec_acoustic)
+ iglob = ibool(i,j,k,ispec)
+
+ ! poroelastic displacement on global point
+ displ_x = displs_poroelastic(1,iglob) + displw_poroelastic(1,iglob)
+ displ_y = displs_poroelastic(2,iglob) + displw_poroelastic(2,iglob)
+ displ_z = displs_poroelastic(3,iglob) + displw_poroelastic(3,iglob)
+
+ ! gets associated normal on GLL point
+ ! (note convention: pointing outwards of acoustic element)
+ nx = coupling_ac_po_normal(1,igll,iface)
+ ny = coupling_ac_po_normal(2,igll,iface)
+ nz = coupling_ac_po_normal(3,igll,iface)
+
+ ! calculates displacement component along normal
+ ! (normal points outwards of acoustic element)
+ displ_n = displ_x*nx + displ_y*ny + displ_z*nz
+
+ ! gets associated, weighted jacobian
+ jacobianw = coupling_ac_po_jacobian2Dw(igll,iface)
+
+ ! continuity of pressure and normal displacement on global point
+ !
+ ! note: newark time scheme together with definition of scalar potential:
+ ! pressure = - chi_dot_dot
+ ! requires that this coupling term uses the updated displacement at time step [t+delta_t],
+ ! which is done at the very beginning of the time loop
+ ! (see e.g. Chaljub & Vilotte, Nissen-Meyer thesis...)
+ ! it also means you have to calculate and update this here first before
+ ! calculating the coupling on the elastic side for the acceleration...
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + jacobianw*displ_n
+
+ enddo ! igll
+
+ endif
+
+ enddo ! iface
+
+end subroutine compute_coupling_acoustic_po
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_coupling_poroelastic_ac.f90 2011-10-06 02:26:04 UTC (rev 19024)
@@ -0,0 +1,146 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! 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
+
+ subroutine compute_coupling_poroelastic_ac(NSPEC_AB,NGLOB_AB, &
+ ibool,accels_poroelastic,accelw_poroelastic, &
+ potential_dot_dot_acoustic, &
+ num_coupling_ac_po_faces, &
+ coupling_ac_po_ispec,coupling_ac_po_ijk, &
+ coupling_ac_po_normal, &
+ coupling_ac_po_jacobian2Dw, &
+ rhoarraystore,phistore,tortstore, &
+ ispec_is_inner,phase_is_inner)
+
+! returns the updated accelerations array: accels_poroelatsic & accelw_poroelastic
+
+ implicit none
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacement and pressure
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: accels_poroelastic,accelw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: potential_dot_dot_acoustic
+
+! global indexing
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+! acoustic-poroelastic coupling surface
+ integer :: num_coupling_ac_po_faces
+ real(kind=CUSTOM_REAL) :: coupling_ac_po_normal(NDIM,NGLLSQUARE,num_coupling_ac_po_faces)
+ real(kind=CUSTOM_REAL) :: coupling_ac_po_jacobian2Dw(NGLLSQUARE,num_coupling_ac_po_faces)
+ integer :: coupling_ac_po_ijk(3,NGLLSQUARE,num_coupling_ac_po_faces)
+ integer :: coupling_ac_po_ispec(num_coupling_ac_po_faces)
+
+! properties
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ phistore,tortstore
+
+! communication overlap
+ logical, dimension(NSPEC_AB) :: ispec_is_inner
+ logical :: phase_is_inner
+
+! local parameters
+ real(kind=CUSTOM_REAL) :: pressure
+ real(kind=CUSTOM_REAL) :: rhol_s,rhol_f,phil,tortl,rhol_bar
+ real(kind=CUSTOM_REAL) :: nx,ny,nz,jacobianw
+
+ integer :: iface,igll,ispec,iglob
+ integer :: i,j,k
+
+! loops on all coupling faces
+ do iface = 1,num_coupling_ac_po_faces
+
+ ! gets corresponding spectral element, constructed such that it is a poroelastic
+ ! element, since we need to have access to porous properties
+ ispec = coupling_ac_po_ispec(iface)
+
+ if( ispec_is_inner(ispec) .eqv. phase_is_inner ) then
+
+ ! loops over common GLL points
+ do igll = 1, NGLLSQUARE
+ i = coupling_ac_po_ijk(1,igll,iface)
+ j = coupling_ac_po_ijk(2,igll,iface)
+ k = coupling_ac_po_ijk(3,igll,iface)
+
+ ! gets global index of this common GLL point
+ ! (note: should be the same as for corresponding i',j',k',ispec_poroelastic or ispec_acoustic )
+ iglob = ibool(i,j,k,ispec)
+
+ ! get poroelastic parameters
+ 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
+
+ ! acoustic pressure on global point
+ pressure = - potential_dot_dot_acoustic(iglob)
+
+ ! gets associated normal on GLL point
+ ! (note convention: pointing outwards of acoustic element)
+ nx = coupling_ac_po_normal(1,igll,iface)
+ ny = coupling_ac_po_normal(2,igll,iface)
+ nz = coupling_ac_po_normal(3,igll,iface)
+
+ ! gets associated, weighted 2D jacobian
+ ! (note: should be the same for poroelastic and acoustic element)
+ jacobianw = coupling_ac_po_jacobian2Dw(igll,iface)
+
+ ! continuity of displacement and pressure on global point
+ !
+ ! note: newark time scheme together with definition of scalar potential:
+ ! pressure = - chi_dot_dot
+ ! requires that this coupling term uses the *UPDATED* pressure (chi_dot_dot), i.e.
+ ! pressure at time step [t + delta_t]
+ ! (see e.g. Chaljub & Vilotte, Nissen-Meyer thesis...)
+ ! it means you have to calculate and update the acoustic pressure first before
+ ! calculating this term...
+! contribution to the solid phase
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + jacobianw*nx*pressure*&
+ (1._CUSTOM_REAL-phil/tortl)
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + jacobianw*ny*pressure*&
+ (1._CUSTOM_REAL-phil/tortl)
+ accels_poroelastic(3,iglob) = accels_poroelastic(3,iglob) + jacobianw*nz*pressure*&
+ (1._CUSTOM_REAL-phil/tortl)
+! contribution to the fluid phase
+ accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) + jacobianw*nx*pressure*&
+ (1._CUSTOM_REAL-rhol_f/rhol_bar)
+ accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + jacobianw*ny*pressure*&
+ (1._CUSTOM_REAL-rhol_f/rhol_bar)
+ accelw_poroelastic(3,iglob) = accelw_poroelastic(3,iglob) + jacobianw*nz*pressure*&
+ (1._CUSTOM_REAL-rhol_f/rhol_bar)
+
+ enddo ! igll
+
+ endif
+
+ enddo ! iface
+
+end subroutine compute_coupling_poroelastic_ac
+
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_fluid.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_fluid.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_fluid.f90 2011-10-06 02:26:04 UTC (rev 19024)
@@ -0,0 +1,578 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine compute_forces_fluid( iphase, &
+ NSPEC_AB,NGLOB_AB,displw_poroelastic,accelw_poroelastic,&
+ velocw_poroelastic,displs_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,ibool,&
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
+ phase_ispec_inner_poroelastic )
+
+! compute forces for the fluid poroelastic part
+
+ use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM, &
+ N_SLS, &
+ ONE_THIRD,FOUR_THIRDS
+
+ implicit none
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displw_poroelastic,accelw_poroelastic,&
+ velocw_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ mustore,etastore,phistore,tortstore,jacobian
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: kappaarraystore
+ real(kind=CUSTOM_REAL), dimension(6,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: permstore
+
+! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: wxgll
+ double precision, dimension(NGLLY) :: wygll
+ double precision, dimension(NGLLZ) :: wzgll
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+! logical,dimension(NSPEC_AB) :: ispec_is_elastic
+ integer :: iphase
+ integer :: num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic
+ integer, dimension(num_phase_ispec_poroelastic,2) :: phase_ispec_inner_poroelastic
+
+
+! adjoint simulations
+ integer :: SIMULATION_TYPE
+ integer :: NSPEC_BOUN
+ integer :: NGLOB_ADJOINT,NSPEC_ADJOINT
+! adjoint wavefields
+! real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT):: b_displs_poroelastic,b_accels_poroelastic
+! real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT):: b_displw_poroelastic
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: &
+! mufr_kl, B_kl
+
+
+!---
+!--- local variables
+!---
+
+ integer :: ispec,i,j,k,l,iglob,num_elements,ispec_p
+
+! spatial derivatives
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ tempx1p,tempx2p,tempx3p,tempy1p,tempy2p,tempy3p,tempz1p,tempz2p,tempz3p
+! b_tempx1,b_tempx2,b_tempx3,b_tempy1,b_tempy2,b_tempy3,b_tempz1,b_tempz2,b_tempz3, &
+! b_tempx1p,b_tempx2p,b_tempx3p,b_tempy1p,b_tempy2p,b_tempy3p,b_tempz1p,b_tempz2p,b_tempz3p
+
+ real(kind=CUSTOM_REAL) :: duxdxl,duydxl,duzdxl,duxdyl,duydyl,duzdyl,duxdzl,duydzl,duzdzl
+ real(kind=CUSTOM_REAL) :: dwxdxl,dwydxl,dwzdxl,dwxdyl,dwydyl,dwzdyl,dwxdzl,dwydzl,dwzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl_plus_duzdzl,dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
+
+ real(kind=CUSTOM_REAL) tempx1ls,tempx2ls,tempx3ls,tempx1lw,tempx2lw,tempx3lw
+ real(kind=CUSTOM_REAL) tempy1ls,tempy2ls,tempy3ls,tempy1lw,tempy2lw,tempy3lw
+ real(kind=CUSTOM_REAL) tempz1ls,tempz2ls,tempz3ls,tempz1lw,tempz2lw,tempz3lw
+
+! real(kind=CUSTOM_REAL) :: b_dux_dxl,b_duz_dxl,b_dux_dzl,b_duz_dzl
+! real(kind=CUSTOM_REAL) :: dsxx,dsxy,dsxz,dsyy,dsyz,dszz
+! real(kind=CUSTOM_REAL) :: b_dsxx,b_dsxy,b_dsxz,b_dsyy,b_dsyz,b_dszz
+! real(kind=CUSTOM_REAL) :: b_dwx_dxl,b_dwz_dxl,b_dwx_dzl,b_dwz_dzl
+! real(kind=CUSTOM_REAL) :: dwxx,dwxy,dwxz,dwyy,dwyz,dwzz
+! real(kind=CUSTOM_REAL) :: b_dwxx,b_dwxy,b_dwxz,b_dwyy,b_dwyz,b_dwzz
+ real(kind=CUSTOM_REAL) :: sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) :: sigmap
+! real(kind=CUSTOM_REAL) :: b_sigma_xx,b_sigma_yy,b_sigma_zz,b_sigma_xy,b_sigma_xz,b_sigma_yz
+! real(kind=CUSTOM_REAL) :: b_sigmap
+! real(kind=CUSTOM_REAL) :: nx,nz,vx,vz,vn,vxf,vzf,vnf,rho_vpI,rho_vpII,rho_vs,tx,tz,weight,xxi,zxi,xgamma,zgamma,jacobian1D
+
+! viscous attenuation (poroelastic media)
+ real(kind=CUSTOM_REAL), dimension(6) :: bl_relaxed
+
+
+! Jacobian matrix and determinant
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+
+! material properties of the poroelastic medium
+! real(kind=CUSTOM_REAL) :: mul_unrelaxed,lambdal_unrelaxed,lambdalplus2mul_unrelaxed
+ real(kind=CUSTOM_REAL) :: kappal_s,rhol_s
+ real(kind=CUSTOM_REAL) :: etal_f,kappal_f,rhol_f
+ real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr,phil,tortl,viscodampx,viscodampy,viscodampz
+ real(kind=CUSTOM_REAL) :: permlxx,permlxy,permlxz,permlyz,permlyy,permlzz,&
+ invpermlxx,invpermlxy,invpermlxz,invpermlyz,invpermlyy,invpermlzz,detk
+ real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,rhol_bar
+
+ real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
+! real(kind=CUSTOM_REAL) :: cpIsquare,cpIIsquare,cssquare,cpIl,cpIIl,csl
+
+! for attenuation
+! real(kind=CUSTOM_REAL) :: Un,Unp1,tauinv,Sn,Snp1,theta_n,theta_np1,tauinvsquare,tauinvcube,tauinvUn
+
+! compute Grad(displs_poroelastic) at time step n for attenuation
+! if(TURN_ATTENUATION_ON) call compute_gradient_attenuation(displs_poroelastic,dux_dxl_n,duz_dxl_n, &
+! dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,npoin)
+
+
+ if( iphase == 1 ) then
+ num_elements = nspec_outer_poroelastic
+ else
+ num_elements = nspec_inner_poroelastic
+ endif
+
+! loop over spectral elements
+ do ispec_p = 1,num_elements
+
+ ispec = phase_ispec_inner_poroelastic(ispec_p,iphase)
+
+! first double loop over GLL points to compute and store gradients
+ do k=1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+!solid properties
+ kappal_s = kappaarraystore(1,i,j,k,ispec)
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+!fluid properties
+ kappal_f = kappaarraystore(2,i,j,k,ispec)
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+!frame properties
+ mul_fr = mustore(i,j,k,ispec)
+ kappal_fr = kappaarraystore(3,i,j,k,ispec)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+!Biot coefficients for the input phi
+ 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)
+!The RHS has the form : div T -phi/c div T_f + phi/ceta_fk^-1.partial t w
+!where T = G:grad u_s + C_biot div w I
+!and T_f = C_biot div u_s I + M_biot div w I
+ mul_G = mul_fr
+ lambdal_G = H_biot - 2._CUSTOM_REAL*mul_fr
+ lambdalplus2mul_G = lambdal_G + 2._CUSTOM_REAL*mul_G
+
+! derivative along x,y,z for u_s and w
+ tempx1ls = 0.
+ tempx2ls = 0.
+ tempx3ls = 0.
+
+ tempy1ls = 0.
+ tempy2ls = 0.
+ tempy3ls = 0.
+
+ tempz1ls = 0.
+ tempz2ls = 0.
+ tempz3ls = 0.
+
+ tempx1lw = 0.
+ tempx2lw = 0.
+ tempx3lw = 0.
+
+ tempy1lw = 0.
+ tempy2lw = 0.
+ tempy3lw = 0.
+
+ tempz1lw = 0.
+ tempz2lw = 0.
+ tempz3lw = 0.
+
+ if (SIMULATION_TYPE == 3) then
+ endif
+! first double loop over GLL points to compute and store gradients
+ do l = 1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1ls = tempx1ls + displs_poroelastic(1,iglob)*hp1
+ tempy1ls = tempy1ls + displs_poroelastic(2,iglob)*hp1
+ tempz1ls = tempz1ls + displs_poroelastic(3,iglob)*hp1
+ tempx1lw = tempx1lw + displw_poroelastic(1,iglob)*hp1
+ tempy1lw = tempy1lw + displw_poroelastic(2,iglob)*hp1
+ tempz1lw = tempz1lw + displw_poroelastic(3,iglob)*hp1
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2ls = tempx2ls + displs_poroelastic(1,iglob)*hp2
+ tempy2ls = tempy2ls + displs_poroelastic(2,iglob)*hp2
+ tempz2ls = tempz2ls + displs_poroelastic(3,iglob)*hp2
+ tempx2lw = tempx2lw + displw_poroelastic(1,iglob)*hp2
+ tempy2lw = tempy2lw + displw_poroelastic(2,iglob)*hp2
+ tempz2lw = tempz2lw + displw_poroelastic(3,iglob)*hp2
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3ls = tempx3ls + displs_poroelastic(1,iglob)*hp3
+ tempy3ls = tempy3ls + displs_poroelastic(2,iglob)*hp3
+ tempz3ls = tempz3ls + displs_poroelastic(3,iglob)*hp3
+ tempx3lw = tempx3lw + displw_poroelastic(1,iglob)*hp3
+ tempy3lw = tempy3lw + displw_poroelastic(2,iglob)*hp3
+ tempz3lw = tempz3lw + displw_poroelastic(3,iglob)*hp3
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ enddo
+
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
+
+! derivatives of displacement
+ duxdxl = xixl*tempx1ls + etaxl*tempx2ls + gammaxl*tempx3ls
+ duxdyl = xiyl*tempx1ls + etayl*tempx2ls + gammayl*tempx3ls
+ duxdzl = xizl*tempx1ls + etazl*tempx2ls + gammazl*tempx3ls
+
+ duydxl = xixl*tempy1ls + etaxl*tempy2ls + gammaxl*tempy3ls
+ duydyl = xiyl*tempy1ls + etayl*tempy2ls + gammayl*tempy3ls
+ duydzl = xizl*tempy1ls + etazl*tempy2ls + gammazl*tempy3ls
+
+ duzdxl = xixl*tempz1ls + etaxl*tempz2ls + gammaxl*tempz3ls
+ duzdyl = xiyl*tempz1ls + etayl*tempz2ls + gammayl*tempz3ls
+ duzdzl = xizl*tempz1ls + etazl*tempz2ls + gammazl*tempz3ls
+
+ dwxdxl = xixl*tempx1lw + etaxl*tempx2lw + gammaxl*tempx3lw
+ dwxdyl = xiyl*tempx1lw + etayl*tempx2lw + gammayl*tempx3lw
+ dwxdzl = xizl*tempx1lw + etazl*tempx2lw + gammazl*tempx3lw
+
+ dwydxl = xixl*tempy1lw + etaxl*tempy2lw + gammaxl*tempy3lw
+ dwydyl = xiyl*tempy1lw + etayl*tempy2lw + gammayl*tempy3lw
+ dwydzl = xizl*tempy1lw + etazl*tempy2lw + gammazl*tempy3lw
+
+ dwzdxl = xixl*tempz1lw + etaxl*tempz2lw + gammaxl*tempz3lw
+ dwzdyl = xiyl*tempz1lw + etayl*tempz2lw + gammayl*tempz3lw
+ dwzdzl = xizl*tempz1lw + etazl*tempz2lw + gammazl*tempz3lw
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl_plus_duzdzl = duxdxl + duydyl + duzdzl
+ dwxdxl_plus_dwydyl_plus_dwzdzl = dwxdxl + dwydyl + dwzdzl
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! compute stress tensor (include attenuation or anisotropy if needed)
+
+! if(VISCOATTENUATION) then
+!chris:check
+
+! Dissipation only controlled by frame share attenuation in poroelastic (see Morency & Tromp, GJI 2008).
+! attenuation is implemented following the memory variable formulation of
+! J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
+! vol. 58(1), p. 110-120 (1993). More details can be found in
+! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a linear
+! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611 (1988).
+
+! else
+
+! no attenuation
+ sigma_xx = lambdalplus2mul_G*duxdxl + lambdal_G*duydyl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_yy = lambdalplus2mul_G*duydyl + lambdal_G*duxdxl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_zz = lambdalplus2mul_G*duzdzl + lambdal_G*duxdxl_plus_duydyl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ sigma_xy = mul_G*duxdyl_plus_duydxl
+ sigma_xz = mul_G*duzdxl_plus_duxdzl
+ sigma_yz = mul_G*duzdyl_plus_duydzl
+
+ sigmap = C_biot*duxdxl_plus_duydyl_plus_duzdzl + M_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ if(SIMULATION_TYPE == 3) then ! kernels calculation
+ endif
+! endif !if(VISCOATTENUATION)
+
+! weak formulation term based on stress tensor (non-symmetric form)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ tempx1p(i,j,k) = jacobianl * sigmap*xixl
+ tempy1p(i,j,k) = jacobianl * sigmap*xiyl
+ tempz1p(i,j,k) = jacobianl * sigmap*xizl
+
+ tempx2p(i,j,k) = jacobianl * sigmap*etaxl
+ tempy2p(i,j,k) = jacobianl * sigmap*etayl
+ tempz2p(i,j,k) = jacobianl * sigmap*etazl
+
+ tempx3p(i,j,k) = jacobianl * sigmap*gammaxl
+ tempy3p(i,j,k) = jacobianl * sigmap*gammayl
+ tempz3p(i,j,k) = jacobianl * sigmap*gammazl
+
+ if(SIMULATION_TYPE == 3) then ! adjoint calculation
+ endif
+
+ enddo
+ enddo
+ enddo
+
+!
+! second double-loop over GLL to compute all the terms
+!
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ tempx1ls = 0.
+ tempy1ls = 0.
+ tempz1ls = 0.
+
+ tempx2ls = 0.
+ tempy2ls = 0.
+ tempz2ls = 0.
+
+ tempx3ls = 0.
+ tempy3ls = 0.
+ tempz3ls = 0.
+
+ tempx1lw = 0.
+ tempy1lw = 0.
+ tempz1lw = 0.
+
+ tempx2lw = 0.
+ tempy2lw = 0.
+ tempz2lw = 0.
+
+ tempx3lw = 0.
+ tempy3lw = 0.
+ tempz3lw = 0.
+
+ if (SIMULATION_TYPE == 3) then
+ endif
+
+ do l=1,NGLLX
+ fac1 = hprimewgll_xx(l,i)
+ tempx1ls = tempx1ls + tempx1(l,j,k)*fac1
+ tempy1ls = tempy1ls + tempy1(l,j,k)*fac1
+ tempz1ls = tempz1ls + tempz1(l,j,k)*fac1
+ tempx1lw = tempx1lw + tempx1p(l,j,k)*fac1
+ tempy1lw = tempy1lw + tempy1p(l,j,k)*fac1
+ tempz1lw = tempz1lw + tempz1p(l,j,k)*fac1
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ fac2 = hprimewgll_yy(l,j)
+ tempx2ls = tempx2ls + tempx2(i,l,k)*fac2
+ tempy2ls = tempy2ls + tempy2(i,l,k)*fac2
+ tempz2ls = tempz2ls + tempz2(i,l,k)*fac2
+ tempx2lw = tempx2lw + tempx2p(i,l,k)*fac2
+ tempy2lw = tempy2lw + tempy2p(i,l,k)*fac2
+ tempz2lw = tempz2lw + tempz2p(i,l,k)*fac2
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ fac3 = hprimewgll_zz(l,k)
+ tempx3ls = tempx3ls + tempx3(i,j,l)*fac3
+ tempy3ls = tempy3ls + tempy3(i,j,l)*fac3
+ tempz3ls = tempz3ls + tempz3(i,j,l)*fac3
+ tempx3lw = tempx3lw + tempx3p(i,j,l)*fac3
+ tempy3lw = tempy3lw + tempy3p(i,j,l)*fac3
+ tempz3lw = tempz3lw + tempz3p(i,j,l)*fac3
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif
+ enddo
+
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec)
+!solid properties
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+!fluid properties
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+!frame properties
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+
+ ! sum contributions from each element to the global mesh
+
+ iglob = ibool(i,j,k,ispec)
+
+
+ accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) + ( fac1*(rhol_f/rhol_bar*tempx1ls - tempx1lw) &
+ + fac2*(rhol_f/rhol_bar*tempx2ls - tempx2lw) + fac3*(rhol_f/rhol_bar*tempx3ls - tempx3lw) )
+
+ accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) + ( fac1*(rhol_f/rhol_bar*tempy1ls - tempy1lw) &
+ + fac2*(rhol_f/rhol_bar*tempy2ls - tempy2lw) + fac3*(rhol_f/rhol_bar*tempy3ls - tempy3lw) )
+
+ 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
+!
+! add + phi/tort eta_f k^-1 dot(w)
+
+ etal_f = etastore(i,j,k,ispec)
+
+ if(etal_f >0.d0) then
+
+ 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)
+
+! calcul of the inverse of k
+ detk = permlxz*(permlxy*permlyz-permlxz*permlyy) &
+ - permlxy*(permlxy*permlzz-permlyz*permlxz) &
+ + permlxx*(permlyy*permlzz-permlyz*permlyz)
+
+ if(detk /= 0.d0) then
+ invpermlxx = (permlyy*permlzz-permlyz*permlyz)/detk
+ invpermlxy = (permlxz*permlyz-permlxy*permlzz)/detk
+ invpermlxz = (permlxy*permlyz-permlxz*permlyy)/detk
+ invpermlyy = (permlxx*permlzz-permlxz*permlxz)/detk
+ invpermlyz = (permlxy*permlxz-permlxx*permlyz)/detk
+ invpermlzz = (permlxx*permlyy-permlxy*permlxy)/detk
+ else
+ stop 'Permeability matrix is not inversible'
+ endif
+
+! relaxed viscous coef
+ bl_relaxed(1) = etal_f*invpermlxx
+ bl_relaxed(2) = etal_f*invpermlxy
+ bl_relaxed(3) = etal_f*invpermlxz
+ bl_relaxed(4) = etal_f*invpermlyy
+ bl_relaxed(5) = etal_f*invpermlyz
+ bl_relaxed(6) = etal_f*invpermlzz
+
+! 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
+! endif
+
+! do k = 1,NGLLZ
+! do j = 1,NGLLY
+! do i = 1,NGLLX
+
+! iglob = ibool(i,j,k,ispec)
+
+! if(VISCOATTENUATION) then
+! compute the viscous damping term with the unrelaxed viscous coef and add memory variable
+! viscodampx = velocw_poroelastic(1,iglob)*bl_unrelaxed(1) + velocw_poroelastic(2,iglob)*bl_unrelaxed(2)&
+! - rx_viscous(i,j,ispec)
+! viscodampz = velocw_poroelastic(1,iglob)*bl_unrelaxed(2) + velocw_poroelastic(2,iglob)*bl_unrelaxed(3)&
+! - rz_viscous(i,j,ispec)
+! else
+
+! no viscous attenuation
+ viscodampx = velocw_poroelastic(1,iglob)*bl_relaxed(1) + velocw_poroelastic(2,iglob)*bl_relaxed(2) + &
+ velocw_poroelastic(3,iglob)*bl_relaxed(3)
+ viscodampy = velocw_poroelastic(1,iglob)*bl_relaxed(2) + velocw_poroelastic(2,iglob)*bl_relaxed(4) + &
+ velocw_poroelastic(3,iglob)*bl_relaxed(5)
+ viscodampz = velocw_poroelastic(1,iglob)*bl_relaxed(3) + velocw_poroelastic(2,iglob)*bl_relaxed(5) + &
+ velocw_poroelastic(3,iglob)*bl_relaxed(6)
+! endif
+
+ accelw_poroelastic(1,iglob) = accelw_poroelastic(1,iglob) - wxgll(i)*wygll(j)*wzgll(k)*jacobian(i,j,k,ispec)*&
+ viscodampx
+ accelw_poroelastic(2,iglob) = accelw_poroelastic(2,iglob) - wxgll(i)*wygll(j)*wzgll(k)*jacobian(i,j,k,ispec)*&
+ viscodampy
+ accelw_poroelastic(3,iglob) = accelw_poroelastic(3,iglob) - wxgll(i)*wygll(j)*wzgll(k)*jacobian(i,j,k,ispec)*&
+ viscodampz
+
+! if isolver == 1 .and. save_forward then b_viscodamp is save in compute_forces_fluid.f90
+! if(isolver == 2) then ! kernels calculation
+! b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + phil/tortl*b_viscodampx(iglob)
+! b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + phil/tortl*b_viscodampz(iglob)
+! endif
+
+! enddo
+! enddo
+! enddo
+
+ endif ! if(etal_f >0.d0) then
+
+ enddo ! second loop over the GLL points
+ enddo
+ enddo
+
+ enddo ! end of loop over all spectral elements
+
+
+ end subroutine compute_forces_fluid
+
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_poroelastic.f90 2011-10-06 02:26:04 UTC (rev 19024)
@@ -0,0 +1,260 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! 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.
+!
+!=====================================================================
+
+! poroelastic solver
+
+subroutine compute_forces_poroelastic()
+
+ use specfem_par
+ use specfem_par_acoustic
+ use specfem_par_elastic
+ use specfem_par_poroelastic
+
+ implicit none
+
+ integer:: iphase
+ logical:: phase_is_inner
+
+! distinguishes two runs: for points on MPI interfaces, and points within the partitions
+ do iphase=1,2
+
+ !first for points on MPI interfaces
+ if( iphase == 1 ) then
+ phase_is_inner = .false.
+ else
+ phase_is_inner = .true.
+ endif
+
+!Note: Contrary to the elastic & acoustic case, the absorbing implementation is within compute_forces
+!due to the number of properties which were needed
+
+! solid phase
+
+ call compute_forces_solid( iphase, &
+ NSPEC_AB,NGLOB_AB,displs_poroelastic,accels_poroelastic,&
+ velocs_poroelastic,displw_poroelastic,velocw_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,ibool,&
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
+ phase_ispec_inner_poroelastic )
+
+! fluid phase
+
+ call compute_forces_fluid( iphase, &
+ NSPEC_AB,NGLOB_AB,displw_poroelastic,accelw_poroelastic,&
+ velocw_poroelastic,displs_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,ibool,&
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ 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 implemented yet'
+ endif
+
+
+! acoustic coupling
+ if( ACOUSTIC_SIMULATION ) then
+ call compute_coupling_poroelastic_ac(NSPEC_AB,NGLOB_AB, &
+ ibool,accels_poroelastic,accelw_poroelastic, &
+ potential_dot_dot_acoustic, &
+ num_coupling_ac_po_faces, &
+ coupling_ac_po_ispec,coupling_ac_po_ijk, &
+ coupling_ac_po_normal, &
+ coupling_ac_po_jacobian2Dw, &
+ rhoarraystore,phistore,tortstore, &
+ ispec_is_inner,phase_is_inner)
+
+ ! adjoint simulations
+ !if( SIMULATION_TYPE == 3 ) &
+! 'adjoint acoustic-poroelastic simulation not implemented yet'
+! call compute_coupling_elastic_ac(NSPEC_ADJOINT,NGLOB_ADJOINT, &
+! ibool,b_accel,b_potential_dot_dot_acoustic, &
+! num_coupling_ac_el_faces, &
+! coupling_ac_el_ispec,coupling_ac_el_ijk, &
+! coupling_ac_el_normal, &
+! coupling_ac_el_jacobian2Dw, &
+! ispec_is_inner,phase_is_inner)
+ endif
+
+! elastic coupling
+!chris: TO DO
+ if( ELASTIC_SIMULATION ) &
+ stop 'elastic-poroelastic simulation not implemented yet'
+! call compute_coupling_poroelastic_el()
+
+ ! adjoint simulations
+!chris: TO DO
+! if( SIMULATION_TYPE == 3 ) &
+! call compute_coupling_elastic_ac(NSPEC_ADJOINT,NGLOB_ADJOINT, &
+! ibool,b_accel,b_potential_dot_dot_acoustic, &
+! num_coupling_ac_el_faces, &
+! coupling_ac_el_ispec,coupling_ac_el_ijk, &
+! coupling_ac_el_normal, &
+! coupling_ac_el_jacobian2Dw, &
+! ispec_is_inner,phase_is_inner)
+! endif
+
+! adds source term (single-force/moment-tensor solution)
+ call compute_add_sources_poroelastic( NSPEC_AB,NGLOB_AB, &
+ accels_poroelastic,accelw_poroelastic,&
+ rhoarraystore,phistore,tortstore,&
+ ibool,ispec_is_inner,phase_is_inner, &
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,&
+ xi_source,eta_source,gamma_source,nu_source, &
+ hdur,hdur_gaussian,tshift_cmt,dt,t0,sourcearrays, &
+ ispec_is_poroelastic,SIMULATION_TYPE,NSTEP,NGLOB_ADJOINT, &
+ nrec,islice_selected_rec,ispec_selected_rec, &
+ nadj_rec_local,adj_sourcearrays)
+
+! assemble all the contributions between slices using MPI
+ if( phase_is_inner .eqv. .false. ) then
+ ! sends accel values to corresponding MPI interface neighbors
+ call assemble_MPI_vector_ext_mesh_po_s(NPROC,NGLOB_AB,accels_poroelastic, &
+ accelw_poroelastic,&
+ buffer_send_vector_ext_mesh_s,buffer_recv_vector_ext_mesh_s, &
+ buffer_send_vector_ext_mesh_w,buffer_recv_vector_ext_mesh_w, &
+ num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
+ my_neighbours_ext_mesh, &
+ request_send_vector_ext_mesh_s,request_recv_vector_ext_mesh_s, &
+ request_send_vector_ext_mesh_w,request_recv_vector_ext_mesh_w)
+
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) then
+! 'adjoint poroelastic simulation not implemented yet'
+! call assemble_MPI_vector_ext_mesh_s(NPROC,NGLOB_ADJOINT,b_accel, &
+! b_buffer_send_vector_ext_mesh,b_buffer_recv_vector_ext_mesh, &
+! num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
+! nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,&
+! my_neighbours_ext_mesh, &
+! b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh)
+ endif !adjoint
+
+ else
+ ! waits for send/receive requests to be completed and assembles values
+! solid phase
+ call assemble_MPI_vector_ext_mesh_po_w(NPROC,NGLOB_AB,accels_poroelastic, &
+ accelw_poroelastic,&
+ buffer_recv_vector_ext_mesh_s,buffer_recv_vector_ext_mesh_w, &
+ num_interfaces_ext_mesh,&
+ max_nibool_interfaces_ext_mesh, &
+ nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+ request_send_vector_ext_mesh_s,request_recv_vector_ext_mesh_s, &
+ request_send_vector_ext_mesh_w,request_recv_vector_ext_mesh_w)
+
+ ! adjoint simulations
+ if( SIMULATION_TYPE == 3 ) then
+! 'adjoint poroelastic simulation not implemented yet'
+! call assemble_MPI_vector_ext_mesh_w(NPROC,NGLOB_ADJOINT,b_accel, &
+! b_buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh,&
+! max_nibool_interfaces_ext_mesh, &
+! nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
+! b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh)
+ 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.
+
+ enddo
+
+! solid phase
+! multiplies with inverse of mass matrix (note: rmass has been inverted already)
+ accels_poroelastic(1,:) = accels_poroelastic(1,:)*rmass_solid_poroelastic(:)
+ accels_poroelastic(2,:) = accels_poroelastic(2,:)*rmass_solid_poroelastic(:)
+ accels_poroelastic(3,:) = accels_poroelastic(3,:)*rmass_solid_poroelastic(:)
+
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+! 'adjoint poroelastic simulation not implemented yet'
+! b_accels_poroelastic(1,:) = b_accels_poroelastic(1,:)*rmass_solid_poroelastic(:)
+! b_accels_poroelastic(2,:) = b_accels_poroelastic(2,:)*rmass_solid_poroelastic(:)
+! b_accels_poroelastic(3,:) = b_accels_poroelastic(3,:)*rmass_solid_poroelastic(:)
+ endif !adjoint
+
+! fluid phase
+! multiplies with inverse of mass matrix (note: rmass has been inverted already)
+ accelw_poroelastic(1,:) = accelw_poroelastic(1,:)*rmass_fluid_poroelastic(:)
+ accelw_poroelastic(2,:) = accelw_poroelastic(2,:)*rmass_fluid_poroelastic(:)
+ accelw_poroelastic(3,:) = accelw_poroelastic(3,:)*rmass_fluid_poroelastic(:)
+
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+! 'adjoint poroelastic simulation not implemented yet'
+! b_accelw_poroelastic(1,:) = b_accelw_poroelastic(1,:)*rmass_fluid_poroelastic(:)
+! b_accelw_poroelastic(2,:) = b_accelw_poroelastic(2,:)*rmass_fluid_poroelastic(:)
+! b_accelw_poroelastic(3,:) = b_accelw_poroelastic(3,:)*rmass_fluid_poroelastic(:)
+ endif !adjoint
+
+! updates velocities
+! Newark finite-difference time scheme with elastic domains:
+! (see e.g. Hughes, 1987; Chaljub et al., 2003)
+!
+! u(t+delta_t) = u(t) + delta_t v(t) + 1/2 delta_t**2 a(t)
+! v(t+delta_t) = v(t) + 1/2 delta_t a(t) + 1/2 delta_t a(t+delta_t)
+! a(t+delta_t) = 1/M_elastic ( -K_elastic u(t+delta) + B_elastic chi_dot_dot(t+delta_t) + f( t+delta_t) )
+!
+! where
+! u, v, a are displacement,velocity & acceleration
+! M is mass matrix, K stiffness matrix and B boundary term for acoustic/elastic domains
+! f denotes a source term (acoustic/elastic)
+! chi_dot_dot is acoustic (fluid) potential ( dotted twice with respect to time)
+!
+! corrector:
+! updates the velocity term which requires a(t+delta)
+! solid phase
+ velocs_poroelastic(:,:) = velocs_poroelastic(:,:) + deltatover2*accels_poroelastic(:,:)
+
+! fluid phase
+ velocw_poroelastic(:,:) = velocw_poroelastic(:,:) + deltatover2*accelw_poroelastic(:,:)
+
+ ! adjoint simulations
+! solid phase
+! if (SIMULATION_TYPE == 3) b_velocs_poroelastic(:,:) = b_velocs_poroelastic(:,:) + &
+! b_deltatover2*b_accels_poroelastic(:,:)
+
+! fluid phase
+! if (SIMULATION_TYPE == 3) b_velocw_poroelastic(:,:) = b_velocw_poroelastic(:,:) + &
+! b_deltatover2*b_accelw_poroelastic(:,:)
+
+end subroutine compute_forces_poroelastic
+
Added: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_solid.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_solid.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_solid.f90 2011-10-06 02:26:04 UTC (rev 19024)
@@ -0,0 +1,599 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine compute_forces_solid( iphase, &
+ NSPEC_AB,NGLOB_AB,displs_poroelastic,accels_poroelastic,&
+ velocs_poroelastic,displw_poroelastic,velocw_poroelastic,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz,&
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,&
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wxgll,wygll,wzgll, &
+ kappaarraystore,rhoarraystore,mustore,etastore,permstore, &
+ phistore,tortstore,jacobian,ibool,&
+ SIMULATION_TYPE,NGLOB_ADJOINT,NSPEC_ADJOINT, &
+ num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic,&
+ phase_ispec_inner_poroelastic )
+
+! compute forces for the solid poroelastic part
+
+ use constants,only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM, &
+ N_SLS, &
+ ONE_THIRD,FOUR_THIRDS
+
+ implicit none
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displs_poroelastic,accels_poroelastic,&
+ velocs_poroelastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displw_poroelastic,velocw_poroelastic
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ mustore,etastore,phistore,tortstore,jacobian
+ real(kind=CUSTOM_REAL), dimension(2,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rhoarraystore
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: kappaarraystore
+ real(kind=CUSTOM_REAL), dimension(6,NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: permstore
+
+! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: wxgll
+ double precision, dimension(NGLLY) :: wygll
+ double precision, dimension(NGLLZ) :: wzgll
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+! logical,dimension(NSPEC_AB) :: ispec_is_elastic
+ integer :: iphase
+ integer :: num_phase_ispec_poroelastic,nspec_inner_poroelastic,nspec_outer_poroelastic
+ integer, dimension(num_phase_ispec_poroelastic,2) :: phase_ispec_inner_poroelastic
+
+
+! adjoint simulations
+ integer :: SIMULATION_TYPE
+ integer :: NSPEC_BOUN
+ integer :: NGLOB_ADJOINT,NSPEC_ADJOINT
+! adjoint wavefields
+! real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT):: b_displs_poroelastic,b_accels_poroelastic
+! real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT):: b_displw_poroelastic
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: &
+! mufr_kl, B_kl
+
+
+!---
+!--- local variables
+!---
+
+ integer :: ispec,i,j,k,l,iglob,num_elements,ispec_p
+
+! spatial derivatives
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+ tempx1p,tempx2p,tempx3p,tempy1p,tempy2p,tempy3p,tempz1p,tempz2p,tempz3p
+! b_tempx1,b_tempx2,b_tempx3,b_tempy1,b_tempy2,b_tempy3,b_tempz1,b_tempz2,b_tempz3, &
+! b_tempx1p,b_tempx2p,b_tempx3p,b_tempy1p,b_tempy2p,b_tempy3p,b_tempz1p,b_tempz2p,b_tempz3p
+
+ real(kind=CUSTOM_REAL) :: duxdxl,duydxl,duzdxl,duxdyl,duydyl,duzdyl,duxdzl,duydzl,duzdzl
+ real(kind=CUSTOM_REAL) :: dwxdxl,dwydxl,dwzdxl,dwxdyl,dwydyl,dwzdyl,dwxdzl,dwydzl,dwzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl_plus_duzdzl,dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
+
+ real(kind=CUSTOM_REAL) tempx1ls,tempx2ls,tempx3ls,tempx1lw,tempx2lw,tempx3lw
+ real(kind=CUSTOM_REAL) tempy1ls,tempy2ls,tempy3ls,tempy1lw,tempy2lw,tempy3lw
+ real(kind=CUSTOM_REAL) tempz1ls,tempz2ls,tempz3ls,tempz1lw,tempz2lw,tempz3lw
+
+! real(kind=CUSTOM_REAL) :: b_dux_dxl,b_duz_dxl,b_dux_dzl,b_duz_dzl
+! real(kind=CUSTOM_REAL) :: dsxx,dsxy,dsxz,dsyy,dsyz,dszz
+! real(kind=CUSTOM_REAL) :: b_dsxx,b_dsxy,b_dsxz,b_dsyy,b_dsyz,b_dszz
+! real(kind=CUSTOM_REAL) :: b_dwx_dxl,b_dwz_dxl,b_dwx_dzl,b_dwz_dzl
+! real(kind=CUSTOM_REAL) :: dwxx,dwxy,dwxz,dwyy,dwyz,dwzz
+! real(kind=CUSTOM_REAL) :: b_dwxx,b_dwxy,b_dwxz,b_dwyy,b_dwyz,b_dwzz
+ real(kind=CUSTOM_REAL) :: sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+ real(kind=CUSTOM_REAL) :: sigmap
+! real(kind=CUSTOM_REAL) :: b_sigma_xx,b_sigma_yy,b_sigma_zz,b_sigma_xy,b_sigma_xz,b_sigma_yz
+! real(kind=CUSTOM_REAL) :: b_sigmap
+! real(kind=CUSTOM_REAL) :: nx,nz,vx,vz,vn,vxf,vzf,vnf,rho_vpI,rho_vpII,rho_vs,tx,tz,weight,xxi,zxi,xgamma,zgamma,jacobian1D
+
+! viscous attenuation (poroelastic media)
+ real(kind=CUSTOM_REAL), dimension(6) :: bl_relaxed
+
+
+! Jacobian matrix and determinant
+ real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+
+! material properties of the poroelastic medium
+! real(kind=CUSTOM_REAL) :: mul_unrelaxed,lambdal_unrelaxed,lambdalplus2mul_unrelaxed
+ real(kind=CUSTOM_REAL) :: kappal_s,rhol_s
+ real(kind=CUSTOM_REAL) :: etal_f,kappal_f,rhol_f
+ real(kind=CUSTOM_REAL) :: mul_fr,kappal_fr,phil,tortl,viscodampx,viscodampy,viscodampz
+ real(kind=CUSTOM_REAL) :: permlxx,permlxy,permlxz,permlyz,permlyy,permlzz,&
+ invpermlxx,invpermlxy,invpermlxz,invpermlyz,invpermlyy,invpermlzz,detk
+ real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,rhol_bar
+
+ real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
+! real(kind=CUSTOM_REAL) :: cpIsquare,cpIIsquare,cssquare,cpIl,cpIIl,csl
+
+! for attenuation
+! real(kind=CUSTOM_REAL) :: Un,Unp1,tauinv,Sn,Snp1,theta_n,theta_np1,tauinvsquare,tauinvcube,tauinvUn
+
+! compute Grad(displs_poroelastic) at time step n for attenuation
+! if(TURN_ATTENUATION_ON) call compute_gradient_attenuation(displs_poroelastic,dux_dxl_n,duz_dxl_n, &
+! dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,poroelastic,hprime_xx,hprime_zz,nspec,npoin)
+
+
+ if( iphase == 1 ) then
+ num_elements = nspec_outer_poroelastic
+ else
+ num_elements = nspec_inner_poroelastic
+ endif
+
+! loop over spectral elements
+ do ispec_p = 1,num_elements
+
+ ispec = phase_ispec_inner_poroelastic(ispec_p,iphase)
+
+! first double loop over GLL points to compute and store gradients
+ do k=1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+! get poroelastic parameters of current local GLL
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+!solid properties
+ kappal_s = kappaarraystore(1,i,j,k,ispec)
+ rhol_s = rhoarraystore(1,i,j,k,ispec)
+!fluid properties
+ kappal_f = kappaarraystore(2,i,j,k,ispec)
+ rhol_f = rhoarraystore(2,i,j,k,ispec)
+!frame properties
+ mul_fr = mustore(i,j,k,ispec)
+ kappal_fr = kappaarraystore(3,i,j,k,ispec)
+ rhol_bar = (1._CUSTOM_REAL - phil)*rhol_s + phil*rhol_f
+!Biot coefficients for the input phi
+ 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)
+!The RHS has the form : div T -phi/c div T_f + phi/ceta_fk^-1.partial t w
+!where T = G:grad u_s + C_biot div w I
+!and T_f = C_biot div u_s I + M_biot div w I
+ mul_G = mul_fr
+ lambdal_G = H_biot - 2._CUSTOM_REAL*mul_fr
+ lambdalplus2mul_G = lambdal_G + 2._CUSTOM_REAL*mul_G
+
+! derivative along x,y,z for u_s and w
+ tempx1ls = 0.
+ tempx2ls = 0.
+ tempx3ls = 0.
+
+ tempy1ls = 0.
+ tempy2ls = 0.
+ tempy3ls = 0.
+
+ tempz1ls = 0.
+ tempz2ls = 0.
+ tempz3ls = 0.
+
+ tempx1lw = 0.
+ tempx2lw = 0.
+ tempx3lw = 0.
+
+ tempy1lw = 0.
+ tempy2lw = 0.
+ tempy3lw = 0.
+
+ tempz1lw = 0.
+ tempz2lw = 0.
+ tempz3lw = 0.
+
+ if (SIMULATION_TYPE == 3) then
+ endif
+
+! first double loop over GLL points to compute and store gradients
+ do l = 1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1ls = tempx1ls + displs_poroelastic(1,iglob)*hp1
+ tempy1ls = tempy1ls + displs_poroelastic(2,iglob)*hp1
+ tempz1ls = tempz1ls + displs_poroelastic(3,iglob)*hp1
+ tempx1lw = tempx1lw + displw_poroelastic(1,iglob)*hp1
+ tempy1lw = tempy1lw + displw_poroelastic(2,iglob)*hp1
+ tempz1lw = tempz1lw + displw_poroelastic(3,iglob)*hp1
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2ls = tempx2ls + displs_poroelastic(1,iglob)*hp2
+ tempy2ls = tempy2ls + displs_poroelastic(2,iglob)*hp2
+ tempz2ls = tempz2ls + displs_poroelastic(3,iglob)*hp2
+ tempx2lw = tempx2lw + displw_poroelastic(1,iglob)*hp2
+ tempy2lw = tempy2lw + displw_poroelastic(2,iglob)*hp2
+ tempz2lw = tempz2lw + displw_poroelastic(3,iglob)*hp2
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3ls = tempx3ls + displs_poroelastic(1,iglob)*hp3
+ tempy3ls = tempy3ls + displs_poroelastic(2,iglob)*hp3
+ tempz3ls = tempz3ls + displs_poroelastic(3,iglob)*hp3
+ tempx3lw = tempx3lw + displw_poroelastic(1,iglob)*hp3
+ tempy3lw = tempy3lw + displw_poroelastic(2,iglob)*hp3
+ tempz3lw = tempz3lw + displw_poroelastic(3,iglob)*hp3
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif ! adjoint
+ enddo
+
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+ jacobianl = jacobian(i,j,k,ispec)
+
+! derivatives of displacement
+ duxdxl = xixl*tempx1ls + etaxl*tempx2ls + gammaxl*tempx3ls
+ duxdyl = xiyl*tempx1ls + etayl*tempx2ls + gammayl*tempx3ls
+ duxdzl = xizl*tempx1ls + etazl*tempx2ls + gammazl*tempx3ls
+
+ duydxl = xixl*tempy1ls + etaxl*tempy2ls + gammaxl*tempy3ls
+ duydyl = xiyl*tempy1ls + etayl*tempy2ls + gammayl*tempy3ls
+ duydzl = xizl*tempy1ls + etazl*tempy2ls + gammazl*tempy3ls
+
+ duzdxl = xixl*tempz1ls + etaxl*tempz2ls + gammaxl*tempz3ls
+ duzdyl = xiyl*tempz1ls + etayl*tempz2ls + gammayl*tempz3ls
+ duzdzl = xizl*tempz1ls + etazl*tempz2ls + gammazl*tempz3ls
+
+ dwxdxl = xixl*tempx1lw + etaxl*tempx2lw + gammaxl*tempx3lw
+ dwxdyl = xiyl*tempx1lw + etayl*tempx2lw + gammayl*tempx3lw
+ dwxdzl = xizl*tempx1lw + etazl*tempx2lw + gammazl*tempx3lw
+
+ dwydxl = xixl*tempy1lw + etaxl*tempy2lw + gammaxl*tempy3lw
+ dwydyl = xiyl*tempy1lw + etayl*tempy2lw + gammayl*tempy3lw
+ dwydzl = xizl*tempy1lw + etazl*tempy2lw + gammazl*tempy3lw
+
+ dwzdxl = xixl*tempz1lw + etaxl*tempz2lw + gammaxl*tempz3lw
+ dwzdyl = xiyl*tempz1lw + etayl*tempz2lw + gammayl*tempz3lw
+ dwzdzl = xizl*tempz1lw + etazl*tempz2lw + gammazl*tempz3lw
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl_plus_duzdzl = duxdxl + duydyl + duzdzl
+ dwxdxl_plus_dwydyl_plus_dwzdzl = dwxdxl + dwydyl + dwzdzl
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! compute stress tensor (include attenuation or anisotropy if needed)
+
+! if(VISCOATTENUATION) then
+!chris:check
+
+! Dissipation only controlled by frame share attenuation in poroelastic (see Morency & Tromp, GJI 2008).
+! attenuation is implemented following the memory variable formulation of
+! J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
+! vol. 58(1), p. 110-120 (1993). More details can be found in
+! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a linear
+! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611 (1988).
+
+! else
+
+! no attenuation
+ sigma_xx = lambdalplus2mul_G*duxdxl + lambdal_G*duydyl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_yy = lambdalplus2mul_G*duydyl + lambdal_G*duxdxl_plus_duzdzl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+ sigma_zz = lambdalplus2mul_G*duzdzl + lambdal_G*duxdxl_plus_duydyl + C_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ sigma_xy = mul_G*duxdyl_plus_duydxl
+ sigma_xz = mul_G*duzdxl_plus_duxdzl
+ sigma_yz = mul_G*duzdyl_plus_duydzl
+
+ sigmap = C_biot*duxdxl_plus_duydyl_plus_duzdzl + M_biot*dwxdxl_plus_dwydyl_plus_dwzdzl
+
+ if(SIMULATION_TYPE == 3) then ! kernels calculation
+ endif
+! endif !if(ATTENUATION)
+
+! 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
+
+
+! weak formulation term based on stress tensor (non-symmetric form)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ tempx1p(i,j,k) = jacobianl * sigmap*xixl
+ tempy1p(i,j,k) = jacobianl * sigmap*xiyl
+ tempz1p(i,j,k) = jacobianl * sigmap*xizl
+
+ tempx2p(i,j,k) = jacobianl * sigmap*etaxl
+ tempy2p(i,j,k) = jacobianl * sigmap*etayl
+ tempz2p(i,j,k) = jacobianl * sigmap*etazl
+
+ tempx3p(i,j,k) = jacobianl * sigmap*gammaxl
+ tempy3p(i,j,k) = jacobianl * sigmap*gammayl
+ tempz3p(i,j,k) = jacobianl * sigmap*gammazl
+
+ if(SIMULATION_TYPE == 3) then ! kernels calculation
+ endif
+
+ enddo
+ enddo
+ enddo
+
+!
+! second double-loop over GLL to compute all the terms
+!
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ tempx1ls = 0.
+ tempy1ls = 0.
+ tempz1ls = 0.
+
+ tempx2ls = 0.
+ tempy2ls = 0.
+ tempz2ls = 0.
+
+ tempx3ls = 0.
+ tempy3ls = 0.
+ tempz3ls = 0.
+
+ tempx1lw = 0.
+ tempy1lw = 0.
+ tempz1lw = 0.
+
+ tempx2lw = 0.
+ tempy2lw = 0.
+ tempz2lw = 0.
+
+ tempx3lw = 0.
+ tempy3lw = 0.
+ tempz3lw = 0.
+
+ if (SIMULATION_TYPE == 3) then
+ endif
+
+ do l=1,NGLLX
+ fac1 = hprimewgll_xx(l,i)
+ tempx1ls = tempx1ls + tempx1(l,j,k)*fac1
+ tempy1ls = tempy1ls + tempy1(l,j,k)*fac1
+ tempz1ls = tempz1ls + tempz1(l,j,k)*fac1
+ tempx1lw = tempx1lw + tempx1p(l,j,k)*fac1
+ tempy1lw = tempy1lw + tempy1p(l,j,k)*fac1
+ tempz1lw = tempz1lw + tempz1p(l,j,k)*fac1
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ fac2 = hprimewgll_yy(l,j)
+ tempx2ls = tempx2ls + tempx2(i,l,k)*fac2
+ tempy2ls = tempy2ls + tempy2(i,l,k)*fac2
+ tempz2ls = tempz2ls + tempz2(i,l,k)*fac2
+ tempx2lw = tempx2lw + tempx2p(i,l,k)*fac2
+ tempy2lw = tempy2lw + tempy2p(i,l,k)*fac2
+ tempz2lw = tempz2lw + tempz2p(i,l,k)*fac2
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+ !!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ fac3 = hprimewgll_zz(l,k)
+ tempx3ls = tempx3ls + tempx3(i,j,l)*fac3
+ tempy3ls = tempy3ls + tempy3(i,j,l)*fac3
+ tempz3ls = tempz3ls + tempz3(i,j,l)*fac3
+ tempx3lw = tempx3lw + tempx3p(i,j,l)*fac3
+ tempy3lw = tempy3lw + tempy3p(i,j,l)*fac3
+ tempz3lw = tempz3lw + tempz3p(i,j,l)*fac3
+ ! adjoint simulations
+ if (SIMULATION_TYPE == 3) then
+ endif
+ enddo
+
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+ phil = phistore(i,j,k,ispec)
+ tortl = tortstore(i,j,k,ispec)
+
+ ! sum contributions from each element to the global mesh
+
+ iglob = ibool(i,j,k,ispec)
+
+
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) - ( fac1*(tempx1ls - phil/tortl*tempx1lw) &
+ + fac2*(tempx2ls - phil/tortl*tempx2lw) + fac3*(tempx3ls - phil/tortl*tempx3lw) )
+
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) - ( fac1*(tempy1ls - phil/tortl*tempy1lw) &
+ + fac2*(tempy2ls - phil/tortl*tempy2lw) + fac3*(tempy3ls - phil/tortl*tempy3lw) )
+
+ 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
+!
+! add + phi/tort eta_f k^-1 dot(w)
+
+ etal_f = etastore(i,j,k,ispec)
+
+ if(etal_f >0.d0) then
+
+ 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)
+
+! calcul of the inverse of k
+ detk = permlxz*(permlxy*permlyz-permlxz*permlyy) &
+ - permlxy*(permlxy*permlzz-permlyz*permlxz) &
+ + permlxx*(permlyy*permlzz-permlyz*permlyz)
+
+ if(detk /= 0.d0) then
+ invpermlxx = (permlyy*permlzz-permlyz*permlyz)/detk
+ invpermlxy = (permlxz*permlyz-permlxy*permlzz)/detk
+ invpermlxz = (permlxy*permlyz-permlxz*permlyy)/detk
+ invpermlyy = (permlxx*permlzz-permlxz*permlxz)/detk
+ invpermlyz = (permlxy*permlxz-permlxx*permlyz)/detk
+ invpermlzz = (permlxx*permlyy-permlxy*permlxy)/detk
+ else
+ stop 'Permeability matrix is not inversible'
+ endif
+
+! relaxed viscous coef
+ bl_relaxed(1) = etal_f*invpermlxx
+ bl_relaxed(2) = etal_f*invpermlxy
+ bl_relaxed(3) = etal_f*invpermlxz
+ bl_relaxed(4) = etal_f*invpermlyy
+ bl_relaxed(5) = etal_f*invpermlyz
+ bl_relaxed(6) = etal_f*invpermlzz
+
+! ifVISCOATTENUATION) 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
+! endif
+
+! do k = 1,NGLLZ
+! do j = 1,NGLLY
+! do i = 1,NGLLX
+
+! iglob = ibool(i,j,k,ispec)
+
+! if(VISCOATTENUATION) then
+! compute the viscous damping term with the unrelaxed viscous coef and add memory variable
+! viscodampx = velocw_poroelastic(1,iglob)*bl_unrelaxed(1) + velocw_poroelastic(2,iglob)*bl_unrelaxed(2)&
+! - rx_viscous(i,j,ispec)
+! viscodampz = velocw_poroelastic(1,iglob)*bl_unrelaxed(2) + velocw_poroelastic(2,iglob)*bl_unrelaxed(3)&
+! - rz_viscous(i,j,ispec)
+! else
+
+! no viscous attenuation
+ viscodampx = velocw_poroelastic(1,iglob)*bl_relaxed(1) + velocw_poroelastic(2,iglob)*bl_relaxed(2) + &
+ velocw_poroelastic(3,iglob)*bl_relaxed(3)
+ viscodampy = velocw_poroelastic(1,iglob)*bl_relaxed(2) + velocw_poroelastic(2,iglob)*bl_relaxed(4) + &
+ velocw_poroelastic(3,iglob)*bl_relaxed(5)
+ viscodampz = velocw_poroelastic(1,iglob)*bl_relaxed(3) + velocw_poroelastic(2,iglob)*bl_relaxed(5) + &
+ velocw_poroelastic(3,iglob)*bl_relaxed(6)
+! endif
+
+ accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) + phil/tortl*wxgll(i)*wygll(j)*wzgll(k)*jacobian(i,j,k,ispec)*&
+ viscodampx
+ accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) + phil/tortl*wxgll(i)*wygll(j)*wzgll(k)*jacobian(i,j,k,ispec)*&
+ viscodampy
+ accels_poroelastic(3,iglob) = accels_poroelastic(3,iglob) + phil/tortl*wxgll(i)*wygll(j)*wzgll(k)*jacobian(i,j,k,ispec)*&
+ viscodampz
+
+! if isolver == 1 .and. save_forward then b_viscodamp is save in compute_forces_fluid.f90
+! if(isolver == 2) then ! kernels calculation
+! b_accels_poroelastic(1,iglob) = b_accels_poroelastic(1,iglob) + phil/tortl*b_viscodampx(iglob)
+! b_accels_poroelastic(2,iglob) = b_accels_poroelastic(2,iglob) + phil/tortl*b_viscodampz(iglob)
+! endif
+
+! enddo
+! enddo
+! enddo
+
+ endif ! if(etal_f >0.d0) then
+
+ enddo ! second loop over the GLL points
+ enddo
+ enddo
+
+ enddo ! end of loop over all spectral elements
+
+
+ end subroutine compute_forces_solid
+
More information about the CIG-COMMITS
mailing list