[cig-commits] r14814 - in seismo/2D/SPECFEM2D/trunk: . DATA
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Apr 29 11:41:10 PDT 2009
Author: dkomati1
Date: 2009-04-29 11:41:10 -0700 (Wed, 29 Apr 2009)
New Revision: 14814
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/compute_energy.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
seismo/2D/SPECFEM2D/trunk/gmat01.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt
Log:
merged more general attenuation (quality factors can now be given for each spectral element rather
than for the whole mesh only) developed by Kasper van Wijk and Dylan Mikesell.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2009-04-29 18:41:10 UTC (rev 14814)
@@ -16,7 +16,7 @@
nproc = 1 # number of processes
partitioning_method = 1 # ascending order = 1, Metis = 2, Scotch = 3
partitioning_strategy = 01110 #b{strat=m{asc=g,low=g,rat=1.0,type=h,vert=4}} #01110 # options concerning partitioning strategy.
-
+
# geometry of the model (origin lower-left corner = 0,0) and mesh description
xmin = 0.d0 # abscissa of left side of the model
xmax = 4000.d0 # abscissa of right side of the model
@@ -30,7 +30,7 @@
# absorbing boundaries parameters
-absorbing_conditions = .true. # absorbing boundary active or not
+absorbing_conditions = .true. # absorbing boundary active or not
absorbbottom = .true.
absorbright = .true.
absorbtop = .false.
@@ -55,9 +55,11 @@
force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
# constants for attenuation
-N_SLS = 2 # number of standard linear solids for attenuation
-Qp_attenuation = 136.4376068115 # quality factor P for attenuation
-Qs_attenuation = 136.4376068115 # quality factor S for attenuation
+N_SLS = 2 # number of standard linear solids for attenuation
+## DK DK Qp and Qs can now vary in each spectral element and are therefore given in the same list
+## DK DK list as rho, Vp and Vs at the end of this file
+#Qp_attenuation = 136.4376068115 # quality factor P for attenuation
+#Qs_attenuation = 136.4376068115 # quality factor S for attenuation
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
# receiver line parameters for seismograms
@@ -93,13 +95,13 @@
# velocity and density models
nbmodels = 4 # nb of different models
-# define models as (model_number,1,rho,vp,vs,0,0) or (model_number,2,rho,c11,c13,c33,c44)
-# set vs to zero to make a given model acoustic
+# define models as (model_number,1,rho,Vp,Vs,0,0,Qp,Qs) or (model_number,2,rho,c11,c13,c33,c44,Qp,Qs)
+# set Vs to zero to make a given model acoustic
# the mesh can contain both acoustic and elastic models simultaneously
-1 1 2700.d0 3000.d0 1732.051d0 0 0
-2 1 2500.d0 2700.d0 0 0 0 #1558.89d0 0 0
-3 1 2200.d0 2500.d0 1443.375d0 0 0
-4 1 2200.d0 2200.d0 1343.375d0 0 0
+1 1 2700.d0 3000.d0 1732.051d0 0 0 10.d0 10.d0
+2 1 2500.d0 2700.d0 0 0 0 136.d0 136.d0 #1558.89d0 0 0 136.d0 136.d0
+3 1 2200.d0 2500.d0 1443.375d0 0 0 136.d0 136.d0
+4 1 2200.d0 2200.d0 1343.375d0 0 0 136.d0 136.d0
# define the different regions of the model in the (nx,nz) spectral element mesh
nbregions = 4 # nb of regions and model number for each
1 80 1 20 1
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2009-04-29 18:41:10 UTC (rev 14814)
@@ -63,7 +63,7 @@
#F90 = ifort
#F90 = mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
#CC = gcc
-#FLAGS_NOCHECK=-O3 -xP -vec-report0 -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz
+#FLAGS_NOCHECK=-O3 -xP -vec-report0 -e95 -std95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz
#FLAGS_CHECK = $(FLAGS_NOCHECK)
# GNU gfortran
@@ -72,7 +72,7 @@
#F90 = /opt/openmpi-1.2.1/gfortran64/bin/mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
CC = gcc
#FLAGS_NOCHECK = -O3 -march=opteron -m64 -mfpmath=sse,387
-FLAGS_NOCHECK = -std=f2003 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -mcmodel=medium
+FLAGS_NOCHECK = -std=f95 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -mcmodel=medium
FLAGS_CHECK = $(FLAGS_NOCHECK) -fbounds-check
# IBM
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2009-04-29 18:41:10 UTC (rev 14814)
@@ -61,7 +61,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: pressure_element
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
- double precision :: Mu_nu1,Mu_nu2
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_acoustic,potential_dot_dot_acoustic
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2009-04-29 18:41:10 UTC (rev 14814)
@@ -83,8 +83,8 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11,e13
- double precision, dimension(N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
- double precision :: Mu_nu1,Mu_nu2
+ double precision, dimension(NGLLX,NGLLZ,nspec,N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL) :: e1_sum,e11_sum,e13_sum
integer :: i_sls
@@ -212,8 +212,8 @@
! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611 (1988).
! compute unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
- lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1 - mul_relaxed * Mu_nu2
- mul_unrelaxed = mul_relaxed * Mu_nu2
+ lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1(i,j,ispec) - mul_relaxed * Mu_nu2(i,j,ispec)
+ mul_unrelaxed = mul_relaxed * Mu_nu2(i,j,ispec)
lambdalplus2mul_unrelaxed = lambdal_unrelaxed + TWO*mul_unrelaxed
! compute the stress using the unrelaxed Lame parameters (Carcione 1993, page 111)
@@ -636,12 +636,12 @@
! evolution e1
Un = e1(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu1(i_sls)
+ tauinv = - inv_tau_sigma_nu1(i,j,ispec,i_sls)
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = theta_n * phi_nu1(i_sls)
- Snp1 = theta_np1 * phi_nu1(i_sls)
+ Sn = theta_n * phi_nu1(i,j,ispec,i_sls)
+ Snp1 = theta_np1 * phi_nu1(i,j,ispec,i_sls)
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
@@ -650,12 +650,12 @@
! evolution e11
Un = e11(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu2(i_sls)
+ tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i_sls)
- Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i_sls)
+ Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
@@ -664,12 +664,12 @@
! evolution e13
Un = e13(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu2(i_sls)
+ tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
tauinvsquare = tauinv * tauinv
tauinvcube = tauinvsquare * tauinv
tauinvUn = tauinv * Un
- Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2(i_sls)
- Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2(i_sls)
+ Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2009-04-29 18:41:10 UTC (rev 14814)
@@ -74,7 +74,7 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
- double precision :: Mu_nu1,Mu_nu2
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
! local variables
integer :: i,j,ispec,iglob
@@ -144,7 +144,7 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
real(kind=CUSTOM_REAL) :: e1_sum,e11_sum
- double precision :: Mu_nu1,Mu_nu2
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
integer :: i_sls
! local variables
@@ -239,8 +239,8 @@
! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611 (1988).
! compute unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
- lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1 - mul_relaxed * Mu_nu2
- mul_unrelaxed = mul_relaxed * Mu_nu2
+ lambdal_unrelaxed = (lambdal_relaxed + mul_relaxed) * Mu_nu1(i,j,ispec) - mul_relaxed * Mu_nu2(i,j,ispec)
+ mul_unrelaxed = mul_relaxed * Mu_nu2(i,j,ispec)
lambdalplus2mul_unrelaxed = lambdal_unrelaxed + TWO*mul_unrelaxed
! compute the stress using the unrelaxed Lame parameters (Carcione 1993, page 111)
Modified: seismo/2D/SPECFEM2D/trunk/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/gmat01.f90 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/gmat01.f90 2009-04-29 18:41:10 UTC (rev 14814)
@@ -40,7 +40,7 @@
!
!========================================================================
- subroutine gmat01(density_array,elastcoef,numat,myrank,ipass)
+ subroutine gmat01(density_array,elastcoef,numat,myrank,ipass,Qp_array,Qs_array)
! read properties of a 2D isotropic or anisotropic linear elastic element
@@ -52,11 +52,11 @@
double precision lambdaplus2mu,kappa
integer numat,myrank,ipass
- double precision density_array(numat),elastcoef(4,numat)
+ double precision density_array(numat),elastcoef(4,numat),Qp_array(numat),Qs_array(numat)
integer in,n,indic
- double precision young,poisson,density,cp,cs,mu,two_mu,lambda
- double precision val1,val2,val3,val4
+ double precision young,poisson,density,cp,cs,mu,two_mu,lambda,Qp,Qs
+ double precision val1,val2,val3,val4,val5,val6
double precision c11,c13,c33,c44
!
@@ -64,13 +64,15 @@
!
density_array(:) = zero
elastcoef(:,:) = zero
+ Qp_array(:) = zero
+ Qs_array(:) = zero
if(myrank == 0 .and. ipass == 1) write(IOUT,100) numat
read(IIN,"(a80)") datlin
do in = 1,numat
- read(IIN,*) n,indic,density,val1,val2,val3,val4
+ read(IIN,*) n,indic,density,val1,val2,val3,val4,val5,val6
if(n<1 .or. n>numat) call exit_MPI('Wrong material set number')
@@ -81,6 +83,10 @@
cp = val1
cs = val2
+! Qp and Qs values
+ Qp = val5
+ Qs = val6
+
! Lam'e parameters
lambdaplus2mu = density*cp*cp
mu = density*cs*cs
@@ -130,6 +136,8 @@
endif
density_array(n) = density
+ Qp_array(n) = Qp
+ Qs_array(n) = Qs
!
!---- check what has been read
@@ -138,12 +146,12 @@
if(indic == 1) then
! material can be acoustic (fluid) or elastic (solid)
if(elastcoef(2,n) > TINYVAL) then
- write(IOUT,200) n,cp,cs,density,poisson,lambda,mu,kappa,young
+ write(IOUT,200) n,cp,cs,density,poisson,lambda,mu,kappa,young,Qp,Qs
else
- write(IOUT,300) n,cp,density,kappa
+ write(IOUT,300) n,cp,density,kappa,Qp,Qs
endif
else
- write(IOUT,400) n,c11,c13,c33,c44,density,sqrt(c33/density),sqrt(c11/density),sqrt(c44/density),sqrt(c44/density)
+ write(IOUT,400) n,c11,c13,c33,c44,density,sqrt(c33/density),sqrt(c11/density),sqrt(c44/density),sqrt(c44/density),Qp,Qs
endif
endif
@@ -167,7 +175,9 @@
'First Lame parameter Lambda. . . (lambda) =',1pe15.8,/5x, &
'Second Lame parameter Mu. . . . . . .(mu) =',1pe15.8,/5x, &
'Bulk modulus Kappa . . . . . . . .(kappa) =',1pe15.8,/5x, &
- 'Young''s modulus E. . . . . . . . .(young) =',1pe15.8)
+ 'Young''s modulus E. . . . . . . . .(young) =',1pe15.8,/5x, &
+ 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
+ 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
300 format(//5x,'-------------------------------',/5x, &
'-- Acoustic (fluid) material --',/5x, &
@@ -175,7 +185,9 @@
'Material set number. . . . . . . . (jmat) =',i6,/5x, &
'P-wave velocity. . . . . . . . . . . (cp) =',1pe15.8,/5x, &
'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
- 'Bulk modulus Kappa . . . . . . . .(kappa) =',1pe15.8)
+ 'Bulk modulus Kappa . . . . . . . .(kappa) =',1pe15.8,/5x, &
+ 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
+ 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
400 format(//5x,'-------------------------------------',/5x, &
'-- Transverse anisotropic material --',/5x, &
@@ -189,7 +201,9 @@
'Velocity of qP along vertical axis. . . . =',1pe15.8,/5x, &
'Velocity of qP along horizontal axis. . . =',1pe15.8,/5x, &
'Velocity of qSV along vertical axis . . . =',1pe15.8,/5x, &
- 'Velocity of qSV along horizontal axis . . =',1pe15.8)
+ 'Velocity of qSV along horizontal axis . . =',1pe15.8,/5x, &
+ 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
+ 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
end subroutine gmat01
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2009-04-29 18:41:10 UTC (rev 14814)
@@ -96,7 +96,7 @@
double precision :: gamma,absx,a00,a01,bot0,top0
! to store density and velocity model
- double precision, dimension(:), allocatable :: rho,cp,cs,aniso3,aniso4
+ double precision, dimension(:), allocatable :: rho,cp,cs,aniso3,aniso4,Qp,Qs
integer, dimension(:), allocatable :: icodemat
integer, dimension(:), allocatable :: num_material
@@ -131,7 +131,7 @@
double precision tang1,tangN,vpregion,vsregion,poisson_ratio
double precision cutsnaps,sizemax_arrows,anglerec,xmin,xmax,deltat
- double precision rhoread,cpread,csread,aniso3read,aniso4read
+ double precision rhoread,cpread,csread,aniso3read,aniso4read,Qpread,Qsread
double precision, dimension(:), allocatable :: xdeb,zdeb,xfin,zfin
@@ -210,8 +210,6 @@
! variables used for attenuation
integer :: N_SLS
- double precision :: Qp_attenuation
- double precision :: Qs_attenuation
double precision :: f0_attenuation
! variables used for tangential detection
@@ -467,8 +465,6 @@
! read constants for attenuation
call read_value_integer(IIN,IGNORE_JUNK,N_SLS)
- call read_value_double_precision(IIN,IGNORE_JUNK,Qp_attenuation)
- call read_value_double_precision(IIN,IGNORE_JUNK,Qs_attenuation)
call read_value_double_precision(IIN,IGNORE_JUNK,f0_attenuation)
! if source is not a Dirac or Heavyside then f0_attenuation is f0
@@ -502,7 +498,7 @@
call read_value_double_precision(IIN,IGNORE_JUNK,zfin(ireceiverlines))
call read_value_logical(IIN,IGNORE_JUNK,enreg_surf(ireceiverlines))
if (read_external_mesh .and. enreg_surf(ireceiverlines)) then
- stop 'Cannot use enreg_surf with external meshes!'
+ stop 'Cannot use enreg_surf with external meshes!'
endif
enddo
@@ -540,6 +536,8 @@
allocate(cs(nb_materials))
allocate(aniso3(nb_materials))
allocate(aniso4(nb_materials))
+ allocate(Qp(nb_materials))
+ allocate(Qs(nb_materials))
allocate(num_material(nelmnts))
icodemat(:) = 0
@@ -548,17 +546,22 @@
cs(:) = 0.d0
aniso3(:) = 0.d0
aniso4(:) = 0.d0
+ Qp(:) = 0.d0
+ Qs(:) = 0.d0
num_material(:) = 0
do imaterial=1,nb_materials
- call read_material_parameters(IIN,DONT_IGNORE_JUNK,i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read)
+ call read_material_parameters(IIN,DONT_IGNORE_JUNK,i,icodematread,rhoread,Qpread,Qsread,cpread,csread,aniso3read,aniso4read)
if(i < 1 .or. i > nb_materials) stop 'Wrong material number!'
icodemat(i) = icodematread
rho(i) = rhoread
cp(i) = cpread
cs(i) = csread
+ Qp(i) = Qpread
+ Qs(i) = Qsread
if(rho(i) <= 0.d0 .or. cp(i) <= 0.d0 .or. cs(i) < 0.d0) stop 'negative value of velocity or density'
+ if(Qp(i) <= 0.d0 .or. Qs(i) <= 0.d0) stop 'negative value of Qp or Qs'
aniso3(i) = aniso3read
aniso4(i) = aniso4read
@@ -570,7 +573,7 @@
do i=1,nb_materials
if(icodemat(i) /= ANISOTROPIC_MATERIAL) then
print *,'Material #',i,' isotropic'
- print *,'rho,cp,cs = ',rho(i),cp(i),cs(i)
+ print *,'rho,cp,cs = ',rho(i),cp(i),cs(i),Qp(i),Qs(i)
if(cs(i) < TINYVAL) then
print *,'Material is fluid'
else
@@ -578,7 +581,7 @@
endif
else
print *,'Material #',i,' anisotropic'
- print *,'rho,c11,c13,c33,c44 = ',rho(i),cp(i),cs(i),aniso3(i),aniso4(i)
+ print *,'rho,c11,c13,c33,c44 = ',rho(i),cp(i),cs(i),aniso3(i),aniso4(i),Qp(i),Qs(i)
endif
print *
enddo
@@ -638,6 +641,8 @@
poisson_ratio = 0.5d0*(vpregion*vpregion-2.d0*vsregion*vsregion) / (vpregion*vpregion-vsregion*vsregion)
print *,'Poisson''s ratio = ',poisson_ratio
if(poisson_ratio <= -1.00001d0 .or. poisson_ratio >= 0.50001d0) stop 'incorrect value of Poisson''s ratio'
+ print *,'Qp = ',Qp(imaterial_number)
+ print *,'Qs = ',Qs(imaterial_number)
else
print *,'Material # ',imaterial_number,' anisotropic'
print *,'c11 = ',cp(imaterial_number)
@@ -645,6 +650,8 @@
print *,'c33 = ',aniso3(imaterial_number)
print *,'c44 = ',aniso4(imaterial_number)
print *,'rho = ',rho(imaterial_number)
+ print *,'Qp = ',Qp(imaterial_number)
+ print *,'Qs = ',Qs(imaterial_number)
endif
print *,' -----'
@@ -1168,20 +1175,17 @@
write(15,*) source_type,time_function_type,xs,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
write(15,*) 'attenuation'
- write(15,*) N_SLS, Qp_attenuation, Qs_attenuation, f0_attenuation
+ write(15,*) N_SLS, f0_attenuation
write(15,*) 'Coordinates of macrobloc mesh (coorg):'
-
call write_glob2loc_nodes_database(15, iproc, npgeo, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, nnodes, 2)
-
write(15,*) 'numat ngnod nspec pointsdisp plot_lowerleft_corner_only'
write(15,*) nb_materials,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only
-
- if ( any_abs ) then
+ if (any_abs) then
call write_abs_merge_database(15, nelemabs_merge, nelemabs_loc, &
abs_surface_char, abs_surface_merge, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
@@ -1191,38 +1195,35 @@
nelemabs_loc = 0
endif
- call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
+ call write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
iproc, glob2loc_elmnts, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 1)
-
call write_fluidsolid_edges_database(15, nedges_coupled, nedges_coupled_loc, &
edges_coupled, glob2loc_elmnts, part, iproc, 1)
write(15,*) 'nelemabs nelem_acoustic_surface num_fluid_solid_edges nnodes_tangential_curve'
write(15,*) nelemabs_loc,nelem_acoustic_surface_loc,nedges_coupled_loc,nnodes_tangential_curve
-
- write(15,*) 'Material sets (num 1 rho vp vs 0 0) or (num 2 rho c11 c13 c33 c44)'
+ write(15,*) 'Material sets (num 1 rho vp vs 0 0 Qp Qs) or (num 2 rho c11 c13 c33 c44 Qp Qs)'
do i=1,nb_materials
- write(15,*) i,icodemat(i),rho(i),cp(i),cs(i),aniso3(i),aniso4(i)
+ write(15,*) i,icodemat(i),rho(i),cp(i),cs(i),aniso3(i),aniso4(i),Qp(i),Qs(i)
enddo
write(15,*) 'Arrays kmato and knods for each bloc:'
-
call write_partition_database(15, iproc, nspec, nelmnts, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 2)
if ( nproc /= 1 ) then
- call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
+ call write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, 1)
write(15,*) 'Interfaces:'
write(15,*) my_ninterface, maxval(my_nb_interfaces)
- call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
+ call write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nproc, iproc, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, 2)
@@ -1242,7 +1243,7 @@
endif
write(15,*) 'List of acoustic free-surface elements:'
- call Write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
+ call write_surface_database(15, nelem_acoustic_surface, acoustic_surface, nelem_acoustic_surface_loc, &
iproc, glob2loc_elmnts, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part, 2)
Modified: seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2009-04-29 18:41:10 UTC (rev 14814)
@@ -139,18 +139,18 @@
!--------------------
- subroutine read_material_parameters(iin,ignore_junk,i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read)
+ subroutine read_material_parameters(iin,ignore_junk,i,icodematread,rhoread,Qpread,Qsread,cpread,csread,aniso3read,aniso4read)
implicit none
integer iin
logical ignore_junk
integer i,icodematread
- double precision rhoread,cpread,csread,aniso3read,aniso4read
+ double precision rhoread,Qpread,Qsread,cpread,csread,aniso3read,aniso4read
character(len=100) string_read
call read_next_line(iin,ignore_junk,string_read)
- read(string_read,*) i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read
+ read(string_read,*) i,icodematread,rhoread,Qpread,Qsread,cpread,csread,aniso3read,aniso4read
end subroutine read_material_parameters
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2009-04-29 18:41:10 UTC (rev 14814)
@@ -249,15 +249,15 @@
! for attenuation
integer :: N_SLS
- double precision :: Qp_attenuation
- double precision :: Qs_attenuation
+ double precision, dimension(:), allocatable :: Qp_attenuation
+ double precision, dimension(:), allocatable :: Qs_attenuation
double precision :: f0_attenuation
integer nspec_allocate
double precision :: deltatsquare,deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
- double precision, dimension(:), allocatable :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
- double precision :: Mu_nu1,Mu_nu2
+ double precision, dimension(:,:,:,:), allocatable :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
+ double precision, dimension(:,:,:) , allocatable :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
@@ -575,7 +575,7 @@
!---- read attenuation information
!
read(IIN,"(a80)") datlin
- read(IIN,*) N_SLS, Qp_attenuation, Qs_attenuation, f0_attenuation
+ read(IIN,*) N_SLS, f0_attenuation
!
!----- check the input
@@ -651,14 +651,16 @@
allocate(Uzinterp(pointsdisp,pointsdisp))
allocate(density(numat))
allocate(elastcoef(4,numat))
+ allocate(Qp_attenuation(numat))
+ allocate(Qs_attenuation(numat))
allocate(kmato(nspec))
allocate(knods(ngnod,nspec))
allocate(ibool(NGLLX,NGLLZ,nspec))
allocate(elastic(nspec))
- allocate(inv_tau_sigma_nu1(N_SLS))
- allocate(inv_tau_sigma_nu2(N_SLS))
- allocate(phi_nu1(N_SLS))
- allocate(phi_nu2(N_SLS))
+ allocate(inv_tau_sigma_nu1(NGLLX,NGLLZ,nspec,N_SLS))
+ 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))
endif
! --- allocate arrays for absorbing boundary conditions
@@ -697,7 +699,7 @@
!
!---- read the material properties
!
- call gmat01(density,elastcoef,numat,myrank,ipass)
+ call gmat01(density,elastcoef,numat,myrank,ipass,Qp_attenuation,Qs_attenuation)
!
!---- read spectral macrobloc data
@@ -758,11 +760,22 @@
allocate(duz_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
+ allocate(Mu_nu1(NGLLX,NGLLZ,nspec))
+ allocate(Mu_nu2(NGLLX,NGLLZ,nspec))
endif
-! define the attenuation constants
- call attenuation_model(N_SLS,Qp_attenuation,Qs_attenuation,f0_attenuation, &
- inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2)
+! define the attenuation quality factors.
+! they can be different for each element.
+!! DK DK if needed in the future, here the quality factor could be different for each point
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ call attenuation_model(N_SLS,Qp_attenuation(kmato(ispec)),Qs_attenuation(kmato(ispec)), &
+ f0_attenuation,inv_tau_sigma_nu1(i,j,ispec,:),phi_nu1(i,j,ispec,:), &
+ inv_tau_sigma_nu2(i,j,ispec,:),phi_nu2(i,j,ispec,:),Mu_nu1(i,j,ispec),Mu_nu2(i,j,ispec))
+ enddo
+ enddo
+ enddo
!
!---- read interfaces data
@@ -1326,7 +1339,7 @@
if ( dist_current < distmin ) then
n1_tangential_detection_curve = i
distmin = dist_current
-
+
endif
enddo
@@ -1356,11 +1369,11 @@
call MPI_send(angleforce,1,MPI_DOUBLE_PRECISION,0,43,MPI_COMM_WORLD,ier)
#endif
end if
-
+
#ifdef USE_MPI
call MPI_bcast(angleforce_recv,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
angleforce = angleforce_recv
-#endif
+#endif
endif
! compute distance from source to receivers following the curve
@@ -1404,19 +1417,19 @@
if ( myrank == 0 ) then
open(unit=11,file='OUTPUT_FILES/dist_rec_tangential_detection_curve', form='formatted', status='unknown')
- end if
+ end if
irecloc = 0
do irec = 1,nrec
-
+
if ( myrank == 0 ) then
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
n1_tangential_detection_curve = rec_tangential_detection_curve(irecloc)
x_final_receiver_dummy = x_final_receiver(irec)
z_final_receiver_dummy = z_final_receiver(irec)
-#ifdef USE_MPI
+#ifdef USE_MPI
else
-
+
call MPI_RECV(n1_tangential_detection_curve,1,MPI_INTEGER,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,request_mpi_status,ier)
call MPI_RECV(x_final_receiver_dummy,1,MPI_DOUBLE_PRECISION,&
@@ -1436,7 +1449,7 @@
call MPI_SEND(z_final_receiver(irec),1,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ier)
end if
#endif
-
+
end if
if ( myrank == 0 ) then
write(11,*) dist_tangential_detection_curve(n1_tangential_detection_curve)
@@ -3600,27 +3613,27 @@
subroutine tri_quad(n, n1, nnodes)
-
+
implicit none
integer :: n1, n2, n3, n4, nnodes
integer, dimension(4) :: n
-
+
n(2) = n1
-
+
if ( n1 == 1 ) then
n(1) = nnodes
else
n(1) = n1-1
- endif
-
+ endif
+
if ( n1 == nnodes ) then
n(3) = 1
else
n(3) = n1+1
endif
-
+
if ( n(3) == nnodes ) then
n(4) = 1
else
@@ -3632,15 +3645,15 @@
subroutine calcul_normale( angle, n1_x, n2_x, n3_x, n4_x, n1_z, n2_z, n3_z, n4_z )
-
+
implicit none
-
+
include 'constants.h'
double precision :: angle
double precision :: n1_x, n2_x, n3_x, n4_x, n1_z, n2_z, n3_z, n4_z
-
- double precision :: theta1, theta2, theta3
+
+ double precision :: theta1, theta2, theta3
double precision :: costheta1, costheta2, costheta3
if ( abs(n2_z - n1_z) < TINYVAL ) then
@@ -3658,12 +3671,12 @@
else
costheta3 = (n4_z - n3_z) / sqrt((n4_x - n3_x)**2 + (n4_z - n3_z)**2)
endif
-
+
theta1 = - sign(1.d0,n2_x - n1_x) * acos(costheta1)
theta2 = - sign(1.d0,n3_x - n2_x) * acos(costheta2)
theta3 = - sign(1.d0,n4_x - n3_x) * acos(costheta3)
-
-
+
+
angle = angle + ( theta1 + theta2 + theta3 ) / 3.d0 + PI/2.d0
!angle = theta2
Modified: seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt 2009-04-29 04:11:16 UTC (rev 14813)
+++ seismo/2D/SPECFEM2D/trunk/todo_list_please_dont_remove.txt 2009-04-29 18:41:10 UTC (rev 14814)
@@ -1,7 +1,14 @@
-April 2009: see list of problems (in particular in the DATA/Par_file format)
-found by Steve Smith below. Pieyre Le Loher will fix them at some point.
+IMPORTANT KNOWN BUG: in compute_forces_elastic.f90 the calculation of
+attenuation (viscoelasticity) is slightly incorrect because the gradients
+are computed twice but *at the same time step* instead of at different
+(staggered) time steps (t_{n-1} and t_n or something like that).
+That's easy to fix but I have no time for now. Let us fix that in the future.
+It will only make a very small difference in the final seismograms therefore
+the current code with the bug can be used without any major problem.
+Dimitri Komatitsch, April 28, 2009.
+
---------------------------
- improving compiling with SCOTCH (issue with header file scotchf.h which is Fortran77 legal). Having our own scotchf.h file (without the comments) is not wise.
@@ -39,6 +46,18 @@
------------------------------
+SOMETHING THAT COULD BE MADE MORE GENERAL:
+
+at line 769 of specfem2D.F90:
+!! DK DK if needed in the future, here the quality factor could be different for each point
+
+i.e. they could be given at each (i,j,ispec) instead of at each (ispec) only
+in the current version. Very easy to do if needed, just that line to change.
+
+Dimitri Komatitsch, April 28, 2009.
+
+------------------------------
+
April 2009: here is the list of problems (in particular in the DATA/Par_file format)
found by Steve Smith below. Pieyre Le Loher will fix them at some point:
@@ -94,7 +113,7 @@
>> From what you tell me in your previous email (below), I understand
>> that:
>>
->> 1) Of all the Par_files included in the current release, *ONLY*
+>> 1) Of all the Par_files included in the current release, *ONLY*
>> DATA/Par_file works. None of the other Par_files/examples work
>> with the code.
>>
More information about the CIG-COMMITS
mailing list