[cig-commits] [commit] devel: added parameter READ_VELOCITIES_AT_f0 and implemented a shift of velocities read from the Par_file if needed when viscoelasticity is on (38c73cd)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed May 21 06:30:14 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branch : devel
Link : https://github.com/geodynamics/specfem2d/compare/967f9b955f19c24208e6a40cf9c157d9ff8b050b...34914ddff686f70de49968a592d69f78c08a30bb
>---------------------------------------------------------------
commit 38c73cd1492a80b13303784999a7955a7a11a5f2
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Wed May 21 15:28:10 2014 +0200
added parameter READ_VELOCITIES_AT_f0 and implemented a shift of velocities read from the Par_file if needed when viscoelasticity is on
>---------------------------------------------------------------
38c73cd1492a80b13303784999a7955a7a11a5f2
DATA/Par_file | 10 +-
EXAMPLES | 2 +-
...s_low_frequency_unrelaxed_is_high_frequency.txt | 1 +
...ion_reference_Specfem2D_fixed_by_Zhinan_Xie.pdf | Bin 39765 -> 42292 bytes
setup/constants.h.in | 3 +
src/meshfem2D/read_parameter_file.F90 | 4 +
src/meshfem2D/save_databases.f90 | 2 +-
src/specfem2D/acoustic_forcing_boundary.f90 | 0
src/specfem2D/attenuation_model.f90 | 105 ++++++++++++++++-----
src/specfem2D/read_databases.F90 | 5 +-
src/specfem2D/read_external_model.f90 | 35 +++++--
src/specfem2D/specfem2D.F90 | 99 ++++++++++++-------
12 files changed, 194 insertions(+), 72 deletions(-)
diff --git a/DATA/Par_file b/DATA/Par_file
index 5da10f1..f6bcdde 100644
--- a/DATA/Par_file
+++ b/DATA/Par_file
@@ -37,9 +37,15 @@ ACOUSTIC_FORCING = .false. # acoustic forcing of an acoust
NSOURCES = 1 # number of sources (source info read from DATA/CMTSOLUTION file)
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
-# constants for attenuation
+# for viscoelastic attenuation
N_SLS = 2 # number of standard linear solids for attenuation
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, otherwise it is f0 the dominant frequency of the source in the CMTSOLUTION file
+# shift (i.e. change) velocities read from the input file to take average physical dispersion into account,
+# i.e. if needed change the reference frequency at which these velocities are defined internally in the code:
+# by default, the seismic velocity values that are read at the end of this Par_file of the code are supposed to be the unrelaxed values,
+# i.e. the velocities at infinite frequency. We may want to change this and impose that the values read are those for a given frequency (here f0_attenuation).
+# (when we do this, velocities will then slightly decrease and waves will thus slightly slow down)
+READ_VELOCITIES_AT_f0 = .false.
# receiver set parameters for seismograms
seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl of displ 6=the fluid potential
@@ -71,7 +77,7 @@ cutsnaps = 1. # minimum amplitude kept in % f
#### for JPEG color images ####
output_color_image = .true. # output JPEG color image of the results every NSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
imagetype_JPEG = 2 # display 1=displ_Ux 2=displ_Uz 3=displ_norm 4=veloc_Vx 5=veloc_Vz 6=veloc_norm 7=accel_Ax 8=accel_Az 9=accel_norm 10=pressure
-factor_subsample_image = 1.0 # (double precision) factor to subsample color images output by the code (useful for very large models)
+factor_subsample_image = 1.0d0 # (double precision) factor to subsample color images output by the code (useful for very large models)
USE_CONSTANT_MAX_AMPLITUDE = .false. # by default the code normalizes each image independently to its maximum; use this option to use the global maximum below instead
CONSTANT_MAX_AMPLITUDE_TO_USE = 1.17d4 # constant maximum amplitude to use for all color images if the above USE_CONSTANT_MAX_AMPLITUDE option is true
POWER_DISPLAY_COLOR = 0.30d0 # non linear display to enhance small amplitudes in JPEG color images
diff --git a/EXAMPLES b/EXAMPLES
index 599ea1a..84746be 160000
--- a/EXAMPLES
+++ b/EXAMPLES
@@ -1 +1 @@
-Subproject commit 599ea1a4656dc2499687591f676dad14cf843f63
+Subproject commit 84746be287e8c81c94cb9db5198858003d9c6ebe
diff --git a/doc/for_attenuation_relaxed_is_low_frequency_unrelaxed_is_high_frequency.txt b/doc/for_attenuation_relaxed_is_low_frequency_unrelaxed_is_high_frequency.txt
new file mode 100644
index 0000000..157d58b
--- /dev/null
+++ b/doc/for_attenuation_relaxed_is_low_frequency_unrelaxed_is_high_frequency.txt
@@ -0,0 +1 @@
+for attenuation, the relaxed state is at low frequency and the unrelaxed state is at high frequency.
diff --git a/doc/old_problem_attenuation_reference_Specfem2D_fixed_by_Zhinan_Xie.pdf b/doc/old_problem_attenuation_reference_Specfem2D_fixed_by_Zhinan_Xie.pdf
index e1a50bf..abc2ccf 100644
Binary files a/doc/old_problem_attenuation_reference_Specfem2D_fixed_by_Zhinan_Xie.pdf and b/doc/old_problem_attenuation_reference_Specfem2D_fixed_by_Zhinan_Xie.pdf differ
diff --git a/setup/constants.h.in b/setup/constants.h.in
index e4f17d5..e3cacfe 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -121,6 +121,9 @@
! pi
double precision, parameter :: PI = 3.141592653589793d0
+! 2/3
+ double precision, parameter :: TWO_THIRDS = 2.d0/3.d0
+
! 4/3
double precision, parameter :: FOUR_THIRDS = 4.d0/3.d0
diff --git a/src/meshfem2D/read_parameter_file.F90 b/src/meshfem2D/read_parameter_file.F90
index 1f45acf..4b6f5fd 100644
--- a/src/meshfem2D/read_parameter_file.F90
+++ b/src/meshfem2D/read_parameter_file.F90
@@ -89,6 +89,7 @@ module parameter_file
! variables used for attenuation
integer :: N_SLS
double precision :: f0_attenuation
+ logical :: READ_VELOCITIES_AT_f0
integer :: seismotype
logical :: use_existing_STATIONS
@@ -275,6 +276,9 @@ contains
call read_value_double_precision_p(f0_attenuation, 'solver.f0_attenuation')
if(err_occurred() /= 0) stop 'error reading parameter 21 in Par_file'
+ call read_value_logical_p(READ_VELOCITIES_AT_f0, 'solver.READ_VELOCITIES_AT_f0')
+ if(err_occurred() /= 0) stop 'error reading parameter 21b in Par_file'
+
! read receiver line parameters
call read_value_integer_p(seismotype, 'solver.seismotype')
if(err_occurred() /= 0) stop 'error reading parameter 22 in Par_file'
diff --git a/src/meshfem2D/save_databases.f90 b/src/meshfem2D/save_databases.f90
index e9c9d64..262c9ed 100644
--- a/src/meshfem2D/save_databases.f90
+++ b/src/meshfem2D/save_databases.f90
@@ -241,7 +241,7 @@
enddo
write(15,*) 'attenuation'
- write(15,*) N_SLS, f0_attenuation
+ write(15,*) N_SLS, f0_attenuation, READ_VELOCITIES_AT_f0
write(15,*) 'Coordinates of macrobloc mesh (coorg):'
diff --git a/src/specfem2D/attenuation_model.f90 b/src/specfem2D/attenuation_model.f90
index 244440d..2ea6580 100644
--- a/src/specfem2D/attenuation_model.f90
+++ b/src/specfem2D/attenuation_model.f90
@@ -42,7 +42,7 @@
!========================================================================
subroutine attenuation_model(N_SLS,QKappa_attenuation,Qmu_attenuation,f0_attenuation, &
- inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2)
+ tau_epsilon_nu1_custom_real,tau_epsilon_nu2_custom_real,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2)
! define the attenuation constants
@@ -52,7 +52,8 @@
integer :: N_SLS
double precision :: QKappa_attenuation,Qmu_attenuation,f0_attenuation
- real(kind=CUSTOM_REAL), dimension(N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
+ real(kind=CUSTOM_REAL), dimension(N_SLS) :: tau_epsilon_nu1_custom_real,tau_epsilon_nu2_custom_real, &
+ inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
real(kind=CUSTOM_REAL) :: Mu_nu1,Mu_nu2
integer :: i_sls
@@ -113,29 +114,14 @@
! tau_epsilon_nu2(2) = 0.0033257d0
! tau_sigma_nu2(2) = 0.0030465d0
-! values for Paul Cristini for fluid-solid ocean acoustics simulations
-
-! for N_SLS = 2
-! frequency range: 1.500000 Hz - 18.000000 Hz
-! central frequency in log scale in Hz = 5.196152422706633
-! target constant attenuation factor Q = 136.4376068115
-! tau sigma evenly spaced in log frequency, do not depend on value of Q
-
-! tau_sigma_nu1(1) = 0.10610329539459699422d0
-! tau_sigma_nu1(2) = 0.00884194128288308401d0
-
-! tau_epsilon_nu1(1) = 0.10754721280605997191d0
-! tau_epsilon_nu1(2) = 0.00895488050110176612d0
-
-! tau_epsilon_nu2(1) = tau_epsilon_nu1(1)
-! tau_epsilon_nu2(2) = tau_epsilon_nu1(2)
-! tau_sigma_nu2(1) = tau_sigma_nu1(1)
-! tau_sigma_nu2(2) = tau_sigma_nu1(2)
-
!
!--- other constants computed from the parameters above, do not modify
!
if(CUSTOM_REAL == SIZE_REAL) then
+
+ tau_epsilon_nu1_custom_real(:) = sngl(tau_epsilon_nu1(:))
+ tau_epsilon_nu2_custom_real(:) = sngl(tau_epsilon_nu2(:))
+
inv_tau_sigma_nu1(:) = sngl(dble(ONE) / tau_sigma_nu1(:))
inv_tau_sigma_nu2(:) = sngl(dble(ONE) / tau_sigma_nu2(:))
@@ -149,7 +135,12 @@
Mu_nu1 = sngl(dble(Mu_nu1) - (dble(ONE) - tau_epsilon_nu1(i_sls)/tau_sigma_nu1(i_sls)))
Mu_nu2 = sngl(dble(Mu_nu2) - (dble(ONE) - tau_epsilon_nu2(i_sls)/tau_sigma_nu2(i_sls)))
enddo
+
else
+
+ tau_epsilon_nu1_custom_real(:) = tau_epsilon_nu1(:)
+ tau_epsilon_nu2_custom_real(:) = tau_epsilon_nu2(:)
+
inv_tau_sigma_nu1(:) = ONE / tau_sigma_nu1(:)
inv_tau_sigma_nu2(:) = ONE / tau_sigma_nu2(:)
@@ -163,7 +154,79 @@
Mu_nu1 = Mu_nu1 - (ONE - tau_epsilon_nu1(i_sls)/tau_sigma_nu1(i_sls))
Mu_nu2 = Mu_nu2 - (ONE - tau_epsilon_nu2(i_sls)/tau_sigma_nu2(i_sls))
enddo
+
endif
end subroutine attenuation_model
+!--------------------------------------------------------------------------------
+
+ subroutine shift_velocities_from_f0(vp,vs,rho,mu,lambda,f0_attenuation, &
+ tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1,inv_tau_sigma_nu2,N_SLS)
+
+! From Emmanuel Chaljub, CNRS Grenoble, France:
+
+! shift (i.e. change) velocities read from the input file to take average physical dispersion into account,
+! i.e. if needed change the reference frequency at which these velocities are defined internally in the code:
+
+! by default, the velocity values that are read in the Par_file of the code are supposed to be the unrelaxed values,
+! i.e. the velocities at infinite frequency.
+! We may want to change this and impose that the values read are those for a given frequency (here f0_attenuation).
+!
+! The unrelaxed values are then defined from the values read at frequency f0_attenuation as follows:
+!
+! mu_unrelaxed = mu (w_ref) / [ 1 - ( Sum_k ak/(1+(w_ref*tau_k)**2) )/( 1 + Sum_k ak ) ]
+! where ak = tau_epsilon_k/tau_sigma_k - 1
+! and the ak are the solutions of the linear system:
+!
+! Sum_k { [ w*tau_k*Q(w) - (w*tau_k)^2 ] / [ 1 + (w*tau_k)^2 ] } ak = 1
+! where tau_k = tau_epsilon_k
+
+ implicit none
+
+ include "constants.h"
+
+! arguments
+ integer, intent(in) :: N_SLS
+ double precision, intent(in) :: rho,f0_attenuation
+ double precision, intent(inout) :: vp,vs
+ double precision, intent(out) :: mu,lambda
+ real(kind=CUSTOM_REAL), dimension(N_SLS), intent(in) :: tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1,inv_tau_sigma_nu2
+
+! local variables
+ integer :: i_sls
+ double precision :: xtmp1_nu1,xtmp1_nu2,xtmp2_nu1,xtmp2_nu2,xtmp_ak_nu1,xtmp_ak_nu2,delta_mu,kappa,delta_kappa
+
+ mu = rho * vs*vs
+ lambda = rho * vp*vp - TWO * mu
+ kappa = lambda + TWO_THIRDS*mu
+
+ xtmp1_nu1 = ZERO
+ xtmp2_nu1 = ONE
+ xtmp1_nu2 = ZERO
+ xtmp2_nu2 = ONE
+
+ do i_sls = 1,N_SLS
+!! DK DK changed this to the pre-computed inverse xtmp_ak_nu2 = tau_epsilon_nu2(i_sls)/tau_sigma_nu2(i_sls) - ONE
+ xtmp_ak_nu2 = tau_epsilon_nu2(i_sls)*inv_tau_sigma_nu2(i_sls) - ONE
+ xtmp1_nu2 = xtmp1_nu2 + xtmp_ak_nu2/(ONE + (TWO * PI * f0_attenuation * tau_epsilon_nu2(i_sls))**2)
+ xtmp2_nu2 = xtmp2_nu2 + xtmp_ak_nu2
+
+!! DK DK changed this to the pre-computed inverse xtmp_ak_nu1 = tau_epsilon_nu1(i_sls)/tau_sigma_nu1(i_sls) - ONE
+ xtmp_ak_nu1 = tau_epsilon_nu1(i_sls)*inv_tau_sigma_nu1(i_sls) - ONE
+ xtmp1_nu1 = xtmp1_nu1 + xtmp_ak_nu1/(ONE + (TWO * PI * f0_attenuation * tau_epsilon_nu1(i_sls))**2)
+ xtmp2_nu1 = xtmp2_nu1 + xtmp_ak_nu1
+ enddo
+
+ delta_mu = ONE - xtmp1_nu2/xtmp2_nu2
+ mu = mu / delta_mu
+
+ delta_kappa = ONE - xtmp1_nu1/xtmp2_nu1
+ kappa = kappa / delta_kappa
+
+ lambda = kappa - TWO_THIRDS*mu
+ vp = dsqrt((lambda + TWO * mu) / rho)
+ vs = dsqrt(mu / rho)
+
+ end subroutine shift_velocities_from_f0
+
diff --git a/src/specfem2D/read_databases.F90 b/src/specfem2D/read_databases.F90
index 1323506..37f4b06 100644
--- a/src/specfem2D/read_databases.F90
+++ b/src/specfem2D/read_databases.F90
@@ -452,7 +452,7 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine read_databases_atten(N_SLS,f0_attenuation)
+ subroutine read_databases_atten(N_SLS,f0_attenuation,READ_VELOCITIES_AT_f0)
! reads attenuation information
@@ -461,12 +461,13 @@
integer :: N_SLS
double precision :: f0_attenuation
+ logical :: READ_VELOCITIES_AT_f0
! local parameters
character(len=80) :: datlin
read(IIN,"(a80)") datlin
- read(IIN,*) N_SLS, f0_attenuation
+ read(IIN,*) N_SLS, f0_attenuation, READ_VELOCITIES_AT_f0
if(N_SLS < 2) stop 'must have N_SLS >= 2 even if attenuation if off because it is used to assign some arrays'
end subroutine read_databases_atten
diff --git a/src/specfem2D/read_external_model.f90 b/src/specfem2D/read_external_model.f90
index 33c795b..699718b 100644
--- a/src/specfem2D/read_external_model.f90
+++ b/src/specfem2D/read_external_model.f90
@@ -44,25 +44,27 @@
subroutine read_external_model(any_acoustic,any_gravitoacoustic,any_elastic,any_poroelastic, &
acoustic,gravitoacoustic,elastic,poroelastic,anisotropic,nspec,nglob,N_SLS,ibool, &
- f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
+ f0_attenuation,READ_VELOCITIES_AT_f0,tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent, &
inv_tau_sigma_nu1,inv_tau_sigma_nu2,phi_nu1,phi_nu2,Mu_nu1,Mu_nu2,&
coord,kmato,rhoext,vpext,vsext,gravityext,Nsqext, &
QKappa_attenuationext,Qmu_attenuationext, &
- c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext,READ_EXTERNAL_SEP_FILE)
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext, &
+ READ_EXTERNAL_SEP_FILE,ATTENUATION_VISCOELASTIC_SOLID)
implicit none
include "constants.h"
integer :: nspec,nglob
double precision :: f0_attenuation
+ logical :: READ_VELOCITIES_AT_f0
! Mesh
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
double precision, dimension(NDIM,nglob) :: coord
! Material properties
- logical :: any_acoustic,any_gravitoacoustic,any_elastic,any_poroelastic,READ_EXTERNAL_SEP_FILE
+ logical :: any_acoustic,any_gravitoacoustic,any_elastic,any_poroelastic,READ_EXTERNAL_SEP_FILE,ATTENUATION_VISCOELASTIC_SOLID
integer, dimension(nspec) :: kmato
logical, dimension(nspec) :: acoustic,gravitoacoustic,elastic,poroelastic
double precision, dimension(NGLLX,NGLLZ,nspec) :: rhoext,vpext,vsext,gravityext,Nsqext
@@ -70,7 +72,7 @@
! for attenuation
integer :: N_SLS
real(kind=CUSTOM_REAL) :: Mu_nu1_sent,Mu_nu2_sent
- real(kind=CUSTOM_REAL), dimension(N_SLS) :: inv_tau_sigma_nu1_sent,phi_nu1_sent, &
+ real(kind=CUSTOM_REAL), dimension(N_SLS) :: tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
inv_tau_sigma_nu2_sent,phi_nu2_sent
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: inv_tau_sigma_nu1,phi_nu1, &
inv_tau_sigma_nu2,phi_nu2
@@ -86,6 +88,8 @@
double precision :: previous_vsext
double precision :: tmp1, tmp2,tmp3
+ double precision :: mu_dummy,lambda_dummy
+
if(READ_EXTERNAL_SEP_FILE) then
write(IOUT,*)
write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
@@ -141,6 +145,7 @@
enddo
enddo
enddo
+
endif
! initializes
@@ -156,13 +161,15 @@
poroelastic(:) = .false.
do ispec = 1,nspec
+
previous_vsext = -1.d0
+
do j = 1,NGLLZ
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
if(.not. (i == 1 .and. j == 1) .and. &
((vsext(i,j,ispec) >= TINYVAL .and. previous_vsext < TINYVAL) .or. &
- (vsext(i,j,ispec) < TINYVAL .and. previous_vsext >= TINYVAL))) &
+ (vsext(i,j,ispec) < TINYVAL .and. previous_vsext >= TINYVAL))) &
call exit_MPI('external velocity model cannot be both fluid and solid inside the same spectral element')
if(c11ext(i,j,ispec) > TINYVAL .or. c13ext(i,j,ispec) > TINYVAL .or. c15ext(i,j,ispec) > TINYVAL .or. &
@@ -173,13 +180,13 @@
any_elastic = .true.
QKappa_attenuationext(i,j,ispec) = 10.d0
Qmu_attenuationext(i,j,ispec) = 10.d0
- else if((vsext(i,j,ispec) < TINYVAL).and.(gravityext(i,j,ispec) < TINYVAL)) then
+ else if((vsext(i,j,ispec) < TINYVAL) .and. (gravityext(i,j,ispec) < TINYVAL)) then
elastic(ispec) = .false.
poroelastic(ispec) = .false.
gravitoacoustic(ispec)=.false.
acoustic(ispec)=.true.
any_acoustic = .true.
- else if((vsext(i,j,ispec) < TINYVAL).and.(gravityext(i,j,ispec) >= TINYVAL)) then
+ else if((vsext(i,j,ispec) < TINYVAL) .and. (gravityext(i,j,ispec) >= TINYVAL)) then
elastic(ispec) = .false.
poroelastic(ispec) = .false.
acoustic(ispec)=.false.
@@ -192,17 +199,27 @@
endif
call attenuation_model(N_SLS,QKappa_attenuationext(i,j,ispec),Qmu_attenuationext(i,j,ispec), &
- f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
- inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent)
+ f0_attenuation,tau_epsilon_nu1,tau_epsilon_nu2, &
+ inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent)
inv_tau_sigma_nu1(i,j,ispec,:) = inv_tau_sigma_nu1_sent(:)
phi_nu1(i,j,ispec,:) = phi_nu1_sent(:)
inv_tau_sigma_nu2(i,j,ispec,:) = inv_tau_sigma_nu2_sent(:)
phi_nu2(i,j,ispec,:) = phi_nu2_sent(:)
Mu_nu1(i,j,ispec) = Mu_nu1_sent
Mu_nu2(i,j,ispec) = Mu_nu2_sent
+
+ if(ATTENUATION_VISCOELASTIC_SOLID .and. READ_VELOCITIES_AT_F0) then
+ if(anisotropic(ispec) .or. poroelastic(ispec) .or. gravitoacoustic(ispec)) stop &
+ 'READ_VELOCITIES_AT_F0 only implemented for non anisotropic, non poroelastic, non gravitoacoustic materials for now'
+ call shift_velocities_from_f0(vpext(i,j,ispec),vsext(i,j,ispec),rhoext(i,j,ispec),mu_dummy,lambda_dummy, &
+ f0_attenuation,tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1_sent,inv_tau_sigma_nu2_sent,N_SLS)
+ endif
+
previous_vsext = vsext(i,j,ispec)
+
enddo
enddo
enddo
end subroutine read_external_model
+
diff --git a/src/specfem2D/specfem2D.F90 b/src/specfem2D/specfem2D.F90
index 95d5524..d69f6ea 100644
--- a/src/specfem2D/specfem2D.F90
+++ b/src/specfem2D/specfem2D.F90
@@ -512,6 +512,7 @@
! poroelastic and elastic coefficients
double precision, dimension(:,:,:), allocatable :: poroelastcoef
+ logical, dimension(:), allocatable :: already_shifted_velocity
! anisotropy parameters
logical :: all_anisotropic
@@ -623,6 +624,7 @@
double precision, dimension(:), allocatable :: QKappa_attenuation
double precision, dimension(:), allocatable :: Qmu_attenuation
double precision :: f0_attenuation
+ logical :: READ_VELOCITIES_AT_f0
integer nspec_allocate
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
@@ -630,7 +632,8 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_initial_rk,e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: e1_force_rk,e11_force_rk,e13_force_rk
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: tau_epsilon_nu1,tau_epsilon_nu2, &
+ inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent
real(kind=CUSTOM_REAL), dimension(:,:,:) , allocatable :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL) :: Mu_nu1_sent,Mu_nu2_sent
@@ -1078,8 +1081,6 @@
logical :: anyabs_glob
integer :: nspec_PML
-!! logical :: backward_simulation for seprate adjoint simulation from backward simulation
-
! acoustic (and gravitoacoustic) forcing parameters
logical, dimension(:,:), allocatable :: codeacforcing
integer, dimension(:), allocatable :: typeacforcing
@@ -1091,6 +1092,9 @@
integer :: nspec_left_acforcing,nspec_right_acforcing,nspec_bottom_acforcing,nspec_top_acforcing
integer, dimension(:), allocatable :: ib_left_acforcing,ib_right_acforcing,ib_bottom_acforcing,ib_top_acforcing
+! for shifting of velocities if needed in the case of viscoelasticity
+ double precision :: vp,vs,rho,mu,lambda
+
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
@@ -1180,7 +1184,7 @@
endif
!---- read attenuation information
- call read_databases_atten(N_SLS,f0_attenuation)
+ call read_databases_atten(N_SLS,f0_attenuation,READ_VELOCITIES_AT_f0)
! if source is not a Dirac or Heavyside then f0_attenuation is f0 of the first source
if(.not. (time_function_type(1) == 4 .or. time_function_type(1) == 5)) then
@@ -1227,6 +1231,7 @@
allocate(tortuosity(numat))
allocate(permeability(3,numat))
allocate(poroelastcoef(4,3,numat))
+ allocate(already_shifted_velocity(numat))
allocate(QKappa_attenuation(numat))
allocate(Qmu_attenuation(numat))
allocate(kmato(nspec))
@@ -1241,11 +1246,15 @@
allocate(inv_tau_sigma_nu2(NGLLX,NGLLZ,nspec,N_SLS))
allocate(phi_nu1(NGLLX,NGLLZ,nspec,N_SLS))
allocate(phi_nu2(NGLLX,NGLLZ,nspec,N_SLS))
+ allocate(tau_epsilon_nu1(N_SLS))
+ allocate(tau_epsilon_nu2(N_SLS))
allocate(inv_tau_sigma_nu1_sent(N_SLS))
allocate(inv_tau_sigma_nu2_sent(N_SLS))
allocate(phi_nu1_sent(N_SLS))
allocate(phi_nu2_sent(N_SLS))
+ already_shifted_velocity(:) = .false.
+
!
!---- read the material properties
!
@@ -1419,7 +1428,7 @@
!! DK DK if needed in the future, here the quality factor could be different for each point
do ispec = 1,nspec
call attenuation_model(N_SLS,QKappa_attenuation(kmato(ispec)),Qmu_attenuation(kmato(ispec)), &
- f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
+ f0_attenuation,tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent)
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -1431,6 +1440,27 @@
Mu_nu2(i,j,ispec) = Mu_nu2_sent
enddo
enddo
+
+ if(ATTENUATION_VISCOELASTIC_SOLID .and. READ_VELOCITIES_AT_F0 .and. .not. assign_external_model) then
+ if(anisotropic(ispec) .or. poroelastic(ispec) .or. gravitoacoustic(ispec)) &
+ stop 'READ_VELOCITIES_AT_F0 only implemented for non anisotropic, non poroelastic, non gravitoacoustic materials for now'
+ n = kmato(ispec)
+ if(.not. already_shifted_velocity(n)) then
+ rho = density(1,n)
+ lambda = poroelastcoef(1,1,n)
+ mu = poroelastcoef(2,1,n)
+ vp = dsqrt((lambda + TWO * mu) / rho)
+ vs = dsqrt(mu / rho)
+ call shift_velocities_from_f0(vp,vs,rho,mu,lambda,f0_attenuation, &
+ tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1_sent,inv_tau_sigma_nu2_sent,N_SLS)
+ poroelastcoef(1,1,n) = lambda
+ poroelastcoef(2,1,n) = mu
+ poroelastcoef(3,1,n) = lambda + TWO*mu
+ already_shifted_velocity(n) = .true.
+ endif
+
+ endif
+
enddo
! allocate memory variables for viscous attenuation (poroelastic media)
@@ -2212,12 +2242,13 @@
if(myrank == 0) write(IOUT,*) 'Assigning an external velocity and density model...'
call read_external_model(any_acoustic,any_gravitoacoustic,any_elastic,any_poroelastic, &
acoustic,gravitoacoustic,elastic,poroelastic,anisotropic,nspec,nglob,N_SLS,ibool, &
- f0_attenuation,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
+ f0_attenuation,READ_VELOCITIES_AT_f0,tau_epsilon_nu1,tau_epsilon_nu2,inv_tau_sigma_nu1_sent,phi_nu1_sent, &
inv_tau_sigma_nu2_sent,phi_nu2_sent,Mu_nu1_sent,Mu_nu2_sent, &
inv_tau_sigma_nu1,inv_tau_sigma_nu2,phi_nu1,phi_nu2,Mu_nu1,Mu_nu2,&
coord,kmato,rhoext,vpext,vsext,gravityext,Nsqext, &
QKappa_attenuationext,Qmu_attenuationext, &
- c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext,READ_EXTERNAL_SEP_FILE)
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,c12ext,c23ext,c25ext, &
+ READ_EXTERNAL_SEP_FILE,ATTENUATION_VISCOELASTIC_SOLID)
endif
!
@@ -2736,21 +2767,21 @@
allocate(rmass_w_inverse_poroelastic(nglob_poroelastic))
if(time_stepping_scheme == 2)then
- allocate(displs_poroelastic_LDDRK(NDIM,nglob_poroelastic))
- allocate(velocs_poroelastic_LDDRK(NDIM,nglob_poroelastic))
- allocate(displw_poroelastic_LDDRK(NDIM,nglob_poroelastic))
- allocate(velocw_poroelastic_LDDRK(NDIM,nglob_poroelastic))
+ allocate(displs_poroelastic_LDDRK(NDIM,nglob_poroelastic))
+ allocate(velocs_poroelastic_LDDRK(NDIM,nglob_poroelastic))
+ allocate(displw_poroelastic_LDDRK(NDIM,nglob_poroelastic))
+ allocate(velocw_poroelastic_LDDRK(NDIM,nglob_poroelastic))
endif
if(time_stepping_scheme == 3)then
- allocate(accels_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
- allocate(velocs_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
- allocate(accelw_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
- allocate(velocw_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
- allocate(displs_poroelastic_initial_rk(NDIM,nglob_poroelastic))
- allocate(velocs_poroelastic_initial_rk(NDIM,nglob_poroelastic))
- allocate(displw_poroelastic_initial_rk(NDIM,nglob_poroelastic))
- allocate(velocw_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ allocate(accels_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(velocs_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(accelw_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(velocw_poroelastic_rk(NDIM,nglob_poroelastic,stage_time_scheme))
+ allocate(displs_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ allocate(velocs_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ allocate(displw_poroelastic_initial_rk(NDIM,nglob_poroelastic))
+ allocate(velocw_poroelastic_initial_rk(NDIM,nglob_poroelastic))
endif
! extra array if adjoint and kernels calculation
@@ -2870,16 +2901,16 @@
allocate(potential_dot_dot_acoustic(nglob_acoustic))
allocate(rmass_inverse_acoustic(nglob_acoustic))
if(time_stepping_scheme == 2) then
- allocate(potential_acoustic_LDDRK(nglob_acoustic))
- allocate(potential_dot_acoustic_LDDRK(nglob_acoustic))
- allocate(potential_dot_acoustic_temp(nglob_acoustic))
+ allocate(potential_acoustic_LDDRK(nglob_acoustic))
+ allocate(potential_dot_acoustic_LDDRK(nglob_acoustic))
+ allocate(potential_dot_acoustic_temp(nglob_acoustic))
endif
if(time_stepping_scheme == 3) then
- allocate(potential_acoustic_init_rk(nglob_acoustic))
- allocate(potential_dot_acoustic_init_rk(nglob_acoustic))
- allocate(potential_dot_dot_acoustic_rk(nglob_acoustic,stage_time_scheme))
- allocate(potential_dot_acoustic_rk(nglob_acoustic,stage_time_scheme))
+ allocate(potential_acoustic_init_rk(nglob_acoustic))
+ allocate(potential_dot_acoustic_init_rk(nglob_acoustic))
+ allocate(potential_dot_dot_acoustic_rk(nglob_acoustic,stage_time_scheme))
+ allocate(potential_dot_acoustic_rk(nglob_acoustic,stage_time_scheme))
endif
if(SIMULATION_TYPE == 3 .and. any_acoustic) then
@@ -2960,7 +2991,7 @@
allocate(PML_interior_interface(4,1))
endif
-! DK DK add support for using pml in mpi mode with external mesh
+! add support for using PML in MPI mode with external mesh
if(read_external_mesh)then
allocate(mask_ibool_pml(nglob))
else
@@ -3095,7 +3126,7 @@
alpha_z_store(:,:,:) = ZERO
endif
- !elastic PML memory variables
+ ! elastic PML memory variables
if (any_elastic .and. nspec_PML>0) then
allocate(rmemory_displ_elastic(2,3,NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_displ_elastic'
@@ -4627,8 +4658,7 @@ if(coupled_elastic_poro) then
do j = 1, NGLLZ
do i = 1, NGLLX
iglob = ibool(i,j,ispec)
- write(504,'(1pe11.3,1pe11.3,2i3,i7)') &
- coord(1,iglob), coord(2,iglob), i, j, ispec
+ write(504,'(1pe11.3,1pe11.3,2i3,i7)') coord(1,iglob), coord(2,iglob), i, j, ispec
enddo
enddo
enddo
@@ -4636,8 +4666,7 @@ if(coupled_elastic_poro) then
open(unit=504,file='OUTPUT_FILES/mesh_glob',status='unknown',action='write')
do iglob = 1, nglob
- write(504,'(1pe11.3,1pe11.3,i7)') &
- coord(1,iglob), coord(2,iglob), iglob
+ write(504,'(1pe11.3,1pe11.3,i7)') coord(1,iglob), coord(2,iglob), iglob
enddo
close(504)
@@ -4645,8 +4674,7 @@ if(coupled_elastic_poro) then
call create_mask_noise(nglob,coord,mask_noise)
open(unit=504,file='OUTPUT_FILES/mask_noise',status='unknown',action='write')
do iglob = 1, nglob
- write(504,'(1pe11.3,1pe11.3,1pe11.3)') &
- coord(1,iglob), coord(2,iglob), mask_noise(iglob)
+ write(504,'(1pe11.3,1pe11.3,1pe11.3)') coord(1,iglob), coord(2,iglob), mask_noise(iglob)
enddo
close(504)
@@ -4658,8 +4686,7 @@ if(coupled_elastic_poro) then
do i = 1, NGLLX
iglob = ibool(i,j,ispec)
write(504,'(1pe11.3,1pe11.3,1pe11.3,1pe11.3,1pe11.3)') &
- coord(1,iglob), coord(2,iglob), &
- rhoext(i,j,ispec), vpext(i,j,ispec), vsext(i,j,ispec)
+ coord(1,iglob), coord(2,iglob), rhoext(i,j,ispec), vpext(i,j,ispec), vsext(i,j,ispec)
enddo
enddo
enddo
More information about the CIG-COMMITS
mailing list