[cig-commits] r8582 - in seismo/2D/SPECFEM2D/trunk: . DATA
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:58:22 PST 2007
Author: walter
Date: 2007-12-07 15:58:21 -0800 (Fri, 07 Dec 2007)
New Revision: 8582
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct
seismo/2D/SPECFEM2D/trunk/attenuation_model.f90
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/constants.h
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
added N_SLS (number of standard linear solids for attenuation) in Par_file, to allow the user to change its value without having to recompile specfem2D; it is no longer in constants.h. Par_file has now changed accordingly.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-07 23:58:21 UTC (rev 8582)
@@ -50,7 +50,8 @@
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
-# constants for attenuation
+# 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
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_unstruct 2007-12-07 23:58:21 UTC (rev 8582)
@@ -51,6 +51,7 @@
factor = 1.d10 # amplification factor
# 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
f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
Modified: seismo/2D/SPECFEM2D/trunk/attenuation_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/attenuation_model.f90 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/attenuation_model.f90 2007-12-07 23:58:21 UTC (rev 8582)
@@ -11,7 +11,7 @@
!
!========================================================================
- subroutine attenuation_model(Qp_attenuation,Qs_attenuation,f0_attenuation, &
+ subroutine 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 constants
@@ -20,6 +20,7 @@
include "constants.h"
+ integer :: N_SLS
double precision :: Qp_attenuation,Qs_attenuation,f0_attenuation
double precision, dimension(N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
double precision :: Mu_nu1,Mu_nu2
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2007-12-07 23:58:21 UTC (rev 8582)
@@ -16,7 +16,7 @@
nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
vpext,vsext,rhoext,wxgll,wzgll,numat, &
pressure_element,vector_field_element,e1,e11, &
- potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute kinetic and potential energy in the solid (acoustic elements are excluded)
@@ -28,11 +28,12 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLX) :: vector_field_element
! pressure in an element
+ integer :: N_SLS
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
-
+
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_acoustic,potential_dot_dot_acoustic
logical :: TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
@@ -172,7 +173,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, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ 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, &
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-12-07 23:58:21 UTC (rev 8582)
@@ -19,7 +19,7 @@
jacobian,vpext,vsext,rhoext,source_time_function,sourcearray,e1,e11, &
e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2, &
+ hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
nspec_inner_outer, ispec_inner_outer_to_glob, num_phase_inner_outer &
)
@@ -50,6 +50,7 @@
real(kind=CUSTOM_REAL), dimension(NSTEP) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLZ) :: sourcearray
+ 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
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2007-12-07 23:58:21 UTC (rev 8582)
@@ -14,7 +14,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, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute pressure in acoustic elements and in elastic elements
@@ -44,6 +44,7 @@
logical :: assign_external_model,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
+ integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
double precision :: Mu_nu1,Mu_nu2
@@ -60,7 +61,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, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ 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
do j = 1,NGLLZ
@@ -81,7 +82,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, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute pressure in acoustic elements and in elastic elements
@@ -113,6 +114,7 @@
logical :: assign_external_model,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
+ 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
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-07 23:58:21 UTC (rev 8582)
@@ -14,9 +14,6 @@
! the code does NOT work if NGLLZ /= NGLLX because it then cannot handle a non-structured mesh
integer, parameter :: NGLLZ = NGLLX
-! number of standard linear solids for attenuation
- integer, parameter :: N_SLS = 2
-
! compute and output acoustic and elastic energy (slows down the code significantly)
logical, parameter :: OUTPUT_ENERGY = .false.
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-07 23:58:21 UTC (rev 8582)
@@ -166,6 +166,7 @@
character(len=256) :: prname
! variables used for attenuation
+ integer :: N_SLS
double precision :: Qp_attenuation
double precision :: Qs_attenuation
double precision :: f0_attenuation
@@ -413,6 +414,7 @@
print *,'Multiplying factor = ',factor
! 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)
@@ -1097,7 +1099,7 @@
write(15,*) source_type,time_function_type,xs,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
write(15,*) 'attenuation'
- write(15,*) Qp_attenuation, Qs_attenuation, f0_attenuation
+ write(15,*) N_SLS, Qp_attenuation, Qs_attenuation, f0_attenuation
write(15,*) 'Coordinates of macrobloc mesh (coorg):'
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-09-12 14:49:41 UTC (rev 8581)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-07 23:58:21 UTC (rev 8582)
@@ -188,6 +188,7 @@
logical, dimension(:,:), allocatable :: codeabs
! for attenuation
+ integer :: N_SLS
double precision :: Qp_attenuation
double precision :: Qs_attenuation
double precision :: f0_attenuation
@@ -195,7 +196,7 @@
double precision :: deltatsquare,deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
- double precision, dimension(N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
+ double precision, dimension(:), allocatable :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
double precision :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
@@ -417,7 +418,7 @@
!---- read attenuation information
!
read(IIN,"(a80)") datlin
- read(IIN,*) Qp_attenuation, Qs_attenuation, f0_attenuation
+ read(IIN,*) N_SLS, Qp_attenuation, Qs_attenuation, f0_attenuation
!
!----- check the input
@@ -484,6 +485,10 @@
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 arrays for absorbing boundary conditions
if(nelemabs <= 0) then
@@ -568,7 +573,7 @@
allocate(dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
! define the attenuation constants
- call attenuation_model(Qp_attenuation,Qs_attenuation,f0_attenuation, &
+ 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)
!
@@ -1913,7 +1918,7 @@
jacobian,vpext,vsext,rhoext,source_time_function,sourcearray, &
e1,e11,e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2, &
+ hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
nspec_outer, ispec_outer_to_glob, .true. &
)
@@ -2012,7 +2017,7 @@
jacobian,vpext,vsext,rhoext,source_time_function,sourcearray, &
e1,e11,e13,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
- hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2, &
+ hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
nspec_inner, ispec_inner_to_glob, .false. &
)
@@ -2049,7 +2054,7 @@
nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
vpext,vsext,rhoext,wxgll,wzgll,numat, &
pressure_element,vector_field_element,e1,e11, &
- potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5) then
@@ -2091,7 +2096,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, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
else if(.not. elastic(ispec)) then
@@ -2281,7 +2286,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, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2)
+ TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
else
call exit_MPI('wrong type for snapshots')
More information about the cig-commits
mailing list