[cig-commits] r9215 - seismo/2D/SPECFEM2D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Feb 3 14:37:15 PST 2008
Author: dkomati1
Date: 2008-02-03 14:37:14 -0800 (Sun, 03 Feb 2008)
New Revision: 9215
Modified:
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
Changed Chi in the fluid from a displacement potential to a potential of (density * displacement) in order to have physical continuity of the potential everywhere including at a fluid-fluid discontinuity and be able to handle heterogeneous acoustic media
with first-order discontinuities (i.e., layers).
Also added test for negative maximum value to detect if the simulation became
unstable, because can happen with some compilers.
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2008-02-03 16:23:04 UTC (rev 9214)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2008-02-03 22:37:14 UTC (rev 9215)
@@ -190,23 +190,23 @@
! for the definition of potential energy in an acoustic fluid, see for instance
! equation (23) of M. Maess et al., Journal of Sound and Vibration 296 (2006) 264-276
-! in case of an acoustic medium, a displacement potential Chi is used as in Chaljub and Valette,
+! in case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette,
! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential
! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002).
! This permits acoustic-elastic coupling based on a non-iterative time scheme.
-! Displacement is then: u = grad(Chi)
-! Velocity is then: v = grad(Chi_dot) (Chi_dot being the time derivative of Chi)
-! and pressure is: p = - rho * Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi).
+! Displacement is then: u = grad(Chi) / rho
+! Velocity is then: v = grad(Chi_dot) / rho (Chi_dot being the time derivative of Chi)
+! and pressure is: p = - Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi).
! compute pressure in this element
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ numat,kmato,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute velocity vector field in this element
call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
! get density of current spectral element
lambdal_relaxed = elastcoef(1,kmato(ispec))
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2008-02-03 16:23:04 UTC (rev 9214)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2008-02-03 22:37:14 UTC (rev 9215)
@@ -45,7 +45,7 @@
assign_external_model,initialfield,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
- vpext,source_time_function,hprime_xx,hprimewgll_xx, &
+ vpext,rhoext,source_time_function,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
@@ -74,7 +74,7 @@
double precision, dimension(numat) :: density
double precision, dimension(4,numat) :: elastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
- double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
real(kind=CUSTOM_REAL), dimension(NSTEP) :: source_time_function
! derivatives of Lagrange polynomials
@@ -106,7 +106,7 @@
real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl,jacobianl
! material properties of the elastic medium
- real(kind=CUSTOM_REAL) :: mul_relaxed,lambdal_relaxed,kappal,cpl
+ real(kind=CUSTOM_REAL) :: mul_relaxed,lambdal_relaxed,kappal,cpl,rhol
! loop over spectral elements
do ispec_inner_outer = 1,nspec_inner_outer
@@ -118,6 +118,8 @@
!---
if(.not. elastic(ispec)) then
+ rhol = density(kmato(ispec))
+
! first double loop over GLL points to compute and store gradients
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -144,10 +146,13 @@
jacobianl = jacobian(i,j,ispec)
+! if external density model
+ if(assign_external_model) rhol = rhoext(i,j,ispec)
+
! for acoustic medium
! also add GLL integration weights
- tempx1(i,j) = wzgll(j)*jacobianl*(xixl*dux_dxl + xizl*dux_dzl)
- tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*dux_dxl + gammazl*dux_dzl)
+ tempx1(i,j) = wzgll(j)*jacobianl*(xixl*dux_dxl + xizl*dux_dzl) / rhol
+ tempx2(i,j) = wxgll(i)*jacobianl*(gammaxl*dux_dxl + gammazl*dux_dzl) / rhol
enddo
enddo
@@ -190,9 +195,9 @@
lambdal_relaxed = elastcoef(1,kmato(ispec))
mul_relaxed = elastcoef(2,kmato(ispec))
kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
- cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_relaxed/3._CUSTOM_REAL)/density(kmato(ispec)))
+ rhol = density(kmato(ispec))
+ cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_relaxed/3._CUSTOM_REAL)/rhol)
-
!--- left absorbing boundary
if(codeabs(ILEFT,ispecabs)) then
@@ -206,7 +211,10 @@
iglob = ibool(i,j,ispec)
! external velocity model
- if(assign_external_model) cpl = vpext(i,j,ispec)
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
@@ -216,7 +224,7 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec)) &
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
enddo
@@ -235,7 +243,10 @@
iglob = ibool(i,j,ispec)
! external velocity model
- if(assign_external_model) cpl = vpext(i,j,ispec)
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
@@ -245,7 +256,7 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec)) &
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
enddo
@@ -268,7 +279,10 @@
iglob = ibool(i,j,ispec)
! external velocity model
- if(assign_external_model) cpl = vpext(i,j,ispec)
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
@@ -278,7 +292,7 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec)) &
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
enddo
@@ -301,7 +315,10 @@
iglob = ibool(i,j,ispec)
! external velocity model
- if(assign_external_model) cpl = vpext(i,j,ispec)
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
@@ -311,7 +328,7 @@
! Sommerfeld condition if acoustic
if(.not. elastic(ispec)) &
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) - potential_dot_acoustic(iglob)*weight/cpl/rhol
enddo
@@ -328,14 +345,15 @@
if (is_proc_source == 1 ) then
! collocated force
! beware, for acoustic medium, source is a pressure source
+! the sign is negative because pressure p = - Chi_dot_dot therefore we need
+! to add minus the source to Chi_dot_dot to get plus the source in pressure
if(source_type == 1) then
if(.not. elastic(ispec_selected_source)) then
- potential_dot_dot_acoustic(iglob_source) = potential_dot_dot_acoustic(iglob_source) + source_time_function(it)
+ potential_dot_dot_acoustic(iglob_source) = potential_dot_dot_acoustic(iglob_source) - source_time_function(it)
endif
! moment tensor
else if(source_type == 2) then
-
if(.not. elastic(ispec_selected_source)) then
call exit_MPI('cannot have moment tensor source in acoustic element')
endif
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2008-02-03 16:23:04 UTC (rev 9214)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2008-02-03 22:37:14 UTC (rev 9215)
@@ -42,7 +42,7 @@
subroutine compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,elastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,e1,e11, &
+ numat,kmato,elastcoef,vpext,vsext,rhoext,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute pressure in acoustic elements and in elastic elements
@@ -56,7 +56,6 @@
integer, dimension(nspec) :: kmato
integer, dimension(NGLLX,NGLLX,nspec) :: ibool
- double precision, dimension(numat) :: density
double precision, dimension(4,numat) :: elastcoef
double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,vsext,rhoext
@@ -89,7 +88,7 @@
! compute pressure in this element
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ numat,kmato,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
! use vector_field_display as temporary storage, store pressure in its second component
@@ -110,7 +109,7 @@
subroutine compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ numat,kmato,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute pressure in acoustic elements and in elastic elements
@@ -124,7 +123,6 @@
integer, dimension(nspec) :: kmato
integer, dimension(NGLLX,NGLLX,nspec) :: ibool
- double precision, dimension(numat) :: density
double precision, dimension(4,numat) :: elastcoef
double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,vsext,rhoext
@@ -161,7 +159,6 @@
real(kind=CUSTOM_REAL) :: sigma_xx,sigma_zz
! material properties of the elastic medium
- integer :: material
real(kind=CUSTOM_REAL) :: mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed,denst
real(kind=CUSTOM_REAL) :: mul_unrelaxed,lambdal_unrelaxed,lambdalplus2mul_unrelaxed,cpl,csl
@@ -286,7 +283,7 @@
enddo
enddo
-! pressure = - rho * Chi_dot_dot if acoustic element
+! pressure = - Chi_dot_dot if acoustic element
else
do j = 1,NGLLZ
@@ -294,12 +291,8 @@
iglob = ibool(i,j,ispec)
- material = kmato(ispec)
- denst = density(material)
- if(assign_external_model) denst = rhoext(i,j,ispec)
-
! store pressure
- pressure_element(i,j) = - denst * potential_dot_dot_acoustic(iglob)
+ pressure_element(i,j) = - potential_dot_dot_acoustic(iglob)
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2008-02-03 16:23:04 UTC (rev 9214)
+++ seismo/2D/SPECFEM2D/trunk/compute_vector_field.f90 2008-02-03 22:37:14 UTC (rev 9215)
@@ -41,7 +41,7 @@
!========================================================================
subroutine compute_vector_whole_medium(potential_acoustic,veloc_elastic,elastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
! compute Grad(potential) in acoustic elements
! and combine with existing velocity vector field in elastic elements
@@ -50,8 +50,16 @@
include "constants.h"
- integer nspec,npoin
+ integer nspec,npoin,numat
+ logical :: assign_external_model
+
+ integer, dimension(nspec) :: kmato
+
+ double precision, dimension(NGLLX,NGLLX,nspec) :: rhoext
+
+ double precision, dimension(numat) :: density
+
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
@@ -76,7 +84,7 @@
! compute vector field in this element
call compute_vector_one_element(vector_field_element,potential_acoustic,veloc_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
! store the result
do j = 1,NGLLZ
@@ -95,7 +103,7 @@
!
subroutine compute_vector_one_element(vector_field_element,potential_acoustic,veloc_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
! compute Grad(potential) if acoustic element or copy existing vector if elastic element
@@ -103,8 +111,16 @@
include "constants.h"
- integer nspec,npoin,ispec
+ integer nspec,npoin,ispec,numat
+ logical :: assign_external_model
+
+ integer, dimension(nspec) :: kmato
+
+ double precision, dimension(NGLLX,NGLLX,nspec) :: rhoext
+
+ double precision, dimension(numat) :: density
+
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
@@ -130,6 +146,9 @@
! jacobian
real(kind=CUSTOM_REAL) xixl,xizl,gammaxl,gammazl
+! material properties of the elastic medium
+ real(kind=CUSTOM_REAL) :: rhol
+
! simple copy of existing vector if elastic element
if(elastic(ispec)) then
@@ -142,8 +161,11 @@
enddo
! compute gradient of potential to calculate vector if acoustic element
+! we then need to divide by density because the potential is a potential of (density * displacement)
else
+ rhol = density(kmato(ispec))
+
! double loop over GLL points to compute and store gradients
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -169,9 +191,11 @@
gammaxl = gammax(i,j,ispec)
gammazl = gammaz(i,j,ispec)
+ if(assign_external_model) rhol = rhoext(i,j,ispec)
+
! derivatives of potential
- vector_field_element(1,i,j) = tempx1l*xixl + tempx2l*gammaxl
- vector_field_element(2,i,j) = tempx1l*xizl + tempx2l*gammazl
+ vector_field_element(1,i,j) = (tempx1l*xixl + tempx2l*gammaxl) / rhol
+ vector_field_element(2,i,j) = (tempx1l*xizl + tempx2l*gammazl) / rhol
enddo
enddo
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-02-03 16:23:04 UTC (rev 9214)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-02-03 22:37:14 UTC (rev 9215)
@@ -109,14 +109,21 @@
! Institut de Physique du Globe de Paris, France
!
-! in case of an acoustic medium, a displacement potential Chi is used as in Chaljub and Valette,
+! in case of an acoustic medium, a potential Chi of (density * displacement) is used as in Chaljub and Valette,
! Geophysical Journal International, vol. 158, p. 131-141 (2004) and *NOT* a velocity potential
! as in Komatitsch and Tromp, Geophysical Journal International, vol. 150, p. 303-318 (2002).
! This permits acoustic-elastic coupling based on a non-iterative time scheme.
-! Displacement is then: u = grad(Chi)
-! Velocity is then: v = grad(Chi_dot) (Chi_dot being the time derivative of Chi)
-! and pressure is: p = - rho * Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi).
+! Displacement is then: u = grad(Chi) / rho
+! Velocity is then: v = grad(Chi_dot) / rho (Chi_dot being the time derivative of Chi)
+! and pressure is: p = - Chi_dot_dot (Chi_dot_dot being the time second derivative of Chi).
! The source in an acoustic element is a pressure source.
+! First-order acoustic-acoustic discontinuities are also handled automatically
+! because pressure is continuous at such an interface, therefore Chi_dot_dot
+! is continuous, therefore Chi is also continuous, which is consistent with
+! the spectral-element basis functions and with the assembling process.
+! This is the reason why a simple displacement potential u = grad(Chi) would
+! not work because it would be discontinuous at such an interface and would
+! therefore not be consistent with the basis functions.
program specfem2D
@@ -174,7 +181,7 @@
double precision :: xixl,xizl,gammaxl,gammazl,jacobianl
! material properties of the elastic medium
- double precision :: mul_relaxed,lambdal_relaxed,cpsquare
+ double precision :: mul_relaxed,lambdal_relaxed,kappal
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_elastic,veloc_elastic,displ_elastic
double precision, dimension(:,:), allocatable :: coord, flagrange,xinterp,zinterp,Uxinterp,Uzinterp,elastcoef,vector_field_display
@@ -634,17 +641,9 @@
!
!---- read interfaces data
!
- print *, 'read the interfaces', myrank
read(IIN,"(a80)") datlin
read(IIN,*) ninterface, max_interface_size
- if ( ninterface == 0 ) then
- !allocate(my_neighbours(1))
- !allocate(my_nelmnts_neighbours(1))
- !allocate(my_interfaces(4,1,1))
- !allocate(ibool_interfaces(NGLLX*1,1,1))
- !allocate(nibool_interfaces(1,1))
-
- else
+ if ( ninterface > 0 ) then
allocate(my_neighbours(ninterface))
allocate(my_nelmnts_neighbours(ninterface))
allocate(my_interfaces(4,max_interface_size,ninterface))
@@ -661,13 +660,10 @@
read(IIN,*) my_interfaces(1,ie,num_interface), my_interfaces(2,ie,num_interface), &
my_interfaces(3,ie,num_interface), my_interfaces(4,ie,num_interface)
- end do
- end do
- print *, 'end read the interfaces', myrank
+ enddo
+ enddo
+ endif
- end if
-
-
!
!---- read absorbing boundary data
!
@@ -1095,18 +1091,20 @@
! if external density model
if(assign_external_model) then
rhol = rhoext(i,j,ispec)
- cpsquare = vpext(i,j,ispec)**2
+ kappal = rhol * vpext(i,j,ispec)**2
else
rhol = density(kmato(ispec))
lambdal_relaxed = elastcoef(1,kmato(ispec))
mul_relaxed = elastcoef(2,kmato(ispec))
- cpsquare = (lambdal_relaxed + 2.d0*mul_relaxed) / rhol
+ kappal = lambdal_relaxed + 2.d0*mul_relaxed
endif
-! for acoustic medium
if(elastic(ispec)) then
+! for elastic medium
rmass_inverse_elastic(iglob) = rmass_inverse_elastic(iglob) + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
else
- rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / cpsquare
+! for acoustic medium
+ rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) + &
+ wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
endif
enddo
enddo
@@ -1958,8 +1956,7 @@
if(assign_external_model) then
vp_display(ibool(i,j,ispec)) = vpext(i,j,ispec)
else
- cpsquare = (lambdal_relaxed + 2.d0*mul_relaxed) / rhol
- vp_display(ibool(i,j,ispec)) = sqrt(cpsquare)
+ vp_display(ibool(i,j,ispec)) = sqrt((lambdal_relaxed + 2.d0*mul_relaxed) / rhol)
endif
enddo
enddo
@@ -2054,12 +2051,11 @@
assign_external_model,initialfield,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
- vpext,source_time_function,hprime_xx,hprimewgll_xx, &
+ vpext,rhoext,source_time_function,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
- nspec_outer, ispec_outer_to_glob, .true. &
- )
+ nspec_outer, ispec_outer_to_glob, .true.)
endif ! end of test if any acoustic element
@@ -2138,8 +2134,7 @@
max_interface_size, max_ibool_interfaces_size_ac,&
ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
tab_requests_send_recv_acoustic, &
- buffer_send_faces_vector_ac &
- )
+ buffer_send_faces_vector_ac)
endif
#endif
@@ -2150,12 +2145,11 @@
assign_external_model,initialfield,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
- vpext,source_time_function,hprime_xx,hprimewgll_xx, &
+ vpext,rhoext,source_time_function,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
- nspec_inner, ispec_inner_to_glob, .false. &
- )
+ nspec_inner, ispec_inner_to_glob, .false.)
endif
! assembling potential_dot_dot for acoustic elements (receive)
@@ -2167,8 +2161,7 @@
max_interface_size, max_ibool_interfaces_size_ac,&
ibool_interfaces_acoustic, nibool_interfaces_acoustic, &
tab_requests_send_recv_acoustic, &
- buffer_recv_faces_vector_ac &
- )
+ buffer_recv_faces_vector_ac)
endif
#endif
@@ -2233,15 +2226,8 @@
j = jvalue_inverse(ipoin1D,iedge_acoustic)
iglob = ibool(i,j,ispec_acoustic)
-! get density of the fluid, depending if external density model
- if(assign_external_model) then
- rhol = rhoext(i,j,ispec_acoustic)
- else
- rhol = density(kmato(ispec_acoustic))
- endif
-
! compute pressure on the fluid/solid edge
- pressure = - rhol * potential_dot_dot_acoustic(iglob)
+ pressure = - potential_dot_dot_acoustic(iglob)
! get point values for the elastic side
i = ivalue(ipoin1D,iedge_elastic)
@@ -2368,7 +2354,10 @@
write(IOUT,*) 'Max norm of vector field in solid = ',displnorm_all_glob
endif
! check stability of the code in solid, exit if unstable
- if(displnorm_all_glob > STABILITY_THRESHOLD) call exit_MPI('code became unstable and blew up in solid')
+! negative values can occur with some compilers when the unstable value is greater
+! than the greatest possible floating-point number of the machine
+ if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
+ call exit_MPI('code became unstable and blew up in solid')
endif
if(any_acoustic_glob) then
@@ -2385,7 +2374,10 @@
write(IOUT,*) 'Max absolute value of scalar field in fluid = ',displnorm_all_glob
endif
! check stability of the code in fluid, exit if unstable
- if(displnorm_all_glob > STABILITY_THRESHOLD) call exit_MPI('code became unstable and blew up in fluid')
+! negative values can occur with some compilers when the unstable value is greater
+! than the greatest possible floating-point number of the machine
+ if(displnorm_all_glob > STABILITY_THRESHOLD .or. displnorm_all_glob < 0) &
+ call exit_MPI('code became unstable and blew up in fluid')
endif
if ( myrank == 0 ) then
write(IOUT,*)
@@ -2404,7 +2396,7 @@
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
+ numat,kmato,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
else if(.not. elastic(ispec)) then
@@ -2412,13 +2404,13 @@
! for acoustic medium, compute vector field from gradient of potential for seismograms
if(seismotype == 1) then
call compute_vector_one_element(vector_field_element,potential_acoustic,displ_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
else if(seismotype == 2) then
call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
else if(seismotype == 3) then
call compute_vector_one_element(vector_field_element,potential_dot_dot_acoustic,accel_elastic,elastic, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
endif
endif
@@ -2500,7 +2492,7 @@
endif
call compute_vector_whole_medium(potential_acoustic,displ_elastic,elastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
it,deltat,coorg,xinterp,zinterp,shape2D_display, &
@@ -2521,7 +2513,7 @@
endif
call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,elastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
it,deltat,coorg,xinterp,zinterp,shape2D_display, &
@@ -2542,7 +2534,7 @@
endif
call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,elastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
call plotpost(vector_field_display,coord,vpext,x_source,z_source,st_xval,st_zval, &
it,deltat,coorg,xinterp,zinterp,shape2D_display, &
@@ -2588,7 +2580,7 @@
endif
call compute_vector_whole_medium(potential_acoustic,displ_elastic,elastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
else if(imagetype == 2) then
@@ -2597,7 +2589,7 @@
endif
call compute_vector_whole_medium(potential_dot_acoustic,veloc_elastic,elastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
else if(imagetype == 3) then
@@ -2606,7 +2598,7 @@
endif
call compute_vector_whole_medium(potential_dot_dot_acoustic,accel_elastic,elastic,vector_field_display, &
- xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin)
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,numat,kmato,density,rhoext,assign_external_model)
else if(imagetype == 4) then
@@ -2616,7 +2608,7 @@
call compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,elastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,elastcoef,vpext,vsext,rhoext,e1,e11, &
+ numat,kmato,elastcoef,vpext,vsext,rhoext,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
else
More information about the cig-commits
mailing list