[cig-commits] [commit] devel: Add the new PML formulation for fluid-solid simulation (64a2a90)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Sat Aug 9 06:24:25 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/1e74fc10b34d28f5c1879591f5c7c639e72fa24c...2e2f3047c3b54ead260dc80b055ddab66a20bb01
>---------------------------------------------------------------
commit 64a2a907f6c9a3ba68d1fd8d69aac7665bcf8d6d
Author: Xie Zhinan <xiezhinan1984 at gmail.com>
Date: Sat Aug 9 20:29:02 2014 +0800
Add the new PML formulation for fluid-solid simulation
>---------------------------------------------------------------
64a2a907f6c9a3ba68d1fd8d69aac7665bcf8d6d
src/specfem2D/Makefile.in | 3 +
...el.f90 => compute_coupling_viscoelastic_ac.f90} | 133 ++++++++++++---------
src/specfem2D/specfem2D.F90 | 127 +++-----------------
3 files changed, 94 insertions(+), 169 deletions(-)
diff --git a/src/specfem2D/Makefile.in b/src/specfem2D/Makefile.in
index 3e1e924..119eb24 100644
--- a/src/specfem2D/Makefile.in
+++ b/src/specfem2D/Makefile.in
@@ -108,6 +108,7 @@ OBJS_SPECFEM2D = \
$O/compute_curl_one_element.o \
$O/compute_energy.o \
$O/compute_coupling_acoustic_el.o \
+ $O/compute_coupling_viscoelastic_ac.o \
$O/compute_forces_gravitoacoustic.o \
$O/compute_forces_acoustic.o \
$O/compute_forces_viscoelastic.o \
@@ -427,6 +428,8 @@ $O/check_stability.o: ${S}/check_stability.F90 ${SETUP}/constants.h
$O/compute_coupling_acoustic_el.o: ${S}/compute_coupling_acoustic_el.f90 ${SETUP}/constants.h
${F90} $(DEF_FFLAGS) -c -o $O/compute_coupling_acoustic_el.o ${S}/compute_coupling_acoustic_el.f90
+$O/compute_coupling_viscoelastic_ac.o: ${S}/compute_coupling_viscoelastic_ac.f90 ${SETUP}/constants.h
+ ${F90} $(DEF_FFLAGS) -c -o $O/compute_coupling_viscoelastic_ac.o ${S}/compute_coupling_viscoelastic_ac.f90
$O/compute_energy.o: ${S}/compute_energy.f90 ${SETUP}/constants.h
${F90} $(DEF_FFLAGS) -c -o $O/compute_energy.o ${S}/compute_energy.f90
diff --git a/src/specfem2D/compute_coupling_acoustic_el.f90 b/src/specfem2D/compute_coupling_viscoelastic_ac.f90
similarity index 59%
copy from src/specfem2D/compute_coupling_acoustic_el.f90
copy to src/specfem2D/compute_coupling_viscoelastic_ac.f90
index edfd4e0..0eb3ddc 100644
--- a/src/specfem2D/compute_coupling_acoustic_el.f90
+++ b/src/specfem2D/compute_coupling_viscoelastic_ac.f90
@@ -3,11 +3,12 @@
! S P E C F E M 2 D Version 7 . 0
! --------------------------------
!
-! Main historical authors: Dimitri Komatitsch and Jeroen Tromp
-! Princeton University, USA
-! and CNRS / University of Marseille, France
-! (there are currently many more authors!)
-! (c) Princeton University and CNRS / University of Marseille, April 2014
+! Copyright CNRS, INRIA and University of Pau, France,
+! and Princeton University / California Institute of Technology, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
!
! This software is a computer program whose purpose is to solve
! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
@@ -16,7 +17,7 @@
! This software is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software. You can use,
! modify and/or redistribute the software under the terms of the CeCILL
-! license as circulated by CEA, CNRS and Inria at the following URL
+! license as circulated by CEA, CNRS and INRIA at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and rights to copy,
@@ -42,23 +43,26 @@
! for acoustic solver
- subroutine compute_coupling_acoustic_el(nspec,nglob_elastic,nglob_acoustic,num_fluid_solid_edges,ibool,wxgll,wzgll,xix,xiz,&
- gammax,gammaz,jacobian,ivalue,jvalue,ivalue_inverse,jvalue_inverse,displ_elastic,displ_elastic_old,&
- potential_dot_dot_acoustic,fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge, &
- fluid_solid_elastic_ispec,fluid_solid_elastic_iedge,&
- AXISYM,nglob,coord,is_on_the_axis,xiglj,wxglj, &
+ subroutine compute_coupling_viscoelastic_ac(SIMULATION_TYPE,nspec,nglob_elastic,nglob_acoustic,num_fluid_solid_edges,&
+ ibool,wxgll,wzgll,xix,xiz,gammax,gammaz,jacobian,ivalue,jvalue,ivalue_inverse,jvalue_inverse,&
+ potential_acoustic,potential_acoustic_old,potential_dot_acoustic,potential_dot_dot_acoustic,&
+ b_potential_dot_dot_acoustic,accel_elastic,b_accel_elastic,fluid_solid_acoustic_ispec, &
+ fluid_solid_acoustic_iedge,fluid_solid_elastic_ispec,fluid_solid_elastic_iedge,&
+ potential_acoustic_adj_coupling,AXISYM,nglob,coord,is_on_the_axis,xiglj,wxglj, &
PML_BOUNDARY_CONDITIONS,nspec_PML,K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,&
- alpha_z_store,is_PML,spec_to_PML,region_CPML,rmemory_fsb_displ_elastic,timeval,deltat)
+ alpha_z_store,is_PML,spec_to_PML,region_CPML,rmemory_sfb_potential_ddot_acoustic,timeval,deltat)
implicit none
include 'constants.h'
- integer :: nspec,nglob_elastic,nglob_acoustic,num_fluid_solid_edges
+ integer :: SIMULATION_TYPE,nspec,nglob_elastic,nglob_acoustic,num_fluid_solid_edges
integer :: nglob
logical :: AXISYM
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll,wzgll
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
+ integer, dimension(NGLLX,NEDGES) :: ivalue,jvalue,ivalue_inverse,jvalue_inverse
! Gauss-Lobatto-Jacobi points and weights
double precision, dimension(NGLJ) :: xiglj
@@ -67,11 +71,10 @@
double precision, dimension(NDIM,nglob), intent(in) :: coord
real(kind=CUSTOM_REAL), dimension(NGLJ,NGLLZ) :: r_xiplus1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
- integer, dimension(NGLLX,NEDGES) :: ivalue,jvalue,ivalue_inverse,jvalue_inverse
-
- real(kind=CUSTOM_REAL),dimension(3,nglob_elastic) :: displ_elastic,displ_elastic_old
- real(kind=CUSTOM_REAL),dimension(nglob_acoustic) :: potential_dot_dot_acoustic
+ real(kind=CUSTOM_REAL),dimension(3,nglob_elastic) :: accel_elastic,b_accel_elastic
+ real(kind=CUSTOM_REAL),dimension(nglob_acoustic) :: potential_acoustic,potential_acoustic_old,potential_dot_acoustic,&
+ potential_dot_dot_acoustic,b_potential_dot_dot_acoustic,&
+ potential_acoustic_adj_coupling
integer, dimension(num_fluid_solid_edges) :: fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge, &
fluid_solid_elastic_ispec,fluid_solid_elastic_iedge
integer :: nspec_PML
@@ -81,20 +84,20 @@
logical, dimension(nspec) :: is_PML
integer, dimension(nspec) :: spec_to_PML
integer, dimension(nspec) :: region_CPML
- real(kind=CUSTOM_REAL),dimension(1,3,NGLLX,NGLLZ,num_fluid_solid_edges) :: rmemory_fsb_displ_elastic
+ real(kind=CUSTOM_REAL),dimension(1,NGLLX,NGLLZ,num_fluid_solid_edges) :: rmemory_sfb_potential_ddot_acoustic
!local variable
- integer :: inum,ispec_acoustic,ispec_elastic,iedge_acoustic,iedge_elastic,ipoin1D,i,j,iglob,&
- ispec_PML,CPML_region_local,singularity_type_xz
- real(kind=CUSTOM_REAL) :: displ_x,displ_z,displ_n,&
+ integer :: inum,ispec_acoustic,ispec_elastic,iedge_acoustic,iedge_elastic,ipoin1D,i,j,iglob,ii2,jj2,&
+ ispec_PML,CPML_region_local,singularity_type
+ real(kind=CUSTOM_REAL) :: pressure,b_pressure,&
xxi,zxi,xgamma,zgamma,jacobian1D,nx,nz,weight
double precision :: timeval,deltat
- double precision :: kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z, &
- A8,A9,A10,bb_xz_1,bb_xz_2,coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2
+ double precision :: kappa_x,kappa_z,d_x,d_z,alpha_x,alpha_z,beta_x,beta_z,&
+ A0,A1,A2,A3,A4,bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2
- ! loop on all the coupling edges
+ ! loop on all the coupling edges
do inum = 1,num_fluid_solid_edges
! get the edge of the acoustic element
@@ -108,15 +111,18 @@
! implement 1D coupling along the edge
do ipoin1D = 1,NGLLX
- ! get point values for the elastic side, which matches our side in the inverse direction
- i = ivalue_inverse(ipoin1D,iedge_elastic)
- j = jvalue_inverse(ipoin1D,iedge_elastic)
- iglob = ibool(i,j,ispec_elastic)
+ ! get point values for the acoustic side, which matches our side in the inverse direction
+ i = ivalue_inverse(ipoin1D,iedge_acoustic)
+ j = jvalue_inverse(ipoin1D,iedge_acoustic)
+ iglob = ibool(i,j,ispec_acoustic)
+
+ ! compute pressure on the fluid/solid edge
+ pressure = - potential_dot_dot_acoustic(iglob)
if(PML_BOUNDARY_CONDITIONS)then
- if(is_PML(ispec_elastic) .and. nspec_PML > 0) then
- ispec_PML = spec_to_PML(ispec_elastic)
- CPML_region_local = region_CPML(ispec_elastic)
+ if(is_PML(ispec_acoustic) .and. nspec_PML > 0) then
+ ispec_PML = spec_to_PML(ispec_acoustic)
+ CPML_region_local = region_CPML(ispec_acoustic)
if(CPML_region_local == CPML_X_ONLY)then
kappa_x = K_x_store(i,j,ispec_PML)
kappa_z = K_z_store(i,j,ispec_PML)
@@ -126,31 +132,36 @@
alpha_z = alpha_z_store(i,j,ispec_PML)
beta_x = alpha_x + d_x / kappa_x
beta_z = alpha_z + d_z / kappa_z
- call lik_parameter_computation(timeval,deltat,kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z,&
- CPML_region_local,13,A8,A9,A10,singularity_type_xz,bb_xz_1,bb_xz_2,&
- coef0_xz_1,coef1_xz_1,coef2_xz_1,coef0_xz_2,coef1_xz_2,coef2_xz_2)
- rmemory_fsb_displ_elastic(1,1,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,1,i,j,inum) + &
- coef1_xz_1 * displ_elastic(1,iglob) + coef2_xz_1 * displ_elastic_old(1,iglob)
- rmemory_fsb_displ_elastic(1,3,i,j,inum) = coef0_xz_1 * rmemory_fsb_displ_elastic(1,3,i,j,inum) + &
- coef1_xz_1 * displ_elastic(3,iglob) + coef2_xz_1 * displ_elastic_old(3,iglob)
- displ_x = A8 * displ_elastic(1,iglob) + A9 * rmemory_fsb_displ_elastic(1,1,i,j,inum)
- displ_z = A8 * displ_elastic(3,iglob) + A9 * rmemory_fsb_displ_elastic(1,3,i,j,inum)
+ call l_parameter_computation(timeval,deltat,kappa_x,beta_x,alpha_x,kappa_z,beta_z,alpha_z, &
+ CPML_region_local,A0,A1,A2,A3,A4,singularity_type,&
+ bb_1,coef0_1,coef1_1,coef2_1,bb_2,coef0_2,coef1_2,coef2_2)
+ rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) = &
+ coef0_1 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum) + &
+ coef1_1 * potential_acoustic(iglob) + coef2_1 * potential_acoustic_old(iglob)
+ pressure = - (A0 * potential_dot_dot_acoustic(iglob) + A1 * potential_dot_acoustic(iglob) + &
+ A2 * potential_acoustic(iglob) + A3 * rmemory_sfb_potential_ddot_acoustic(1,i,j,inum))
+
else
- stop 'PML currently does not support a fluid-solid boundary located in a PML that is not CPML_X_ONLY'
+ stop 'PML do not support the fluid-solid boundary not inside CPML_X_ONLY'
endif
else
- displ_x = displ_elastic(1,iglob)
- displ_z = displ_elastic(3,iglob)
+ pressure = - potential_dot_dot_acoustic(iglob)
endif
else
- displ_x = displ_elastic(1,iglob)
- displ_z = displ_elastic(3,iglob)
+ pressure = - potential_dot_dot_acoustic(iglob)
endif
- ! get point values for the acoustic side
- i = ivalue(ipoin1D,iedge_acoustic)
- j = jvalue(ipoin1D,iedge_acoustic)
- iglob = ibool(i,j,ispec_acoustic)
+ if(SIMULATION_TYPE == 3) then
+ b_pressure = - b_potential_dot_dot_acoustic(iglob)
+ !<YANGL
+ ! new definition of adjoint displacement and adjoint potential
+ pressure = potential_acoustic_adj_coupling(iglob)
+ !>YANGL
+ endif
+ ! get point values for the elastic side
+ ii2 = ivalue(ipoin1D,iedge_elastic)
+ jj2 = jvalue(ipoin1D,iedge_elastic)
+ iglob = ibool(ii2,jj2,ispec_elastic)
! compute the 1D Jacobian and the normal to the edge: for their expression see for instance
! O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics,
@@ -159,11 +170,11 @@
! Blackwell Science, page 110, equation (4.60).
if (AXISYM) then
- if (abs(coord(1,iglob)) < TINYVAL) then
+ if (abs(coord(1,ibool(i,j,ispec_acoustic))) < TINYVAL) then
xxi = + gammaz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
r_xiplus1(i,j) = xxi
else if (is_on_the_axis(ispec_acoustic)) then
- r_xiplus1(i,j) = coord(1,iglob)/(xiglj(i)+ONE)
+ r_xiplus1(i,j) = coord(1,ibool(i,j,ispec_acoustic))/(xiglj(i)+ONE)
endif
endif
@@ -177,7 +188,7 @@
if (is_on_the_axis(ispec_acoustic)) then
weight = jacobian1D * wxglj(i) * r_xiplus1(i,j)
else
- weight = jacobian1D * wxgll(i) * coord(1,iglob)
+ weight = jacobian1D * wxgll(i) * coord(1,ibool(i,j,ispec_acoustic))
endif
else
weight = jacobian1D * wxgll(i)
@@ -192,12 +203,12 @@
if (is_on_the_axis(ispec_acoustic)) then
weight = jacobian1D * wxglj(i) * r_xiplus1(i,j)
else
- weight = jacobian1D * wxgll(i) * coord(1,iglob)
+ weight = jacobian1D * wxgll(i) * coord(1,ibool(i,j,ispec_acoustic))
endif
else
weight = jacobian1D * wxgll(i)
endif
- else if(iedge_acoustic ==ILEFT)then
+ else if(iedge_acoustic == ILEFT)then
xgamma = - xiz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
zgamma = + xix(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
jacobian1D = sqrt(xgamma**2 + zgamma**2)
@@ -213,14 +224,18 @@
weight = jacobian1D * wzgll(j)
endif
- ! compute dot product
- displ_n = displ_x*nx + displ_z*nz
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + weight*nx*pressure
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + weight*nz*pressure
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) + weight*displ_n
+ if(SIMULATION_TYPE == 3) then
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + weight*nx*b_pressure
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + weight*nz*b_pressure
+ endif !if(SIMULATION_TYPE == 3) then
enddo
enddo
- end subroutine compute_coupling_acoustic_el
+ end subroutine compute_coupling_viscoelastic_ac
+
diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90
index e2117d9..bf38e33 100644
--- a/src/specfem2D/specfem2D.F90
+++ b/src/specfem2D/specfem2D.F90
@@ -1046,6 +1046,7 @@
rmemory_dux_dx,rmemory_duz_dx,rmemory_dux_dz,rmemory_duz_dz
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_fsb_displ_elastic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rmemory_sfb_potential_ddot_acoustic
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
rmemory_dux_dx_prime,rmemory_duz_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dz_prime
@@ -3162,6 +3163,8 @@
if(any_acoustic .and. num_fluid_solid_edges > 0)then
allocate(rmemory_fsb_displ_elastic(1,3,NGLLX,NGLLZ,num_fluid_solid_edges),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_fsb_displ_elastic'
+ allocate(rmemory_sfb_potential_ddot_acoustic(1,NGLLX,NGLLZ,num_fluid_solid_edges),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_sfb_potential_ddot_acoustic'
endif
if(ROTATE_PML_ACTIVATE)then
@@ -3211,6 +3214,7 @@
if(any_acoustic .and. num_fluid_solid_edges > 0)then
rmemory_fsb_displ_elastic(:,:,:,:,:) = ZERO
+ rmemory_sfb_potential_ddot_acoustic(:,:,:,:) = ZERO
endif
if(ROTATE_PML_ACTIVATE)then
@@ -3237,6 +3241,7 @@
allocate(rmemory_duz_dz(1,1,1,1))
if(any_acoustic .and. num_fluid_solid_edges > 0)then
allocate(rmemory_fsb_displ_elastic(1,3,NGLLX,NGLLZ,1))
+ allocate(rmemory_sfb_potential_ddot_acoustic(1,NGLLX,NGLLZ,1))
endif
allocate(rmemory_dux_dx_prime(1,1,1,1))
@@ -3292,6 +3297,7 @@
allocate(rmemory_duz_dx(1,1,1,1))
allocate(rmemory_duz_dz(1,1,1,1))
allocate(rmemory_fsb_displ_elastic(1,3,NGLLX,NGLLZ,1))
+ allocate(rmemory_sfb_potential_ddot_acoustic(1,NGLLX,NGLLZ,1))
allocate(rmemory_dux_dx_prime(1,1,1,1))
allocate(rmemory_dux_dz_prime(1,1,1,1))
@@ -3328,6 +3334,7 @@
! if so we need to allocate a dummy version in order to be able to use that array as an argument
! in some subroutine calls below
if(.not. allocated(rmemory_fsb_displ_elastic)) allocate(rmemory_fsb_displ_elastic(1,3,NGLLX,NGLLZ,1))
+ if(.not. allocated(rmemory_sfb_potential_ddot_acoustic)) allocate(rmemory_sfb_potential_ddot_acoustic(1,NGLLX,NGLLZ,1))
! Test compatibility with axisymmetric formulation
if(AXISYM) then
@@ -5390,8 +5397,6 @@ if(coupled_elastic_poro) then
if(SIMULATION_TYPE == 3)then
-! Zhinan Xie's new fluid-solid coupling formulation for adjoint runs
-! involves a change of sign compared to Yang Luo's initial formulation
accel_elastic_adj_coupling2 = - accel_elastic_adj_coupling
call compute_coupling_acoustic_el(nspec,nglob_elastic,nglob_acoustic,num_fluid_solid_edges,ibool,wxgll,wzgll,xix,xiz,&
@@ -5407,7 +5412,7 @@ if(coupled_elastic_poro) then
b_potential_dot_dot_acoustic,fluid_solid_acoustic_ispec,fluid_solid_acoustic_iedge, &
fluid_solid_elastic_ispec,fluid_solid_elastic_iedge,&
AXISYM,nglob,coord,is_on_the_axis,xiglj,wxglj, &
- PML_BOUNDARY_CONDITIONS,nspec_PML,K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,&
+ .false.,nspec_PML,K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,&
alpha_z_store,is_PML,spec_to_PML,region_CPML,rmemory_fsb_displ_elastic,timeval,deltat)
endif
@@ -5669,7 +5674,7 @@ if(coupled_elastic_poro) then
+ deltat * potential_dot_dot_acoustic
potential_acoustic_LDDRK = alpha_LDDRK(i_stage) * potential_acoustic_LDDRK &
- +deltat*potential_dot_acoustic
+ + deltat*potential_dot_acoustic
if(i_stage==1 .and. it == 1 .and. (.not. initialfield))then
!! DK DK this should be vectorized
@@ -6498,112 +6503,14 @@ if(coupled_elastic_poro) then
! *********************************************************
if(coupled_acoustic_elastic) then
-
- ! loop on all the coupling edges
- do inum = 1,num_fluid_solid_edges
-
- ! get the edge of the acoustic element
- ispec_acoustic = fluid_solid_acoustic_ispec(inum)
- iedge_acoustic = fluid_solid_acoustic_iedge(inum)
-
- ! get the corresponding edge of the elastic element
- ispec_elastic = fluid_solid_elastic_ispec(inum)
- iedge_elastic = fluid_solid_elastic_iedge(inum)
-
- ! implement 1D coupling along the edge
- do ipoin1D = 1,NGLLX
-
- ! get point values for the acoustic side, which matches our side in the inverse direction
- i = ivalue_inverse(ipoin1D,iedge_acoustic)
- j = jvalue_inverse(ipoin1D,iedge_acoustic)
- iglob = ibool(i,j,ispec_acoustic)
-
- ! compute pressure on the fluid/solid edge
- pressure = - potential_dot_dot_acoustic(iglob)
- if(SIMULATION_TYPE == 3) then
- b_pressure = - b_potential_dot_dot_acoustic(iglob)
- !<YANGL
- ! new definition of adjoint displacement and adjoint potential
- pressure = potential_acoustic_adj_coupling(iglob)
- !>YANGL
- endif
- ! get point values for the elastic side
- ii2 = ivalue(ipoin1D,iedge_elastic)
- jj2 = jvalue(ipoin1D,iedge_elastic)
- iglob = ibool(ii2,jj2,ispec_elastic)
-
- ! compute the 1D Jacobian and the normal to the edge: for their expression see for instance
- ! O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method for Solid and Structural Mechanics,
- ! Sixth Edition, electronic version, www.amazon.com, p. 204 and Figure 7.7(a),
- ! or Y. K. Cheung, S. H. Lo and A. Y. T. Leung, Finite Element Implementation,
- ! Blackwell Science, page 110, equation (4.60).
-
- if (AXISYM) then
- if (abs(coord(1,ibool(i,j,ispec_acoustic))) < TINYVAL) then
- xxi = + gammaz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- r_xiplus1(i,j) = xxi
- else if (is_on_the_axis(ispec_acoustic)) then
- r_xiplus1(i,j) = coord(1,ibool(i,j,ispec_acoustic))/(xiglj(i)+ONE)
- endif
- endif
-
- if(iedge_acoustic == ITOP)then
- xxi = + gammaz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- zxi = - gammax(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- jacobian1D = sqrt(xxi**2 + zxi**2)
- nx = - zxi / jacobian1D
- nz = + xxi / jacobian1D
- if (AXISYM) then
- if (is_on_the_axis(ispec_acoustic)) then
- weight = jacobian1D * wxglj(i) * r_xiplus1(i,j)
- else
- weight = jacobian1D * wxgll(i) * coord(1,ibool(i,j,ispec_acoustic))
- endif
- else
- weight = jacobian1D * wxgll(i)
- endif
- else if(iedge_acoustic == IBOTTOM)then
- xxi = + gammaz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- zxi = - gammax(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- jacobian1D = sqrt(xxi**2 + zxi**2)
- nx = + zxi / jacobian1D
- nz = - xxi / jacobian1D
- if (AXISYM) then
- if (is_on_the_axis(ispec_acoustic)) then
- weight = jacobian1D * wxglj(i) * r_xiplus1(i,j)
- else
- weight = jacobian1D * wxgll(i) * coord(1,ibool(i,j,ispec_acoustic))
- endif
- else
- weight = jacobian1D * wxgll(i)
- endif
- else if(iedge_acoustic == ILEFT)then
- xgamma = - xiz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- zgamma = + xix(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- jacobian1D = sqrt(xgamma**2 + zgamma**2)
- nx = - zgamma / jacobian1D
- nz = + xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
- else if(iedge_acoustic == IRIGHT)then
- xgamma = - xiz(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- zgamma = + xix(i,j,ispec_acoustic) * jacobian(i,j,ispec_acoustic)
- jacobian1D = sqrt(xgamma**2 + zgamma**2)
- nx = + zgamma / jacobian1D
- nz = - xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
- endif
-
- accel_elastic(1,iglob) = accel_elastic(1,iglob) + weight*nx*pressure
- accel_elastic(3,iglob) = accel_elastic(3,iglob) + weight*nz*pressure
-
- if(SIMULATION_TYPE == 3) then
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + weight*nx*b_pressure
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + weight*nz*b_pressure
- endif !if(SIMULATION_TYPE == 3) then
-
- enddo
-
- enddo
+ call compute_coupling_viscoelastic_ac(SIMULATION_TYPE,nspec,nglob_elastic,nglob_acoustic,num_fluid_solid_edges,&
+ ibool,wxgll,wzgll,xix,xiz,gammax,gammaz,jacobian,ivalue,jvalue,ivalue_inverse,jvalue_inverse,&
+ potential_acoustic,potential_acoustic_old,potential_dot_acoustic,potential_dot_dot_acoustic,&
+ b_potential_dot_dot_acoustic,accel_elastic,b_accel_elastic,fluid_solid_acoustic_ispec, &
+ fluid_solid_acoustic_iedge,fluid_solid_elastic_ispec,fluid_solid_elastic_iedge,&
+ potential_acoustic_adj_coupling,AXISYM,nglob,coord,is_on_the_axis,xiglj,wxglj, &
+ PML_BOUNDARY_CONDITIONS,nspec_PML,K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,&
+ alpha_z_store,is_PML,spec_to_PML,region_CPML,rmemory_sfb_potential_ddot_acoustic,timeval,deltat)
endif
More information about the CIG-COMMITS
mailing list