[cig-commits] r16405 - in seismo/2D/SPECFEM2D/trunk: . DATA
pieyre at geodynamics.org
pieyre at geodynamics.org
Fri Mar 12 12:23:27 PST 2010
Author: pieyre
Date: 2010-03-12 12:23:22 -0800 (Fri, 12 Mar 2010)
New Revision: 16405
Added:
seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90
seismo/2D/SPECFEM2D/trunk/read_external_model.f90
seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/SOURCE
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
seismo/2D/SPECFEM2D/trunk/checkgrid.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/define_external_model.f90
seismo/2D/SPECFEM2D/trunk/gmat01.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
seismo/2D/SPECFEM2D/trunk/plotpost.F90
seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
modified anisotropy, added attenuation and anisotropy for external models, fixed some bugs, put some parts of specfem2D.F90 in 3 new routines and added an ifort flag to avoid optimisation loss.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2010-03-12 20:23:22 UTC (rev 16405)
@@ -25,7 +25,6 @@
initialfield = .false. # use a plane wave as source or not
add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
assign_external_model = .false. # define external earth model or not
-TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off for solid medium
TURN_ATTENUATION_ON = .false. # turn attenuation on or off for solid medium
TURN_VISCATTENUATION_ON = .false. # turn viscous attenuation on or off
Q0 = 1 # quality factor for viscous attenuation
@@ -89,19 +88,21 @@
outputgrid = .false. # save the grid in a text file or not
# velocity and density models
-nbmodels = 4 # nb of different models
-# define models as I: (model_number,1,rho,Vp,Vs,0,0,Qp,Qs) or II: (model_number,2,rho,c11,c13,c33,c44,Qp,Qs) or III: (model_number,3,rhos,rhof,phi,c,kxx,kxz,kzz,Ks,Kf,Kfr,etaf,mufr,Qs).
+nbmodels = 5 # nb of different models
+# define models as I: (model_number,1,rho,Vp,Vs,0,0,Qp,Qs) or II: (model_number,2,rho,c11,c13,c15,c33,c35,c55) or III: (model_number,3,rhos,rhof,phi,c,kxx,kxz,kzz,Ks,Kf,Kfr,etaf,mufr,Qs).
# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, & poroelastic models simultaneously
1 1 1000.d0 1500.d0 0.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
2 1 2500.d0 5000.d0 2500.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
3 1 2200.d0 2500.d0 1443.375d0 0 0 10.d0 10.d0 0 0 0 0 0 0
4 3 2200.d0 786.3d0 0.4 2.0 1d-11 0.0 1d-11 5.341d9 2d9 3d9 0.0d-4 3.204d9 10.d0
+5 2 2500.d0 169.d9 122.d9 0.d0 169.d9 0.d0 75.3d9 0 0 0 0 0 0
#4 1 2200.d0 2200.d0 1343.375d0 0 0 10.d0 10.d0 0 0 0 0 0 0
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 5 # nb of regions and model number for each
+nbregions = 6 # nb of regions and model number for each
1 80 1 20 1
1 80 21 40 4
1 80 41 60 3
60 80 21 40 4
30 40 50 60 2
+35 40 50 60 5
Modified: seismo/2D/SPECFEM2D/trunk/DATA/SOURCE
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/SOURCE 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/DATA/SOURCE 2010-03-12 20:23:22 UTC (rev 16405)
@@ -1,4 +1,4 @@
-#source 1
+#source 1. The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
source_surf = .false. # source inside the medium or at the surface
xs = 0. # source location x in meters
zs = -100. # source location z in meters
@@ -11,7 +11,3 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
-
-# The components of a moment tensor source must be given in N.m,
-# not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
-
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2010-03-12 20:23:22 UTC (rev 16405)
@@ -65,7 +65,7 @@
F90 = ifort
#F90 = mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
CC = gcc
-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_NOCHECK=-O3 -xP -override-limits -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
@@ -92,6 +92,7 @@
# uncomment this to use Metis on MareNostrum in Barcelona
#LIB = /home/hpce08/hpce08548/utils/metis-4.0/libmetis.a
+
LINK = $(F90)
OBJS_MESHFEM2D = $O/part_unstruct.o $O/meshfem2D.o $O/read_value_parameters.o $O/spline_routines.o
@@ -100,15 +101,16 @@
$O/compute_forces_acoustic.o $O/compute_forces_elastic.o\
$O/compute_forces_solid.o $O/compute_forces_fluid.o\
$O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivation_matrices.o\
- $O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o\
+ $O/plotpost.o $O/locate_receivers.o $O/locate_source_force.o $O/compute_gradient_attenuation.o $O/setup_sources_receivers.o\
$O/specfem2D.o $O/write_seismograms.o $O/define_external_model.o $O/createnum_fast.o $O/createnum_slow.o\
$O/define_shape_functions.o $O/attenuation_model.o $O/create_color_image.o $O/compute_vector_field.o $O/compute_pressure.o\
$O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/netlib_specfun_erf.o\
$O/construct_acoustic_surface.o $O/assemble_MPI.o $O/compute_energy.o $O/compute_curl_one_element.o\
$O/attenuation_compute_param.o $O/compute_Bielak_conditions.o $O/paco_beyond_critical.o\
- $O/paco_convolve_fft.o $O/is_in_convex_quadrilateral.o $O/get_perm_cuthill_mckee.o
+ $O/paco_convolve_fft.o $O/is_in_convex_quadrilateral.o $O/get_perm_cuthill_mckee.o\
+ $O/read_external_model.o $O/invert_mass_matrix.o
-default: clean meshfem2D specfem2D convolve_source_timefunction
+default: clean meshfem2D specfem2D convolve_source_timefunction
all: default
@@ -267,3 +269,11 @@
$O/get_perm_cuthill_mckee.o: get_perm_cuthill_mckee.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/get_perm_cuthill_mckee.o get_perm_cuthill_mckee.f90
+$O/read_external_model.o: read_external_model.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/read_external_model.o read_external_model.f90
+
+$O/setup_sources_receivers.o: setup_sources_receivers.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/setup_sources_receivers.o setup_sources_receivers.f90
+
+$O/invert_mass_matrix.o: invert_mass_matrix.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/invert_mass_matrix.o invert_mass_matrix.f90
Modified: seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt
===================================================================
--- seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/README_MANUAL.txt 2010-03-12 20:23:22 UTC (rev 16405)
@@ -251,7 +251,7 @@
Three types of models:
I: (model_number,1,rho,Vp,Vs,0,0,Qp,Qs,0,0,0,0,0,0), for isotropic elastic/acoustic
material
-or II: (model_number,2,rho,c11,c13,c33,c44,Qp,Qs,0,0,0,0,0,0), for anisotropic material
+or II: (model_number,2,rho,Vp,Vs,c11,c13,c15,c33,c35,c55,Qp,Qs,0,0), for anisotropic material
or III: (model_number,3,rhos,rhof,phi,c,kxx,kxz,kzz,Ks,Kf,Kfr,etaf,mufr,Qs),
for isotropic poroelastic material
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -45,7 +45,7 @@
subroutine checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,ibool,kmato,coord,npoin, &
vpImin,vpImax,vpIImin,vpIImax,assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat, &
f0,t0,initialfield,time_function_type,coorg,xinterp,zinterp,shapeint,knods,simulation_title, &
- npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,myrank,nproc,NSOURCE,poroelastic)
+ npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic,myrank,nproc,NSOURCE,poroelastic)
! check the mesh, stability and number of points per wavelength
@@ -94,7 +94,7 @@
lambdaSmin,lambdaSmax
double precision deltat,distance_1,distance_2,distance_3,distance_4
double precision, dimension(NSOURCE) :: f0,t0
- logical assign_external_model,initialfield,any_elastic,any_poroelastic
+ logical assign_external_model,initialfield,any_elastic,any_poroelastic,all_anisotropic
! for the stability condition
! maximum polynomial degree for which we can compute the stability condition
@@ -141,6 +141,7 @@
! title of the plot
character(len=60) simulation_title
+
if(UPPER_LIMIT_DISPLAY > nspec) stop 'cannot have UPPER_LIMIT_DISPLAY > nspec in checkgrid.F90'
#ifndef USE_MPI
@@ -1359,9 +1360,8 @@
! Bisque
red(236) = 1.00000000000000
green(236) = 0.894117647058824
- blue(236) = 0.768627450980392
+ blue(236) = 0.768627450980392
-
!---- compute parameters for the spectral elements
vpImin = HUGEVAL
@@ -1578,28 +1578,29 @@
#endif
if ( myrank == 0 ) then
- write(IOUT,*)
- write(IOUT,*) '********'
- write(IOUT,*) 'Model: P (or PI) velocity min,max = ',vpImin,vpImax
- write(IOUT,*) 'Model: PII velocity min,max = ',vpIImin,vpIImax
- write(IOUT,*) 'Model: S velocity min,max = ',vsmin,vsmax
- write(IOUT,*) 'Model: density min,max = ',densmin,densmax
- write(IOUT,*) '********'
- write(IOUT,*)
+ if(.not. all_anisotropic) then
+ write(IOUT,*)
+ write(IOUT,*) '********'
+ write(IOUT,*) 'Model: P (or PI) velocity min,max = ',vpImin,vpImax
+ write(IOUT,*) 'Model: PII velocity min,max = ',vpIImin,vpIImax
+ write(IOUT,*) 'Model: S velocity min,max = ',vsmin,vsmax
+ write(IOUT,*) 'Model: density min,max = ',densmin,densmax
+ write(IOUT,*) '********'
+ write(IOUT,*)
- write(IOUT,*)
- write(IOUT,*) '*********************************************'
- write(IOUT,*) '*** Verification of simulation parameters ***'
- write(IOUT,*) '*********************************************'
- write(IOUT,*)
- write(IOUT,*) '*** Max grid size = ',distance_max
- write(IOUT,*) '*** Min grid size = ',distance_min
- write(IOUT,*) '*** Max/min ratio = ',distance_max/distance_min
- write(IOUT,*)
- write(IOUT,*) '*** Max stability for P wave velocity = ',courant_stability_number_max
- write(IOUT,*)
+ write(IOUT,*)
+ write(IOUT,*) '*********************************************'
+ write(IOUT,*) '*** Verification of simulation parameters ***'
+ write(IOUT,*) '*********************************************'
+ write(IOUT,*)
+ write(IOUT,*) '*** Max grid size = ',distance_max
+ write(IOUT,*) '*** Min grid size = ',distance_min
+ write(IOUT,*) '*** Max/min ratio = ',distance_max/distance_min
+ write(IOUT,*)
+ write(IOUT,*) '*** Max stability for P wave velocity = ',courant_stability_number_max
+ write(IOUT,*)
+ end if
-
! only if time source is not a Dirac or Heaviside (otherwise maximum frequency of spectrum undefined)
! and if source is not an initial field, for the same reason
if(.not. initialfield) then
@@ -1690,7 +1691,7 @@
!
!---- write PostScript header
!
- write(24,10) simulation_title
+ write(24,*) simulation_title
write(24,*) '/CM {28.5 mul} def'
write(24,*) '/LR {rlineto} def'
write(24,*) '/LT {lineto} def'
@@ -1805,7 +1806,7 @@
zinterp(i,j) = zinterp(i,j) + shapeint(in,i,j)*coorg(2,nnum)
enddo
enddo
- enddo
+ enddo
is = 1
ir = 1
@@ -2681,7 +2682,7 @@
endif
else
x1 = 0.5d0
- endif
+ endif
! rescale to avoid very dark gray levels
x1 = x1*0.7 + 0.2
@@ -3003,3 +3004,4 @@
end subroutine checkgrid
+
Modified: seismo/2D/SPECFEM2D/trunk/compute_energy.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/compute_energy.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -47,9 +47,9 @@
xix,xiz,gammax,gammaz,jacobian,ibool,elastic,poroelastic,hprime_xx,hprime_zz, &
nspec,npoin,assign_external_model,it,deltat,t0,kmato,elastcoef,density, &
porosity,tortuosity, &
- vpext,vsext,rhoext,wxgll,wzgll,numat, &
+ vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,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,N_SLS)
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute kinetic and potential energy in the solid (acoustic elements are excluded)
@@ -69,7 +69,7 @@
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_acoustic,potential_dot_dot_acoustic
- logical :: TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
+ logical :: TURN_ATTENUATION_ON
integer :: nspec,npoin,numat
@@ -78,7 +78,7 @@
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
- logical, dimension(nspec) :: elastic,poroelastic
+ logical, dimension(nspec) :: elastic,poroelastic,anisotropic
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
@@ -88,8 +88,10 @@
double precision, dimension(2,numat) :: density
double precision, dimension(numat) :: porosity,tortuosity
+ double precision, dimension(6,numat) :: anisotropy
double precision, dimension(4,3,numat) :: elastcoef
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displ_elastic,veloc_elastic
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,velocs_poroelastic
@@ -324,8 +326,9 @@
! 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,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
+ numat,kmato,elastcoef,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
+ ispec,e1,e11, &
+ TURN_ATTENUATION_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 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -42,25 +42,26 @@
!
!========================================================================
- subroutine compute_forces_elastic(p_sv,npoin,nspec,myrank,nelemabs,numat, &
- ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
- source_type,it,NSTEP,anyabs,assign_external_model, &
- initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
- deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
- accel_elastic,veloc_elastic,displ_elastic,b_accel_elastic,b_displ_elastic, &
- density,elastcoef,xix,xiz,gammax,gammaz, &
- jacobian,vpext,vsext,rhoext,source_time_function,sourcearray,adj_sourcearrays,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,N_SLS, &
- deltat,coord,add_Bielak_conditions, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset,f0, &
- v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,&
- nleft,nright,nbot,over_critical_angle,NSOURCE,nrec,isolver,save_forward,b_absorb_elastic_left,&
- b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top,nspec_xmin,nspec_xmax,&
- nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,mu_k,kappa_k)
+subroutine compute_forces_elastic(p_sv,npoin,nspec,myrank,nelemabs,numat, &
+ ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
+ source_type,it,NSTEP,anyabs,assign_external_model, &
+ initialfield,TURN_ATTENUATION_ON,angleforce,deltatcube, &
+ deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
+ accel_elastic,veloc_elastic,displ_elastic,b_accel_elastic,b_displ_elastic, &
+ density,elastcoef,xix,xiz,gammax,gammaz, &
+ jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
+ source_time_function,sourcearray,adj_sourcearrays,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,N_SLS, &
+ deltat,coord,add_Bielak_conditions, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset,f0, &
+ v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,&
+ nleft,nright,nbot,over_critical_angle,NSOURCE,nrec,isolver,save_forward,b_absorb_elastic_left,&
+ b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top,nspec_xmin,nspec_xmax,&
+ nspec_zmin,nspec_zmax,ib_xmin,ib_xmax,ib_zmin,ib_zmax,mu_k,kappa_k)
-! compute forces for the elastic elements
+ ! compute forces for the elastic elements
implicit none
@@ -79,7 +80,7 @@
integer, dimension(nspec_zmin) :: ib_zmin
integer, dimension(nspec_zmax) :: ib_zmax
- logical :: anyabs,assign_external_model,initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,add_Bielak_conditions
+ logical :: anyabs,assign_external_model,initialfield,TURN_ATTENUATION_ON,add_Bielak_conditions
logical :: save_forward
@@ -90,14 +91,17 @@
integer, dimension(nspec) :: kmato
integer, dimension(nelemabs) :: numabs
- logical, dimension(nspec) :: elastic
+ logical, dimension(nspec) :: elastic,anisotropic
logical, dimension(4,nelemabs) :: codeabs
real(kind=CUSTOM_REAL), dimension(3,npoin) :: accel_elastic,veloc_elastic,displ_elastic
double precision, dimension(2,numat) :: density
double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(6,numat) :: anisotropy
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,vsext,rhoext
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
+
real(kind=CUSTOM_REAL), dimension(NSOURCE,NSTEP) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
@@ -117,24 +121,24 @@
integer :: i_sls
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
- dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
+ dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
-! derivatives of Lagrange polynomials
+ ! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
-! Gauss-Lobatto-Legendre weights
+ ! Gauss-Lobatto-Legendre weights
real(kind=CUSTOM_REAL), dimension(NGLLX) :: wxgll
real(kind=CUSTOM_REAL), dimension(NGLLZ) :: wzgll
-!---
-!--- local variables
-!---
+ !---
+ !--- local variables
+ !---
integer :: ispec,i,j,k,iglob,ispecabs,ibegin,iend,irec,irec_local
-! spatial derivatives
+ ! spatial derivatives
real(kind=CUSTOM_REAL) :: dux_dxi,dux_dgamma,duy_dxi,duy_dgamma,duz_dxi,duz_dgamma
real(kind=CUSTOM_REAL) :: dux_dxl,duy_dxl,duz_dxl,dux_dzl,duy_dzl,duz_dzl
real(kind=CUSTOM_REAL) :: b_dux_dxi,b_dux_dgamma,b_duy_dxi,b_duy_dgamma,b_duz_dxi,b_duz_dgamma
@@ -148,22 +152,25 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: tempx1,tempx2,tempy1,tempy2,tempz1,tempz2
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: b_tempx1,b_tempx2,b_tempy1,b_tempy2,b_tempz1,b_tempz2
-! Jacobian matrix and determinant
+ ! Jacobian matrix and determinant
real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl,jacobianl
-! material properties of the elastic medium
+ ! material properties of the elastic medium
real(kind=CUSTOM_REAL) :: mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed,kappal,cpl,csl,rhol, &
- lambdal_unrelaxed,mul_unrelaxed,lambdalplus2mul_unrelaxed
+ lambdal_unrelaxed,mul_unrelaxed,lambdalplus2mul_unrelaxed
-! for attenuation
+ ! for attenuation
real(kind=CUSTOM_REAL) :: Un,Unp1,tauinv,Sn,Snp1,theta_n,theta_np1,tauinvsquare,tauinvcube,tauinvUn
-! for analytical initial plane wave for Bielak's conditions
+ ! for anisotropy
+ double precision :: c11,c15,c13,c33,c35,c55
+
+ ! for analytical initial plane wave for Bielak's conditions
double precision :: veloc_horiz,veloc_vert,dxUx,dzUx,dxUz,dzUz,traction_x_t0,traction_z_t0,deltat
double precision, dimension(NDIM,npoin), intent(in) :: coord
double precision x0_source, z0_source, angleforce_refl, c_inc, c_refl, time_offset, f0
double precision, dimension(NDIM) :: A_plane, B_plane, C_plane
-!over critical angle
+ !over critical angle
logical :: over_critical_angle
integer :: nleft, nright, nbot
double precision, dimension(nleft) :: v0x_left,v0z_left,t0x_left,t0z_left
@@ -173,785 +180,809 @@
integer :: ifirstelem,ilastelem
-! compute Grad(displ_elastic) at time step n for attenuation
+ ! compute Grad(displ_elastic) at time step n for attenuation
if(TURN_ATTENUATION_ON) then
- call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
- dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,npoin)
+ call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
+ dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,npoin)
endif
ifirstelem = 1
ilastelem = nspec
-! loop over spectral elements
+ ! loop over spectral elements
do ispec = ifirstelem,ilastelem
- tempx1(:,:) = ZERO
- tempy1(:,:) = ZERO
- tempz1(:,:) = ZERO
- tempx2(:,:) = ZERO
- tempy2(:,:) = ZERO
- tempz2(:,:) = ZERO
- if(isolver ==2)then
- b_tempx1(:,:) = ZERO
- b_tempy1(:,:) = ZERO
- b_tempz1(:,:) = ZERO
- b_tempx2(:,:) = ZERO
- b_tempy2(:,:) = ZERO
- b_tempz2(:,:) = ZERO
- endif
+ tempx1(:,:) = ZERO
+ tempy1(:,:) = ZERO
+ tempz1(:,:) = ZERO
+ tempx2(:,:) = ZERO
+ tempy2(:,:) = ZERO
+ tempz2(:,:) = ZERO
+ if(isolver ==2)then
+ b_tempx1(:,:) = ZERO
+ b_tempy1(:,:) = ZERO
+ b_tempz1(:,:) = ZERO
+ b_tempx2(:,:) = ZERO
+ b_tempy2(:,:) = ZERO
+ b_tempz2(:,:) = ZERO
+ endif
-!---
-!--- elastic spectral element
-!---
- if(elastic(ispec)) then
+ !---
+ !--- elastic spectral element
+ !---
+ if(elastic(ispec)) then
-! get relaxed elastic parameters of current spectral element
- lambdal_relaxed = elastcoef(1,1,kmato(ispec))
- mul_relaxed = elastcoef(2,1,kmato(ispec))
- lambdalplus2mul_relaxed = elastcoef(3,1,kmato(ispec))
+ ! get relaxed elastic parameters of current spectral element
+ lambdal_relaxed = elastcoef(1,1,kmato(ispec))
+ mul_relaxed = elastcoef(2,1,kmato(ispec))
+ lambdalplus2mul_relaxed = elastcoef(3,1,kmato(ispec))
-! first double loop over GLL points to compute and store gradients
- do j = 1,NGLLZ
- do i = 1,NGLLX
+ ! first double loop over GLL points to compute and store gradients
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
-!--- if external medium, get elastic parameters of current grid point
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- csl = vsext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- mul_relaxed = rhol*csl*csl
- lambdal_relaxed = rhol*cpl*cpl - TWO*mul_relaxed
- lambdalplus2mul_relaxed = lambdal_relaxed + TWO*mul_relaxed
- endif
+ !--- if external medium, get elastic parameters of current grid point
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ mul_relaxed = rhol*csl*csl
+ lambdal_relaxed = rhol*cpl*cpl - TWO*mul_relaxed
+ lambdalplus2mul_relaxed = lambdal_relaxed + TWO*mul_relaxed
+ endif
-! derivative along x and along z
- dux_dxi = ZERO
- duy_dxi = ZERO
- duz_dxi = ZERO
+ ! derivative along x and along z
+ dux_dxi = ZERO
+ duy_dxi = ZERO
+ duz_dxi = ZERO
- dux_dgamma = ZERO
- duy_dgamma = ZERO
- duz_dgamma = ZERO
+ dux_dgamma = ZERO
+ duy_dgamma = ZERO
+ duz_dgamma = ZERO
- if(isolver == 2) then ! Adjoint calculation, backward wavefield
- b_dux_dxi = ZERO
- b_duy_dxi = ZERO
- b_duz_dxi = ZERO
+ if(isolver == 2) then ! Adjoint calculation, backward wavefield
+ b_dux_dxi = ZERO
+ b_duy_dxi = ZERO
+ b_duz_dxi = ZERO
- b_dux_dgamma = ZERO
- b_duy_dgamma = ZERO
- b_duz_dgamma = ZERO
- endif
+ b_dux_dgamma = ZERO
+ b_duy_dgamma = ZERO
+ b_duz_dgamma = ZERO
+ endif
-! first double loop over GLL points to compute and store gradients
-! we can merge the two loops because NGLLX == NGLLZ
- do k = 1,NGLLX
- dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
- duy_dxi = duy_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
- duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
- dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- duy_dgamma = duy_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
- duz_dgamma = duz_dgamma + displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
+ duy_dxi = duy_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ duz_dxi = duz_dxi + displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
+ dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
+ duy_dgamma = duy_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
- if(isolver == 2) then ! Adjoint calculation, backward wavefield
- b_dux_dxi = b_dux_dxi + b_displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
- b_duy_dxi = b_duy_dxi + b_displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
- b_duz_dxi = b_duz_dxi + b_displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
- b_dux_dgamma = b_dux_dgamma + b_displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
- b_duy_dgamma = b_duy_dgamma + b_displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
- b_duz_dgamma = b_duz_dgamma + b_displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
- endif
- enddo
+ if(isolver == 2) then ! Adjoint calculation, backward wavefield
+ b_dux_dxi = b_dux_dxi + b_displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
+ b_duy_dxi = b_duy_dxi + b_displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ b_duz_dxi = b_duz_dxi + b_displ_elastic(3,ibool(k,j,ispec))*hprime_xx(i,k)
+ b_dux_dgamma = b_dux_dgamma + b_displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
+ b_duy_dgamma = b_duy_dgamma + b_displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ b_duz_dgamma = b_duz_dgamma + b_displ_elastic(3,ibool(i,k,ispec))*hprime_zz(j,k)
+ endif
+ enddo
- xixl = xix(i,j,ispec)
- xizl = xiz(i,j,ispec)
- gammaxl = gammax(i,j,ispec)
- gammazl = gammaz(i,j,ispec)
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
-! derivatives of displacement
- dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
- dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+ ! derivatives of displacement
+ dux_dxl = dux_dxi*xixl + dux_dgamma*gammaxl
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
- duy_dxl = duy_dxi*xixl + duy_dgamma*gammaxl
- duy_dzl = duy_dxi*xizl + duy_dgamma*gammazl
+ duy_dxl = duy_dxi*xixl + duy_dgamma*gammaxl
+ duy_dzl = duy_dxi*xizl + duy_dgamma*gammazl
- duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
- duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
+ duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
+ duz_dzl = duz_dxi*xizl + duz_dgamma*gammazl
- if(isolver == 2) then ! Adjoint calculation, backward wavefield
- b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
- b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
+ if(isolver == 2) then ! Adjoint calculation, backward wavefield
+ b_dux_dxl = b_dux_dxi*xixl + b_dux_dgamma*gammaxl
+ b_dux_dzl = b_dux_dxi*xizl + b_dux_dgamma*gammazl
- b_duy_dxl = b_duy_dxi*xixl + b_duy_dgamma*gammaxl
- b_duy_dzl = b_duy_dxi*xizl + b_duy_dgamma*gammazl
+ b_duy_dxl = b_duy_dxi*xixl + b_duy_dgamma*gammaxl
+ b_duy_dzl = b_duy_dxi*xizl + b_duy_dgamma*gammazl
- b_duz_dxl = b_duz_dxi*xixl + b_duz_dgamma*gammaxl
- b_duz_dzl = b_duz_dxi*xizl + b_duz_dgamma*gammazl
- endif
+ b_duz_dxl = b_duz_dxi*xixl + b_duz_dgamma*gammaxl
+ b_duz_dzl = b_duz_dxi*xizl + b_duz_dgamma*gammazl
+ endif
-! compute stress tensor (include attenuation or anisotropy if needed)
+ ! compute stress tensor (include attenuation or anisotropy if needed)
- if(TURN_ATTENUATION_ON) then
+ if(TURN_ATTENUATION_ON) then
-! attenuation is implemented following the memory variable formulation of
-! J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
-! vol. 58(1), p. 110-120 (1993). More details can be found in
-! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a linear
-! viscoelastic medium, Geophysical Journal International, vol. 95, p. 597-611 (1988).
+ ! attenuation is implemented following the memory variable formulation of
+ ! J. M. Carcione, Seismic modeling in viscoelastic media, Geophysics,
+ ! vol. 58(1), p. 110-120 (1993). More details can be found in
+ ! J. M. Carcione, D. Kosloff and R. Kosloff, Wave propagation simulation in a linear
+ ! 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(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 unrelaxed elastic coefficients from formulas in Carcione 1993 page 111
+ 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)
- sigma_xx = lambdalplus2mul_unrelaxed*dux_dxl + lambdal_unrelaxed*duz_dzl
- sigma_xz = mul_unrelaxed*(duz_dxl + dux_dzl)
- sigma_zz = lambdalplus2mul_unrelaxed*duz_dzl + lambdal_unrelaxed*dux_dxl
+ ! compute the stress using the unrelaxed Lame parameters (Carcione 1993, page 111)
+ sigma_xx = lambdalplus2mul_unrelaxed*dux_dxl + lambdal_unrelaxed*duz_dzl
+ sigma_xz = mul_unrelaxed*(duz_dxl + dux_dzl)
+ sigma_zz = lambdalplus2mul_unrelaxed*duz_dzl + lambdal_unrelaxed*dux_dxl
-! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
-! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
- e1_sum = 0._CUSTOM_REAL
- e11_sum = 0._CUSTOM_REAL
- e13_sum = 0._CUSTOM_REAL
+ ! add the memory variables using the relaxed parameters (Carcione 1993, page 111)
+ ! beware: there is a bug in Carcione's equation (2c) for sigma_zz, we fixed it in the code below
+ e1_sum = 0._CUSTOM_REAL
+ e11_sum = 0._CUSTOM_REAL
+ e13_sum = 0._CUSTOM_REAL
- do i_sls = 1,N_SLS
- e1_sum = e1_sum + e1(i,j,ispec,i_sls)
- e11_sum = e11_sum + e11(i,j,ispec,i_sls)
- e13_sum = e13_sum + e13(i,j,ispec,i_sls)
- enddo
+ do i_sls = 1,N_SLS
+ e1_sum = e1_sum + e1(i,j,ispec,i_sls)
+ e11_sum = e11_sum + e11(i,j,ispec,i_sls)
+ e13_sum = e13_sum + e13(i,j,ispec,i_sls)
+ enddo
- sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum + TWO * mul_relaxed * e11_sum
- sigma_xz = sigma_xz + mul_relaxed * e13_sum
- sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum - TWO * mul_relaxed * e11_sum
+ sigma_xx = sigma_xx + (lambdal_relaxed + mul_relaxed) * e1_sum + TWO * mul_relaxed * e11_sum
+ sigma_xz = sigma_xz + mul_relaxed * e13_sum
+ sigma_zz = sigma_zz + (lambdal_relaxed + mul_relaxed) * e1_sum - TWO * mul_relaxed * e11_sum
- else
+ else
-! no attenuation
- sigma_xx = lambdalplus2mul_relaxed*dux_dxl + lambdal_relaxed*duz_dzl
- sigma_xy = mul_relaxed*duy_dxl
- sigma_xz = mul_relaxed*(duz_dxl + dux_dzl)
- sigma_zy = mul_relaxed*duy_dzl
- sigma_zz = lambdalplus2mul_relaxed*duz_dzl + lambdal_relaxed*dux_dxl
+ ! no attenuation
+ sigma_xx = lambdalplus2mul_relaxed*dux_dxl + lambdal_relaxed*duz_dzl
+ sigma_xy = mul_relaxed*duy_dxl
+ sigma_xz = mul_relaxed*(duz_dxl + dux_dzl)
+ sigma_zy = mul_relaxed*duy_dzl
+ sigma_zz = lambdalplus2mul_relaxed*duz_dzl + lambdal_relaxed*dux_dxl
- if(isolver == 2) then ! Adjoint calculation, backward wavefield
- b_sigma_xx = lambdalplus2mul_relaxed*b_dux_dxl + lambdal_relaxed*b_duz_dzl
- b_sigma_xy = mul_relaxed*b_duy_dxl
- b_sigma_xz = mul_relaxed*(b_duz_dxl + b_dux_dzl)
- b_sigma_zy = mul_relaxed*b_duy_dzl
- b_sigma_zz = lambdalplus2mul_relaxed*b_duz_dzl + lambdal_relaxed*b_dux_dxl
- endif
+ if(isolver == 2) then ! Adjoint calculation, backward wavefield
+ b_sigma_xx = lambdalplus2mul_relaxed*b_dux_dxl + lambdal_relaxed*b_duz_dzl
+ b_sigma_xy = mul_relaxed*b_duy_dxl
+ b_sigma_xz = mul_relaxed*(b_duz_dxl + b_dux_dzl)
+ b_sigma_zy = mul_relaxed*b_duy_dzl
+ b_sigma_zz = lambdalplus2mul_relaxed*b_duz_dzl + lambdal_relaxed*b_dux_dxl
+ endif
- endif
+ endif
-! full anisotropy
- if(TURN_ANISOTROPY_ON) then
+ ! full anisotropy
+ if(anisotropic(ispec)) then
+ if(assign_external_model) then
+ c11 = c11ext(i,j,ispec)
+ c13 = c13ext(i,j,ispec)
+ c15 = c15ext(i,j,ispec)
+ c33 = c33ext(i,j,ispec)
+ c35 = c35ext(i,j,ispec)
+ c55 = c55ext(i,j,ispec)
+ else
+ c11 = anisotropy(1,kmato(ispec))
+ c13 = anisotropy(2,kmato(ispec))
+ c15 = anisotropy(3,kmato(ispec))
+ c33 = anisotropy(4,kmato(ispec))
+ c35 = anisotropy(5,kmato(ispec))
+ c55 = anisotropy(6,kmato(ispec))
+ end if
+
+ ! implement anisotropy in 2D
+ sigma_xx = c11*dux_dxl + c15*(duz_dxl + dux_dzl) + c13*duz_dzl
+ sigma_zz = c13*dux_dxl + c35*(duz_dxl + dux_dzl) + c33*duz_dzl
+ sigma_xz = c15*dux_dxl + c55*(duz_dxl + dux_dzl) + c35*duz_dzl
-! implement anisotropy in 2D
- sigma_xx = c11val*dux_dxl + c15val*(duz_dxl + dux_dzl) + c13val*duz_dzl
- sigma_zz = c13val*dux_dxl + c35val*(duz_dxl + dux_dzl) + c33val*duz_dzl
- sigma_xz = c15val*dux_dxl + c55val*(duz_dxl + dux_dzl) + c35val*duz_dzl
+ endif
- endif
+ ! Pre-kernels calculation
+ if(isolver == 2) then
+ iglob = ibool(i,j,ispec)
+ if(p_sv)then !P-SV waves
+ dsxx = dux_dxl
+ dsxz = HALF * (duz_dxl + dux_dzl)
+ dszz = duz_dzl
-! Pre-kernels calculation
- if(isolver == 2) then
- iglob = ibool(i,j,ispec)
- if(p_sv)then !P-SV waves
- dsxx = dux_dxl
- dsxz = HALF * (duz_dxl + dux_dzl)
- dszz = duz_dzl
+ b_dsxx = b_dux_dxl
+ b_dsxz = HALF * (b_duz_dxl + b_dux_dzl)
+ b_dszz = b_duz_dzl
- b_dsxx = b_dux_dxl
- b_dsxz = HALF * (b_duz_dxl + b_dux_dzl)
- b_dszz = b_duz_dzl
+ kappa_k(iglob) = (dux_dxl + duz_dzl) * (b_dux_dxl + b_duz_dzl)
+ mu_k(iglob) = dsxx * b_dsxx + dszz * b_dszz + &
+ 2._CUSTOM_REAL * dsxz * b_dsxz - 1._CUSTOM_REAL/3._CUSTOM_REAL * kappa_k(iglob)
+ else !SH (membrane) waves
+ mu_k(iglob) = duy_dxl * b_duy_dxl + duy_dzl * b_duy_dzl
+ endif
+ endif
- kappa_k(iglob) = (dux_dxl + duz_dzl) * (b_dux_dxl + b_duz_dzl)
- mu_k(iglob) = dsxx * b_dsxx + dszz * b_dszz + &
- 2._CUSTOM_REAL * dsxz * b_dsxz - 1._CUSTOM_REAL/3._CUSTOM_REAL * kappa_k(iglob)
- else !SH (membrane) waves
- mu_k(iglob) = duy_dxl * b_duy_dxl + duy_dzl * b_duy_dzl
- endif
- endif
+ jacobianl = jacobian(i,j,ispec)
- jacobianl = jacobian(i,j,ispec)
+ ! weak formulation term based on stress tensor (non-symmetric form)
+ ! also add GLL integration weights
+ tempx1(i,j) = wzgll(j)*jacobianl*(sigma_xx*xixl+sigma_xz*xizl)
+ !tempy1(i,j) = wzgll(j)*jacobianl*(sigma_xy*xixl+sigma_zy*xizl)
+ tempz1(i,j) = wzgll(j)*jacobianl*(sigma_xz*xixl+sigma_zz*xizl)
-! weak formulation term based on stress tensor (non-symmetric form)
-! also add GLL integration weights
- tempx1(i,j) = wzgll(j)*jacobianl*(sigma_xx*xixl+sigma_xz*xizl)
- tempy1(i,j) = wzgll(j)*jacobianl*(sigma_xy*xixl+sigma_zy*xizl)
- tempz1(i,j) = wzgll(j)*jacobianl*(sigma_xz*xixl+sigma_zz*xizl)
+ tempx2(i,j) = wxgll(i)*jacobianl*(sigma_xx*gammaxl+sigma_xz*gammazl)
+ !tempy2(i,j) = wxgll(i)*jacobianl*(sigma_xy*gammaxl+sigma_zy*gammazl)
+ tempz2(i,j) = wxgll(i)*jacobianl*(sigma_xz*gammaxl+sigma_zz*gammazl)
- tempx2(i,j) = wxgll(i)*jacobianl*(sigma_xx*gammaxl+sigma_xz*gammazl)
- tempy2(i,j) = wxgll(i)*jacobianl*(sigma_xy*gammaxl+sigma_zy*gammazl)
- tempz2(i,j) = wxgll(i)*jacobianl*(sigma_xz*gammaxl+sigma_zz*gammazl)
+ if(isolver == 2) then ! Adjoint calculation, backward wavefield
+ b_tempx1(i,j) = wzgll(j)*jacobianl*(b_sigma_xx*xixl+b_sigma_xz*xizl)
+ !b_tempy1(i,j) = wzgll(j)*jacobianl*(b_sigma_xy*xixl+b_sigma_zy*xizl)
+ b_tempz1(i,j) = wzgll(j)*jacobianl*(b_sigma_xz*xixl+b_sigma_zz*xizl)
- if(isolver == 2) then ! Adjoint calculation, backward wavefield
- b_tempx1(i,j) = wzgll(j)*jacobianl*(b_sigma_xx*xixl+b_sigma_xz*xizl)
- b_tempy1(i,j) = wzgll(j)*jacobianl*(b_sigma_xy*xixl+b_sigma_zy*xizl)
- b_tempz1(i,j) = wzgll(j)*jacobianl*(b_sigma_xz*xixl+b_sigma_zz*xizl)
+ b_tempx2(i,j) = wxgll(i)*jacobianl*(b_sigma_xx*gammaxl+b_sigma_xz*gammazl)
+ !b_tempy2(i,j) = wxgll(i)*jacobianl*(b_sigma_xy*gammaxl+b_sigma_zy*gammazl)
+ b_tempz2(i,j) = wxgll(i)*jacobianl*(b_sigma_xz*gammaxl+b_sigma_zz*gammazl)
+ endif
- b_tempx2(i,j) = wxgll(i)*jacobianl*(b_sigma_xx*gammaxl+b_sigma_xz*gammazl)
- b_tempy2(i,j) = wxgll(i)*jacobianl*(b_sigma_xy*gammaxl+b_sigma_zy*gammazl)
- b_tempz2(i,j) = wxgll(i)*jacobianl*(b_sigma_xz*gammaxl+b_sigma_zz*gammazl)
- endif
-
+ enddo
enddo
- enddo
-!
-! second double-loop over GLL to compute all the terms
-!
- do j = 1,NGLLZ
- do i = 1,NGLLX
+ !
+ ! second double-loop over GLL to compute all the terms
+ !
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
- iglob = ibool(i,j,ispec)
+ iglob = ibool(i,j,ispec)
-! along x direction and z direction
-! and assemble the contributions
-! we can merge the two loops because NGLLX == NGLLZ
- do k = 1,NGLLX
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tempx1(k,j)*hprimewgll_xx(k,i) + tempx2(i,k)*hprimewgll_zz(k,j))
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tempy1(k,j)*hprimewgll_xx(k,i) + tempy2(i,k)*hprimewgll_zz(k,j))
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tempz1(k,j)*hprimewgll_xx(k,i) + tempz2(i,k)*hprimewgll_zz(k,j))
+ ! along x direction and z direction
+ ! and assemble the contributions
+ ! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tempx1(k,j)*hprimewgll_xx(k,i) + tempx2(i,k)*hprimewgll_zz(k,j))
+ ! accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tempy1(k,j)*hprimewgll_xx(k,i) + tempy2(i,k)*hprimewgll_zz(k,j))
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tempz1(k,j)*hprimewgll_xx(k,i) + tempz2(i,k)*hprimewgll_zz(k,j))
- if(isolver == 2) then ! Adjoint calculation, backward wavefield
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
+ if(isolver == 2) then ! Adjoint calculation, backward wavefield
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
(b_tempx1(k,j)*hprimewgll_xx(k,i) + b_tempx2(i,k)*hprimewgll_zz(k,j))
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
(b_tempy1(k,j)*hprimewgll_xx(k,i) + b_tempy2(i,k)*hprimewgll_zz(k,j))
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
(b_tempz1(k,j)*hprimewgll_xx(k,i) + b_tempz2(i,k)*hprimewgll_zz(k,j))
- endif
- enddo
+ endif
+ enddo
- enddo ! second loop over the GLL points
- enddo
+ enddo ! second loop over the GLL points
+ enddo
- endif ! end of test if elastic element
+ endif ! end of test if elastic element
- enddo ! end of loop over all spectral elements
+ enddo ! end of loop over all spectral elements
-!
-!--- absorbing boundaries
-!
+ !
+ !--- absorbing boundaries
+ !
if(anyabs) then
- count_left=1
- count_right=1
- count_bot=1
+ count_left=1
+ count_right=1
+ count_bot=1
- do ispecabs = 1,nelemabs
+ do ispecabs = 1,nelemabs
- ispec = numabs(ispecabs)
+ ispec = numabs(ispecabs)
-! get elastic parameters of current spectral element
- lambdal_relaxed = elastcoef(1,1,kmato(ispec))
- mul_relaxed = elastcoef(2,1,kmato(ispec))
- rhol = density(1,kmato(ispec))
- kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
- cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_relaxed/3._CUSTOM_REAL)/rhol)
- csl = sqrt(mul_relaxed/rhol)
+ ! get elastic parameters of current spectral element
+ lambdal_relaxed = elastcoef(1,1,kmato(ispec))
+ mul_relaxed = elastcoef(2,1,kmato(ispec))
+ rhol = density(1,kmato(ispec))
+ kappal = lambdal_relaxed + TWO*mul_relaxed/3._CUSTOM_REAL
+ cpl = sqrt((kappal + 4._CUSTOM_REAL*mul_relaxed/3._CUSTOM_REAL)/rhol)
+ csl = sqrt(mul_relaxed/rhol)
-!--- left absorbing boundary
- if(codeabs(ILEFT,ispecabs)) then
+ !--- left absorbing boundary
+ if(codeabs(ILEFT,ispecabs)) then
- i = 1
+ i = 1
- do j = 1,NGLLZ
+ do j = 1,NGLLZ
- iglob = ibool(i,j,ispec)
+ iglob = ibool(i,j,ispec)
-! for analytical initial plane wave for Bielak's conditions
-! left or right edge, horizontal normal vector
- if(add_Bielak_conditions .and. initialfield) then
- if (.not.over_critical_angle) then
- call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
- c_inc, c_refl, time_offset,f0)
- traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
- traction_z_t0 = mul_relaxed*(dxUz + dzUx)
- else
- veloc_horiz=v0x_left(count_left)
- veloc_vert=v0z_left(count_left)
- traction_x_t0=t0x_left(count_left)
- traction_z_t0=t0z_left(count_left)
- count_left=count_left+1
- end if
- else
- veloc_horiz = 0
- veloc_vert = 0
- traction_x_t0 = 0
- traction_z_t0 = 0
- endif
+ ! for analytical initial plane wave for Bielak's conditions
+ ! left or right edge, horizontal normal vector
+ if(add_Bielak_conditions .and. initialfield) then
+ if (.not.over_critical_angle) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ c_inc, c_refl, time_offset,f0)
+ traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
+ traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ else
+ veloc_horiz=v0x_left(count_left)
+ veloc_vert=v0z_left(count_left)
+ traction_x_t0=t0x_left(count_left)
+ traction_z_t0=t0z_left(count_left)
+ count_left=count_left+1
+ end if
+ else
+ veloc_horiz = 0
+ veloc_vert = 0
+ traction_x_t0 = 0
+ traction_z_t0 = 0
+ endif
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- csl = vsext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
- rho_vp = rhol*cpl
- rho_vs = rhol*csl
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
- xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
- zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xgamma**2 + zgamma**2)
- nx = - zgamma / jacobian1D
- nz = + xgamma / jacobian1D
+ xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
+ zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = - zgamma / jacobian1D
+ nz = + xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
+ weight = jacobian1D * wzgll(j)
-! Clayton-Engquist condition if elastic
- if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob) - veloc_horiz
- vy = veloc_elastic(2,iglob)
- vz = veloc_elastic(3,iglob) - veloc_vert
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ ! vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
- vn = nx*vx+nz*vz
+ vn = nx*vx+nz*vz
- tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
- ty = rho_vs*vy
- tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ! ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
+ ! accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
- if(save_forward .and. isolver ==1) then
- if(p_sv)then !P-SV waves
- b_absorb_elastic_left(1,j,ib_xmin(ispecabs),it) = tx*weight
- b_absorb_elastic_left(3,j,ib_xmin(ispecabs),it) = tz*weight
- else !SH (membrane) waves
- b_absorb_elastic_left(2,j,ib_xmin(ispecabs),it) = ty*weight
- endif
- elseif(isolver == 2) then
- if(p_sv)then !P-SV waves
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_left(1,j,ib_xmin(ispecabs),NSTEP-it+1)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_left(3,j,ib_xmin(ispecabs),NSTEP-it+1)
- else !SH (membrane) waves
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_left(2,j,ib_xmin(ispecabs),NSTEP-it+1)
- endif
- endif
+ if(save_forward .and. isolver ==1) then
+ if(p_sv)then !P-SV waves
+ b_absorb_elastic_left(1,j,ib_xmin(ispecabs),it) = tx*weight
+ b_absorb_elastic_left(3,j,ib_xmin(ispecabs),it) = tz*weight
+ else !SH (membrane) waves
+ b_absorb_elastic_left(2,j,ib_xmin(ispecabs),it) = ty*weight
+ endif
+ elseif(isolver == 2) then
+ if(p_sv)then !P-SV waves
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
+ b_absorb_elastic_left(1,j,ib_xmin(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
+ b_absorb_elastic_left(3,j,ib_xmin(ispecabs),NSTEP-it+1)
+ else !SH (membrane) waves
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ b_absorb_elastic_left(2,j,ib_xmin(ispecabs),NSTEP-it+1)
+ endif
+ endif
- endif
+ endif
- enddo
+ enddo
- endif ! end of left absorbing boundary
+ endif ! end of left absorbing boundary
-!--- right absorbing boundary
- if(codeabs(IRIGHT,ispecabs)) then
+ !--- right absorbing boundary
+ if(codeabs(IRIGHT,ispecabs)) then
- i = NGLLX
+ i = NGLLX
- do j = 1,NGLLZ
+ do j = 1,NGLLZ
- iglob = ibool(i,j,ispec)
+ iglob = ibool(i,j,ispec)
-! for analytical initial plane wave for Bielak's conditions
-! left or right edge, horizontal normal vector
- if(add_Bielak_conditions .and. initialfield) then
- if (.not.over_critical_angle) then
- call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
- c_inc, c_refl, time_offset,f0)
- traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
- traction_z_t0 = mul_relaxed*(dxUz + dzUx)
- else
- veloc_horiz=v0x_right(count_right)
- veloc_vert=v0z_right(count_right)
- traction_x_t0=t0x_right(count_right)
- traction_z_t0=t0z_right(count_right)
- count_right=count_right+1
- end if
- else
- veloc_horiz = 0
- veloc_vert = 0
- traction_x_t0 = 0
- traction_z_t0 = 0
- endif
+ ! for analytical initial plane wave for Bielak's conditions
+ ! left or right edge, horizontal normal vector
+ if(add_Bielak_conditions .and. initialfield) then
+ if (.not.over_critical_angle) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ c_inc, c_refl, time_offset,f0)
+ traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
+ traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ else
+ veloc_horiz=v0x_right(count_right)
+ veloc_vert=v0z_right(count_right)
+ traction_x_t0=t0x_right(count_right)
+ traction_z_t0=t0z_right(count_right)
+ count_right=count_right+1
+ end if
+ else
+ veloc_horiz = 0
+ veloc_vert = 0
+ traction_x_t0 = 0
+ traction_z_t0 = 0
+ endif
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- csl = vsext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
- rho_vp = rhol*cpl
- rho_vs = rhol*csl
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
- xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
- zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xgamma**2 + zgamma**2)
- nx = + zgamma / jacobian1D
- nz = - xgamma / jacobian1D
+ xgamma = - xiz(i,j,ispec) * jacobian(i,j,ispec)
+ zgamma = + xix(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = + zgamma / jacobian1D
+ nz = - xgamma / jacobian1D
- weight = jacobian1D * wzgll(j)
+ weight = jacobian1D * wzgll(j)
-! Clayton-Engquist condition if elastic
- if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob) - veloc_horiz
- vy = veloc_elastic(2,iglob)
- vz = veloc_elastic(3,iglob) - veloc_vert
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
- vn = nx*vx+nz*vz
+ vn = nx*vx+nz*vz
- tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
- ty = rho_vs*vy
- tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
- if(save_forward .and. isolver ==1) then
- if(p_sv)then !P-SV waves
- b_absorb_elastic_right(1,j,ib_xmax(ispecabs),it) = tx*weight
- b_absorb_elastic_right(3,j,ib_xmax(ispecabs),it) = tz*weight
- else! SH (membrane) waves
- b_absorb_elastic_right(2,j,ib_xmax(ispecabs),it) = ty*weight
- endif
- elseif(isolver == 2) then
- if(p_sv)then !P-SV waves
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_right(1,j,ib_xmax(ispecabs),NSTEP-it+1)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_right(3,j,ib_xmax(ispecabs),NSTEP-it+1)
- else! SH (membrane) waves
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_right(2,j,ib_xmax(ispecabs),NSTEP-it+1)
- endif
- endif
+ if(save_forward .and. isolver ==1) then
+ if(p_sv)then !P-SV waves
+ b_absorb_elastic_right(1,j,ib_xmax(ispecabs),it) = tx*weight
+ b_absorb_elastic_right(3,j,ib_xmax(ispecabs),it) = tz*weight
+ else! SH (membrane) waves
+ b_absorb_elastic_right(2,j,ib_xmax(ispecabs),it) = ty*weight
+ endif
+ elseif(isolver == 2) then
+ if(p_sv)then !P-SV waves
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
+ b_absorb_elastic_right(1,j,ib_xmax(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
+ b_absorb_elastic_right(3,j,ib_xmax(ispecabs),NSTEP-it+1)
+ else! SH (membrane) waves
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ b_absorb_elastic_right(2,j,ib_xmax(ispecabs),NSTEP-it+1)
+ endif
+ endif
- endif
+ endif
- enddo
+ enddo
- endif ! end of right absorbing boundary
+ endif ! end of right absorbing boundary
-!--- bottom absorbing boundary
- if(codeabs(IBOTTOM,ispecabs)) then
+ !--- bottom absorbing boundary
+ if(codeabs(IBOTTOM,ispecabs)) then
- j = 1
+ j = 1
-! exclude corners to make sure there is no contradiction on the normal
- ibegin = 1
- iend = NGLLX
- if(codeabs(ILEFT,ispecabs)) ibegin = 2
- if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+ ! exclude corners to make sure there is no contradiction on the normal
+ ibegin = 1
+ iend = NGLLX
+ if(codeabs(ILEFT,ispecabs)) ibegin = 2
+ if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
- do i = ibegin,iend
+ do i = ibegin,iend
- iglob = ibool(i,j,ispec)
+ iglob = ibool(i,j,ispec)
-! for analytical initial plane wave for Bielak's conditions
-! top or bottom edge, vertical normal vector
- if(add_Bielak_conditions .and. initialfield) then
- if (.not.over_critical_angle) then
- call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
- c_inc, c_refl, time_offset,f0)
- traction_x_t0 = mul_relaxed*(dxUz + dzUx)
- traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
- else
- veloc_horiz=v0x_bot(count_bot)
- veloc_vert=v0z_bot(count_bot)
- traction_x_t0=t0x_bot(count_bot)
- traction_z_t0=t0z_bot(count_bot)
- count_bot=count_bot+1
- end if
- else
- veloc_horiz = 0
- veloc_vert = 0
- traction_x_t0 = 0
- traction_z_t0 = 0
- endif
+ ! for analytical initial plane wave for Bielak's conditions
+ ! top or bottom edge, vertical normal vector
+ if(add_Bielak_conditions .and. initialfield) then
+ if (.not.over_critical_angle) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ c_inc, c_refl, time_offset,f0)
+ traction_x_t0 = mul_relaxed*(dxUz + dzUx)
+ traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
+ else
+ veloc_horiz=v0x_bot(count_bot)
+ veloc_vert=v0z_bot(count_bot)
+ traction_x_t0=t0x_bot(count_bot)
+ traction_z_t0=t0z_bot(count_bot)
+ count_bot=count_bot+1
+ end if
+ else
+ veloc_horiz = 0
+ veloc_vert = 0
+ traction_x_t0 = 0
+ traction_z_t0 = 0
+ endif
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- csl = vsext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
- rho_vp = rhol*cpl
- rho_vs = rhol*csl
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
- xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
- zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xxi**2 + zxi**2)
- nx = + zxi / jacobian1D
- nz = - xxi / jacobian1D
+ xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
+ zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = + zxi / jacobian1D
+ nz = - xxi / jacobian1D
- weight = jacobian1D * wxgll(i)
+ weight = jacobian1D * wxgll(i)
-! Clayton-Engquist condition if elastic
- if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob) - veloc_horiz
- vy = veloc_elastic(2,iglob)
- vz = veloc_elastic(3,iglob) - veloc_vert
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
- vn = nx*vx+nz*vz
+ vn = nx*vx+nz*vz
- tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
- ty = rho_vs*vy
- tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz + traction_z_t0)*weight
- if(save_forward .and. isolver ==1) then
- if(p_sv)then !P-SV waves
- b_absorb_elastic_bottom(1,i,ib_zmin(ispecabs),it) = tx*weight
- b_absorb_elastic_bottom(3,i,ib_zmin(ispecabs),it) = tz*weight
- else!SH (membrane) waves
- b_absorb_elastic_bottom(2,i,ib_zmin(ispecabs),it) = ty*weight
- endif
- elseif(isolver == 2) then
- if(p_sv)then !P-SV waves
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_bottom(1,i,ib_zmin(ispecabs),NSTEP-it+1)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_bottom(3,i,ib_zmin(ispecabs),NSTEP-it+1)
- else!SH (membrane) waves
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_bottom(2,i,ib_zmin(ispecabs),NSTEP-it+1)
- endif
- endif
+ if(save_forward .and. isolver ==1) then
+ if(p_sv)then !P-SV waves
+ b_absorb_elastic_bottom(1,i,ib_zmin(ispecabs),it) = tx*weight
+ b_absorb_elastic_bottom(3,i,ib_zmin(ispecabs),it) = tz*weight
+ else!SH (membrane) waves
+ b_absorb_elastic_bottom(2,i,ib_zmin(ispecabs),it) = ty*weight
+ endif
+ elseif(isolver == 2) then
+ if(p_sv)then !P-SV waves
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
+ b_absorb_elastic_bottom(1,i,ib_zmin(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
+ b_absorb_elastic_bottom(3,i,ib_zmin(ispecabs),NSTEP-it+1)
+ else!SH (membrane) waves
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ b_absorb_elastic_bottom(2,i,ib_zmin(ispecabs),NSTEP-it+1)
+ endif
+ endif
- endif
+ endif
- enddo
+ enddo
- endif ! end of bottom absorbing boundary
+ endif ! end of bottom absorbing boundary
-!--- top absorbing boundary
- if(codeabs(ITOP,ispecabs)) then
+ !--- top absorbing boundary
+ if(codeabs(ITOP,ispecabs)) then
- j = NGLLZ
+ j = NGLLZ
-! exclude corners to make sure there is no contradiction on the normal
- ibegin = 1
- iend = NGLLX
- if(codeabs(ILEFT,ispecabs)) ibegin = 2
- if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+ ! exclude corners to make sure there is no contradiction on the normal
+ ibegin = 1
+ iend = NGLLX
+ if(codeabs(ILEFT,ispecabs)) ibegin = 2
+ if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
- do i = ibegin,iend
+ do i = ibegin,iend
- iglob = ibool(i,j,ispec)
+ iglob = ibool(i,j,ispec)
-! for analytical initial plane wave for Bielak's conditions
-! top or bottom edge, vertical normal vector
- if(add_Bielak_conditions .and. initialfield) then
- call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
- c_inc, c_refl, time_offset,f0)
- traction_x_t0 = mul_relaxed*(dxUz + dzUx)
- traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
- else
- veloc_horiz = 0
- veloc_vert = 0
- traction_x_t0 = 0
- traction_z_t0 = 0
- endif
+ ! for analytical initial plane wave for Bielak's conditions
+ ! top or bottom edge, vertical normal vector
+ if(add_Bielak_conditions .and. initialfield) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ c_inc, c_refl, time_offset,f0)
+ traction_x_t0 = mul_relaxed*(dxUz + dzUx)
+ traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
+ else
+ veloc_horiz = 0
+ veloc_vert = 0
+ traction_x_t0 = 0
+ traction_z_t0 = 0
+ endif
-! external velocity model
- if(assign_external_model) then
- cpl = vpext(i,j,ispec)
- csl = vsext(i,j,ispec)
- rhol = rhoext(i,j,ispec)
- endif
+ ! external velocity model
+ if(assign_external_model) then
+ cpl = vpext(i,j,ispec)
+ csl = vsext(i,j,ispec)
+ rhol = rhoext(i,j,ispec)
+ endif
- rho_vp = rhol*cpl
- rho_vs = rhol*csl
+ rho_vp = rhol*cpl
+ rho_vs = rhol*csl
- xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
- zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
- jacobian1D = sqrt(xxi**2 + zxi**2)
- nx = - zxi / jacobian1D
- nz = + xxi / jacobian1D
+ xxi = + gammaz(i,j,ispec) * jacobian(i,j,ispec)
+ zxi = - gammax(i,j,ispec) * jacobian(i,j,ispec)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = - zxi / jacobian1D
+ nz = + xxi / jacobian1D
- weight = jacobian1D * wxgll(i)
+ weight = jacobian1D * wxgll(i)
-! Clayton-Engquist condition if elastic
- if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob) - veloc_horiz
- vy = veloc_elastic(2,iglob)
- vz = veloc_elastic(3,iglob) - veloc_vert
+ ! Clayton-Engquist condition if elastic
+ if(elastic(ispec)) then
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ vy = veloc_elastic(2,iglob)
+ vz = veloc_elastic(3,iglob) - veloc_vert
- vn = nx*vx+nz*vz
+ vn = nx*vx+nz*vz
- tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
- ty = rho_vs*vy
- tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
+ tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
+ ty = rho_vs*vy
+ tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
- accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - ty*weight
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) - (tz - traction_z_t0)*weight
- if(save_forward .and. isolver ==1) then
- if(p_sv)then !P-SV waves
- b_absorb_elastic_top(1,i,ib_zmax(ispecabs),it) = tx*weight
- b_absorb_elastic_top(3,i,ib_zmax(ispecabs),it) = tz*weight
- else!SH (membrane) waves
- b_absorb_elastic_top(2,i,ib_zmax(ispecabs),it) = ty*weight
- endif
- elseif(isolver == 2) then
- if(p_sv)then !P-SV waves
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_top(1,i,ib_zmax(ispecabs),NSTEP-it+1)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_top(3,i,ib_zmax(ispecabs),NSTEP-it+1)
- else!SH (membrane) waves
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_top(2,i,ib_zmax(ispecabs),NSTEP-it+1)
- endif
- endif
+ if(save_forward .and. isolver ==1) then
+ if(p_sv)then !P-SV waves
+ b_absorb_elastic_top(1,i,ib_zmax(ispecabs),it) = tx*weight
+ b_absorb_elastic_top(3,i,ib_zmax(ispecabs),it) = tz*weight
+ else!SH (membrane) waves
+ b_absorb_elastic_top(2,i,ib_zmax(ispecabs),it) = ty*weight
+ endif
+ elseif(isolver == 2) then
+ if(p_sv)then !P-SV waves
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_top(1,i,ib_zmax(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_top(3,i,ib_zmax(ispecabs),NSTEP-it+1)
+ else!SH (membrane) waves
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_top(2,i,ib_zmax(ispecabs),NSTEP-it+1)
+ endif
+ endif
- endif
+ endif
- enddo
+ enddo
- endif ! end of top absorbing boundary
+ endif ! end of top absorbing boundary
- enddo
+ enddo
endif ! end of absorbing boundaries
-! --- add the source if it is a moment tensor
+ ! --- add the source if it is a moment tensor
if(.not. initialfield) then
- do i_source=1,NSOURCE
-! if this processor carries the source and the source element is elastic
- if (is_proc_source(i_source) == 1 .and. elastic(ispec_selected_source(i_source))) then
+ do i_source=1,NSOURCE
+ ! if this processor carries the source and the source element is elastic
+ if (is_proc_source(i_source) == 1 .and. elastic(ispec_selected_source(i_source))) then
-! moment tensor
- if(source_type(i_source) == 2) then
+ ! moment tensor
+ if(source_type(i_source) == 2) then
- if(.not.p_sv) call exit_MPI('cannot have moment tensor source in SH (membrane) waves calculation')
+ if(.not.p_sv) call exit_MPI('cannot have moment tensor source in SH (membrane) waves calculation')
- if(isolver == 1) then ! forward wavefield
-! add source array
- do j=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,ispec_selected_source(i_source))
- accel_elastic(1,iglob) = accel_elastic(1,iglob) + &
- sourcearray(i_source,1,i,j)*source_time_function(i_source,it)
- accel_elastic(3,iglob) = accel_elastic(3,iglob) + &
- sourcearray(i_source,2,i,j)*source_time_function(i_source,it)
- enddo
- enddo
- else ! backward wavefield
- do j=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,ispec_selected_source(i_source))
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + &
- sourcearray(i_source,1,i,j)*source_time_function(i_source,NSTEP-it+1)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + &
- sourcearray(i_source,2,i,j)*source_time_function(i_source,NSTEP-it+1)
- enddo
- enddo
- endif !endif isolver == 1
+ if(isolver == 1) then ! forward wavefield
+ ! add source array
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + &
+ sourcearray(i_source,1,i,j)*source_time_function(i_source,it)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + &
+ sourcearray(i_source,2,i,j)*source_time_function(i_source,it)
+ enddo
+ enddo
+ else ! backward wavefield
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec_selected_source(i_source))
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) + &
+ sourcearray(i_source,1,i,j)*source_time_function(i_source,NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + &
+ sourcearray(i_source,2,i,j)*source_time_function(i_source,NSTEP-it+1)
+ enddo
+ enddo
+ endif !endif isolver == 1
- endif !if(source_type(i_source) == 2)
+ endif !if(source_type(i_source) == 2)
- endif ! if this processor carries the source and the source element is elastic
- enddo ! do i_source=1,NSOURCE
+ endif ! if this processor carries the source and the source element is elastic
+ enddo ! do i_source=1,NSOURCE
- if(isolver == 2) then ! adjoint wavefield
+ if(isolver == 2) then ! adjoint wavefield
- irec_local = 0
- do irec = 1,nrec
-! add the source (only if this proc carries the source)
- if(myrank == which_proc_receiver(irec)) then
+ irec_local = 0
+ do irec = 1,nrec
+ ! add the source (only if this proc carries the source)
+ if(myrank == which_proc_receiver(irec)) then
- irec_local = irec_local + 1
- if(elastic(ispec_selected_rec(irec))) then
-! add source array
- do j=1,NGLLZ
- do i=1,NGLLX
- iglob = ibool(i,j,ispec_selected_rec(irec))
- if(p_sv)then !P-SH waves
- accel_elastic(1,iglob) = accel_elastic(1,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
- accel_elastic(3,iglob) = accel_elastic(3,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,3,i,j)
- else !SH (membrane) waves
- accel_elastic(2,iglob) = accel_elastic(2,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,2,i,j)
- endif
- enddo
- enddo
- endif ! if element is elastic
+ irec_local = irec_local + 1
+ if(elastic(ispec_selected_rec(irec))) then
+ ! add source array
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec_selected_rec(irec))
+ if(p_sv)then !P-SH waves
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,3,i,j)
+ else !SH (membrane) waves
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,2,i,j)
+ endif
+ enddo
+ enddo
+ endif ! if element is elastic
- endif ! if this processor carries the adjoint source and the source element is elastic
- enddo ! irec = 1,nrec
+ endif ! if this processor carries the adjoint source and the source element is elastic
+ enddo ! irec = 1,nrec
- endif ! if isolver == 2 adjoint wavefield
+ endif ! if isolver == 2 adjoint wavefield
endif ! if not using an initial field
-! implement attenuation
+ ! implement attenuation
if(TURN_ATTENUATION_ON) then
-! compute Grad(displ_elastic) at time step n+1 for attenuation
- call compute_gradient_attenuation(displ_elastic,dux_dxl_np1,duz_dxl_np1, &
- dux_dzl_np1,duz_dzl_np1,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,npoin)
+ ! compute Grad(displ_elastic) at time step n+1 for attenuation
+ call compute_gradient_attenuation(displ_elastic,dux_dxl_np1,duz_dxl_np1, &
+ dux_dzl_np1,duz_dzl_np1,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,npoin)
-! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
-! loop over spectral elements
- do ispec = 1,nspec
+ ! update memory variables with fourth-order Runge-Kutta time scheme for attenuation
+ ! loop over spectral elements
+ do ispec = 1,nspec
- do j=1,NGLLZ
- do i=1,NGLLX
+ do j=1,NGLLZ
+ do i=1,NGLLX
- theta_n = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
- theta_np1 = dux_dxl_np1(i,j,ispec) + duz_dzl_np1(i,j,ispec)
+ theta_n = dux_dxl_n(i,j,ispec) + duz_dzl_n(i,j,ispec)
+ theta_np1 = dux_dxl_np1(i,j,ispec) + duz_dzl_np1(i,j,ispec)
-! loop on all the standard linear solids
- do i_sls = 1,N_SLS
+ ! loop on all the standard linear solids
+ do i_sls = 1,N_SLS
-! evolution e1
- Un = e1(i,j,ispec,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,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) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e1(i,j,ispec,i_sls) = Unp1
+ ! evolution e1
+ Un = e1(i,j,ispec,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,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) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e1(i,j,ispec,i_sls) = Unp1
-! evolution e11
- Un = e11(i,j,ispec,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,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) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e11(i,j,ispec,i_sls) = Unp1
+ ! evolution e11
+ Un = e11(i,j,ispec,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,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) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e11(i,j,ispec,i_sls) = Unp1
-! evolution e13
- Un = e13(i,j,ispec,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,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) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e13(i,j,ispec,i_sls) = Unp1
+ ! evolution e13
+ Un = e13(i,j,ispec,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,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) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e13(i,j,ispec,i_sls) = Unp1
- enddo
+ enddo
- enddo
- enddo
- enddo
+ enddo
+ enddo
+ enddo
endif ! end of test on attenuation
- end subroutine compute_forces_elastic
+end subroutine compute_forces_elastic
Modified: seismo/2D/SPECFEM2D/trunk/compute_pressure.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/compute_pressure.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -45,8 +45,9 @@
subroutine compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext,e1,e11, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
+ numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,e1,e11, &
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute pressure in acoustic elements and in elastic elements
@@ -56,17 +57,20 @@
integer :: nspec,npoin,numat
+
integer, dimension(nspec) :: kmato
integer, dimension(NGLLX,NGLLX,nspec) :: ibool
double precision, dimension(2,numat) :: density
double precision, dimension(numat) :: porosity,tortuosity
double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(6,numat) :: anisotropy
double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,vsext,rhoext
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
- logical, dimension(nspec) :: elastic,poroelastic
+ logical, dimension(nspec) :: elastic,poroelastic,anisotropic
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic
real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
@@ -76,7 +80,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
- logical :: assign_external_model,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
+ logical :: assign_external_model,TURN_ATTENUATION_ON
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
@@ -95,8 +99,9 @@
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,&
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
+ numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,ispec,e1,e11, &
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
! use vector_field_display as temporary storage, store pressure in its second component
do j = 1,NGLLZ
@@ -117,8 +122,9 @@
subroutine compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,&
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
+ numat,kmato,density,porosity,tortuosity,elastcoef,vpext,vsext,rhoext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,ispec,e1,e11, &
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
! compute pressure in acoustic elements and in elastic elements
@@ -134,14 +140,16 @@
double precision, dimension(2,numat) :: density
double precision, dimension(numat) :: porosity,tortuosity
double precision, dimension(4,3,numat) :: elastcoef
+ double precision, dimension(6,numat) :: anisotropy
double precision, dimension(NGLLX,NGLLX,nspec) :: vpext,vsext,rhoext
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c15ext,c13ext,c33ext,c35ext,c55ext
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
! pressure in this element
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: pressure_element
- logical, dimension(nspec) :: elastic,poroelastic
+ logical, dimension(nspec) :: elastic,poroelastic,anisotropic
real(kind=CUSTOM_REAL), dimension(npoin) :: potential_dot_dot_acoustic
real(kind=CUSTOM_REAL), dimension(3,npoin) :: displ_elastic
real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displs_poroelastic,displw_poroelastic
@@ -150,7 +158,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
- logical :: assign_external_model,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
+ logical :: assign_external_model,TURN_ATTENUATION_ON
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11
@@ -181,6 +189,9 @@
real(kind=CUSTOM_REAL) :: D_biot,H_biot,C_biot,M_biot,rhol_bar
real(kind=CUSTOM_REAL) :: mul_G,lambdal_G,lambdalplus2mul_G
+! for anisotropy
+ double precision :: c11,c15,c13,c33,c35,c55
+
! if elastic element
!
! from L. S. Bennethum, Compressibility Moduli for Porous Materials Incorporating Volume Fraction,
@@ -203,6 +214,8 @@
! sigma_yy = lambda (epsilon_xx + epsilon_yy) + 2 mu epsilon_yy
! pressure = - trace(sigma) / 2 = - (lambda + mu) trace(epsilon)
!
+
+
if(elastic(ispec)) then
! get relaxed elastic parameters of current spectral element
@@ -220,7 +233,7 @@
denst = rhoext(i,j,ispec)
mul_relaxed = denst*csl*csl
lambdal_relaxed = denst*cpl*cpl - TWO*mul_relaxed
- endif
+ endif
! derivative along x and along z
dux_dxi = ZERO
@@ -288,11 +301,26 @@
endif
! full anisotropy
- if(TURN_ANISOTROPY_ON) then
+ if(anisotropic(ispec)) then
+ if(assign_external_model) then
+ c11 = c11ext(i,j,ispec)
+ c15 = c15ext(i,j,ispec)
+ c13 = c13ext(i,j,ispec)
+ c33 = c33ext(i,j,ispec)
+ c35 = c35ext(i,j,ispec)
+ c55 = c55ext(i,j,ispec)
+ else
+ c11 = anisotropy(1,kmato(ispec))
+ c13 = anisotropy(2,kmato(ispec))
+ c15 = anisotropy(3,kmato(ispec))
+ c33 = anisotropy(4,kmato(ispec))
+ c35 = anisotropy(5,kmato(ispec))
+ c55 = anisotropy(6,kmato(ispec))
+ endif
! implement anisotropy in 2D
- sigma_xx = c11val*dux_dxl + c15val*(duz_dxl + dux_dzl) + c13val*duz_dzl
- sigma_zz = c13val*dux_dxl + c35val*(duz_dxl + dux_dzl) + c33val*duz_dzl
+ sigma_xx = c11*dux_dxl + c15*(duz_dxl + dux_dzl) + c13*duz_dzl
+ sigma_zz = c13*dux_dxl + c35*(duz_dxl + dux_dzl) + c33*duz_dzl
endif
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2010-03-12 20:23:22 UTC (rev 16405)
@@ -62,6 +62,9 @@
!--- end of Nicolas Le Goff's constants for an unstructured CUBIT/METIS/SCOTCH mesh
+! Read external model from file DATA/model_velocity.dat_input or use routine
+ logical, parameter :: READ_EXTERNAL_MODEL_FILE = .false.
+
! compute and output acoustic and elastic energy (slows down the code significantly)
logical, parameter :: OUTPUT_ENERGY = .false.
Modified: seismo/2D/SPECFEM2D/trunk/define_external_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/define_external_model.f90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/define_external_model.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -41,8 +41,10 @@
! The full text of the license is available in file "LICENSE".
!
!========================================================================
+
- subroutine define_external_model(x,y,iflag_element,myrank,rho,vp,vs)
+ subroutine define_external_model(x,y,iflag_element,myrank,rho,vp,vs,Qp_attenuation,&
+ Qs_attenuation,c11,c13,c15,c33,c35,c55 )
implicit none
@@ -56,17 +58,34 @@
double precision, intent(in) :: x,y
double precision, intent(out) :: rho,vp,vs
+ double precision, intent(out) :: Qp_attenuation,Qs_attenuation
+ double precision, intent(out) :: c11,c15,c13,c33,c35,c55
! dummy routine here, just to demonstrate how the model can be assigned
- if(myrank == 0 .and. iflag_element == 1 .or. x < 1700.d0 .or. y >= 2300.d0) then
- rho = 2000.d0
- vp = 3000.d0
- vs = vp / sqrt(3.d0)
- else
- rho = 2500.d0
- vp = 3600.d0
- vs = vp / 2.d0
- endif
+ if(myrank == 0 .and. iflag_element == 1 .or. x < 1700.d0 .or. y >= 2300.d0) then
+ rho = 2000.d0
+ vp = 3000.d0
+ vs = vp / sqrt(3.d0)
+ Qp_attenuation = 0
+ Qs_attenuation = 0
+ c11 = 169.d9
+ c13 = 122.d9
+ c15 = 0.d0
+ c33 = c11
+ c35 = 0.d0
+ c55 = 75.3d9
+ else
+ rho = 2500.d0
+ vp = 3600.d0
+ vs = vp / 2.d0
+ Qp_attenuation = 60
+ Qs_attenuation = 60
+ c11 = 0.d0
+ c13 = 0.d0
+ c15 = 0.d0
+ c33 = 0.d0
+ c35 = 0.d0
+ c55 = 0.d0
+ endif
end subroutine define_external_model
-
Modified: seismo/2D/SPECFEM2D/trunk/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/gmat01.f90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/gmat01.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -42,10 +42,10 @@
!
!========================================================================
- subroutine gmat01(density_array,porosity_array,tortuosity_array,permeability,poroelastcoef,&
- numat,myrank,ipass,Qp_array,Qs_array)
+subroutine gmat01(density_array,porosity_array,tortuosity_array,aniso_array,permeability,poroelastcoef,&
+ numat,myrank,ipass,Qp_array,Qs_array)
-! read properties of a 2D isotropic or anisotropic linear elastic element
+ ! read properties of a 2D isotropic or anisotropic linear elastic element
implicit none
@@ -56,7 +56,7 @@
integer numat,myrank,ipass
double precision density_array(2,numat),poroelastcoef(4,3,numat),porosity_array(numat)
- double precision tortuosity_array(numat),permeability(3,numat)
+ double precision aniso_array(6,numat),tortuosity_array(numat),permeability(3,numat)
double precision Qp_array(numat),Qs_array(numat)
integer in,n,indic
@@ -66,15 +66,16 @@
double precision cpIsquare,cpIIsquare,cssquare,mu_s,mu_fr,eta_f,lambda_s,lambda_fr
double precision val1,val2,val3,val4,val5,val6
double precision val7,val8,val9,val10,val11,val12,val0
- double precision c11,c13,c33,c44
+ double precision :: c11,c13,c15,c33,c35,c55
double precision afactor,bfactor,cfactor,D_biot,H_biot,C_biot,M_biot,density_bar
-!
-!---- loop over the different material sets
-!
+ !
+ !---- loop over the different material sets
+ !
density_array(:,:) = zero
porosity_array(:) = zero
tortuosity_array(:) = zero
+ aniso_array(:,:) = zero
permeability(:,:) = zero
poroelastcoef(:,:,:) = zero
Qp_array(:) = zero
@@ -87,275 +88,306 @@
read(IIN,"(a80)") datlin
do in = 1,numat
- read(IIN,*) n,indic,val0,val1,val2,val3,val4,val5,val6,val7,val8,val9,val10,val11,val12
+ read(IIN,*) n,indic,val0,val1,val2,val3,val4,val5,val6,val7,val8,val9,val10,val11,val12
- if(n<1 .or. n>numat) call exit_MPI('Wrong material set number')
+ if(n<1 .or. n>numat) call exit_MPI('Wrong material set number')
-!---- isotropic material, P and S velocities given, allows for declaration of elastic/acoustic material
-!---- elastic (cs/=0) and acoustic (cs=0)
- if(indic == 1) then
- density(1) = val0
+ !---- isotropic material, P and S velocities given, allows for declaration of elastic/acoustic material
+ !---- elastic (cs/=0) and acoustic (cs=0)
+ if(indic == 1) then
+ density(1) = val0
-! P and S velocity
- cp = val1
- cs = val2
+ ! P and S velocity
+ cp = val1
+ cs = val2
-! Qp and Qs values
- Qp = val5
- Qs = val6
+ ! Qp and Qs values
+ Qp = val5
+ Qs = val6
-! Lam'e parameters
- lambdaplus2mu = density(1)*cp*cp
- mu = density(1)*cs*cs
- two_mu = 2.d0*mu
- lambda = lambdaplus2mu - two_mu
+ ! Lam'e parameters
+ lambdaplus2mu = density(1)*cp*cp
+ mu = density(1)*cs*cs
+ two_mu = 2.d0*mu
+ lambda = lambdaplus2mu - two_mu
-! bulk modulus Kappa
- kappa = lambda + two_mu/3.d0
+ ! bulk modulus Kappa
+ kappa = lambda + two_mu/3.d0
-! Young modulus
- young = 9.d0*kappa*mu/(3.d0*kappa + mu)
+ ! Young modulus
+ young = 9.d0*kappa*mu/(3.d0*kappa + mu)
-! Poisson's ratio
- poisson = half*(3.d0*kappa-two_mu)/(3.d0*kappa+mu)
+ ! Poisson's ratio
+ poisson = half*(3.d0*kappa-two_mu)/(3.d0*kappa+mu)
-! Poisson's ratio must be between -1 and +1/2
- if (poisson < -1.d0 .or. poisson > 0.5d0) call exit_MPI('Poisson''s ratio out of range')
+ ! Poisson's ratio must be between -1 and +1/2
+ if (poisson < -1.d0 .or. poisson > 0.5d0) call exit_MPI('Poisson''s ratio out of range')
-!---- anisotropic material, c11, c13, c33 and c44 given in Pascal
- else if (indic == 2) then
+ !---- anisotropic material, c11, c13, c33 and c44 given in Pascal
+ else if (indic == 2) then
- density(1) =val0
- c11 = val1
- c13 = val2
- c33 = val3
- c44 = val4
+ density(1) =val0
-!---- isotropic material, moduli are given, allows for declaration of poroelastic material
-!---- poroelastic (<0phi<1)
- else if (indic == 3) then
-! Qs values
- Qs = val12
+ ! P and S velocity
+ cp = 20
+ cs = 10
- density(1) =val0
- density(2) =val1
+ ! Anisotropy parameters
+ c11 = val1
+ c13 = val2
+ c15 = val3
+ c33 = val4
+ c35 = val5
+ c55 = val6
-! Solid properties
- kappa_s = val7
- mu_s = val11
-! Fluid properties
- kappa_f = val8
- eta_f = val10
-! Frame properties
- kappa_fr = val9
- mu_fr = val11
-! Lam'e parameters for the solid phase and the frame
- lambdaplus2mu_s = kappa_s + FOUR_THIRDS*mu_s
- lambda_s = lambdaplus2mu_s - 2.d0*mu_s
- lambdaplus2mu_fr = kappa_fr + FOUR_THIRDS*mu_fr
- lambda_fr = lambdaplus2mu_fr - 2.d0*mu_fr
- phi = val2
- tortuosity = val3
+ ! Qp and Qs values
+ !Qp = val9
+ !Qs = val10
-! Biot coefficients for the input phi
- D_biot = kappa_s*(1.d0 + phi*(kappa_s/kappa_f - 1.d0))
- H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + FOUR_THIRDS*mu_fr
- C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
- M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
-! Approximated velocities (no viscous dissipation)
- density_bar = (1.d0 - phi)*density(1) + phi*density(2)
- afactor = density_bar - phi/tortuosity*density(2)
- bfactor = H_biot + phi*density_bar/(tortuosity*density(2))*M_biot - 2.d0*phi/tortuosity*C_biot
- cfactor = phi/(tortuosity*density(2))*(H_biot*M_biot - C_biot*C_biot)
- cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
- cssquare = val11/afactor
+ ! Lam'e parameters
+ lambdaplus2mu = density(1)*cp*cp
+ mu = density(1)*cs*cs
+ two_mu = 2.d0*mu
+ lambda = lambdaplus2mu - two_mu
- porosity_array(n) = val2
- tortuosity_array(n) = val3
- permeability(1,n) = val4
- permeability(2,n) = val5
- permeability(3,n) = val6
+ ! bulk modulus Kappa
+ kappa = lambda + two_mu/3.d0
-! Young modulus for the solid phase
- young_s = 9.d0*kappa_s*mu_s/(3.d0*kappa_s + mu_s)
+ ! Young modulus
+ young = 9.d0*kappa*mu/(3.d0*kappa + mu)
-! Poisson's ratio for the solid phase
- poisson_s = HALF*(3.d0*kappa_s- 2.d0*mu_s)/(3.d0*kappa_s+mu_s)
+ ! Poisson's ratio
+ poisson = half*(3.d0*kappa-two_mu)/(3.d0*kappa+mu)
-! Poisson's ratio must be between -1 and +1/2
- if (poisson_s < -1.d0 .or. poisson_s > 0.5d0) stop 'Poisson''s ratio for the solid phase out of range'
+ !---- isotropic material, moduli are given, allows for declaration of poroelastic material
+ !---- poroelastic (<0phi<1)
+ else if (indic == 3) then
+ ! Qs values
+ Qs = val12
- else
- call exit_MPI('wrong model flag read')
+ density(1) =val0
+ density(2) =val1
- endif
+ ! Solid properties
+ kappa_s = val7
+ mu_s = val11
+ ! Fluid properties
+ kappa_f = val8
+ eta_f = val10
+ ! Frame properties
+ kappa_fr = val9
+ mu_fr = val11
+ ! Lam'e parameters for the solid phase and the frame
+ lambdaplus2mu_s = kappa_s + FOUR_THIRDS*mu_s
+ lambda_s = lambdaplus2mu_s - 2.d0*mu_s
+ lambdaplus2mu_fr = kappa_fr + FOUR_THIRDS*mu_fr
+ lambda_fr = lambdaplus2mu_fr - 2.d0*mu_fr
+ phi = val2
+ tortuosity = val3
-!
-!---- set elastic coefficients and density
-!
-! Isotropic : lambda, mu, K (= lambda + 2*mu), zero
-! Transverse anisotropic : c11, c13, c33, c44
-!
- if(indic == 1) then
- density_array(1,n) = density(1)
- poroelastcoef(1,1,n) = lambda
- poroelastcoef(2,1,n) = mu
- poroelastcoef(3,1,n) = lambdaplus2mu
- poroelastcoef(4,1,n) = zero
- Qp_array(n) = Qp
- Qs_array(n) = Qs
- if(mu > TINYVAL) then
- porosity_array(n) = 0.d0
+ ! Biot coefficients for the input phi
+ D_biot = kappa_s*(1.d0 + phi*(kappa_s/kappa_f - 1.d0))
+ H_biot = (kappa_s - kappa_fr)*(kappa_s - kappa_fr)/(D_biot - kappa_fr) + kappa_fr + FOUR_THIRDS*mu_fr
+ C_biot = kappa_s*(kappa_s - kappa_fr)/(D_biot - kappa_fr)
+ M_biot = kappa_s*kappa_s/(D_biot - kappa_fr)
+ ! Approximated velocities (no viscous dissipation)
+ density_bar = (1.d0 - phi)*density(1) + phi*density(2)
+ afactor = density_bar - phi/tortuosity*density(2)
+ bfactor = H_biot + phi*density_bar/(tortuosity*density(2))*M_biot - 2.d0*phi/tortuosity*C_biot
+ cfactor = phi/(tortuosity*density(2))*(H_biot*M_biot - C_biot*C_biot)
+ cpIsquare = (bfactor + sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
+ cpIIsquare = (bfactor - sqrt(bfactor*bfactor - 4.d0*afactor*cfactor))/(2.d0*afactor)
+ cssquare = val11/afactor
+
+ porosity_array(n) = val2
+ tortuosity_array(n) = val3
+ permeability(1,n) = val4
+ permeability(2,n) = val5
+ permeability(3,n) = val6
+
+ ! Young modulus for the solid phase
+ young_s = 9.d0*kappa_s*mu_s/(3.d0*kappa_s + mu_s)
+
+ ! Poisson's ratio for the solid phase
+ poisson_s = HALF*(3.d0*kappa_s- 2.d0*mu_s)/(3.d0*kappa_s+mu_s)
+
+ ! Poisson's ratio must be between -1 and +1/2
+ if (poisson_s < -1.d0 .or. poisson_s > 0.5d0) stop 'Poisson''s ratio for the solid phase out of range'
+
else
- porosity_array(n) = 1.d0
+ call exit_MPI('wrong model flag read')
+
endif
- elseif(indic == 2) then
- density_array(1,n) = density(1)
- poroelastcoef(1,1,n) = c11
- poroelastcoef(2,1,n) = c13
- poroelastcoef(3,1,n) = c33
- poroelastcoef(4,1,n) = c44
- Qp_array(n) = Qp
- Qs_array(n) = Qs
- porosity_array(n) = 0.d0
- else
- density_array(1,n) = density(1)
- density_array(2,n) = density(2)
- poroelastcoef(1,1,n) = lambda_s
- poroelastcoef(2,1,n) = mu_s ! = mu_fr
- poroelastcoef(3,1,n) = lambdaplus2mu_s
- poroelastcoef(4,1,n) = zero
- poroelastcoef(1,2,n) = kappa_f
- poroelastcoef(2,2,n) = eta_f
- poroelastcoef(3,2,n) = zero
- poroelastcoef(4,2,n) = zero
+ !
+ !---- set elastic coefficients and density
+ !
+ ! Isotropic : lambda, mu, K (= lambda + 2*mu), zero
+ ! Transverse anisotropic : c11, c13, c33, c44
+ !
+ if(indic == 1) then
+ density_array(1,n) = density(1)
+ poroelastcoef(1,1,n) = lambda
+ poroelastcoef(2,1,n) = mu
+ poroelastcoef(3,1,n) = lambdaplus2mu
+ poroelastcoef(4,1,n) = zero
+ Qp_array(n) = Qp
+ Qs_array(n) = Qs
+ if(mu > TINYVAL) then
+ porosity_array(n) = 0.d0
+ else
+ porosity_array(n) = 1.d0
+ endif
+ elseif(indic == 2) then
+ density_array(1,n) = density(1)
+! dummy poroelastcoef values, trick to avoid floating invalid
+ poroelastcoef(1,1,n) = lambda
+ poroelastcoef(2,1,n) = mu
+ poroelastcoef(3,1,n) = lambdaplus2mu
+ poroelastcoef(4,1,n) = zero
+ aniso_array(1,n) = c11
+ aniso_array(2,n) = c13
+ aniso_array(3,n) = c15
+ aniso_array(4,n) = c33
+ aniso_array(5,n) = c35
+ aniso_array(6,n) = c55
+! dummy Q values, trick to avoid a bug in attenuation_model
+ Qp_array(n) = 15
+ Qs_array(n) = 15
+ porosity_array(n) = 0.d0
+ else
+ density_array(1,n) = density(1)
+ density_array(2,n) = density(2)
+ poroelastcoef(1,1,n) = lambda_s
+ poroelastcoef(2,1,n) = mu_s ! = mu_fr
+ poroelastcoef(3,1,n) = lambdaplus2mu_s
+ poroelastcoef(4,1,n) = zero
- poroelastcoef(1,3,n) = lambda_fr
- poroelastcoef(2,3,n) = mu_fr
- poroelastcoef(3,3,n) = lambdaplus2mu_fr
- poroelastcoef(4,3,n) = zero
- Qp_array(n) = 10.d0 ! dummy for attenuation_model
- Qs_array(n) = Qs
- endif
+ poroelastcoef(1,2,n) = kappa_f
+ poroelastcoef(2,2,n) = eta_f
+ poroelastcoef(3,2,n) = zero
+ poroelastcoef(4,2,n) = zero
-!
-!---- check what has been read
-!
- if(myrank == 0 .and. ipass == 1) then
- if(indic == 1) then
-! material can be acoustic (fluid) or elastic (solid)
- if(poroelastcoef(2,1,n) > TINYVAL) then ! elastic
- write(IOUT,200) n,cp,cs,density(1),poisson,lambda,mu,kappa,young,Qp,Qs
- else ! acoustic
- write(IOUT,300) n,cp,density(1),kappa,Qp,Qs
- endif
- elseif(indic == 2) then ! elastic (anisotropic)
- write(IOUT,400) n,c11,c13,c33,c44,density,sqrt(c33/density),sqrt(c11/density),sqrt(c44/density),sqrt(c44/density),Qp,Qs
- else
-! material is poroelastic (solid/fluid)
- write(iout,500) n,sqrt(cpIsquare),sqrt(cpIIsquare),sqrt(cssquare)
- write(iout,600) density(1),poisson_s,lambda_s,mu_s,kappa_s,young_s
- write(iout,700) density(2),kappa_f,eta_f
- write(iout,800) lambda_fr,mu_fr,kappa_fr,porosity_array(n),tortuosity_array(n),&
- permeability(1,n),permeability(2,n),permeability(3,n),Qs
- write(iout,900) D_biot,H_biot,C_biot,M_biot
- endif
- endif
+ poroelastcoef(1,3,n) = lambda_fr
+ poroelastcoef(2,3,n) = mu_fr
+ poroelastcoef(3,3,n) = lambdaplus2mu_fr
+ poroelastcoef(4,3,n) = zero
+ Qp_array(n) = 10.d0 ! dummy for attenuation_model
+ Qs_array(n) = Qs
+ endif
+ !
+ !---- check what has been read
+ !
+ if(myrank == 0 .and. ipass == 1) then
+ if(indic == 1) then
+ ! material can be acoustic (fluid) or elastic (solid)
+ if(poroelastcoef(2,1,n) > TINYVAL) then ! elastic
+ write(IOUT,200) n,cp,cs,density(1),poisson,lambda,mu,kappa,young,Qp,Qs
+ else ! acoustic
+ write(IOUT,300) n,cp,density(1),kappa,Qp,Qs
+ endif
+ elseif(indic == 2) then ! elastic (anisotropic)
+ write(IOUT,400) n,density(1),c11,c13,c15,c33,c35,c55
+ else
+ ! material is poroelastic (solid/fluid)
+ write(iout,500) n,sqrt(cpIsquare),sqrt(cpIIsquare),sqrt(cssquare)
+ write(iout,600) density(1),poisson_s,lambda_s,mu_s,kappa_s,young_s
+ write(iout,700) density(2),kappa_f,eta_f
+ write(iout,800) lambda_fr,mu_fr,kappa_fr,porosity_array(n),tortuosity_array(n),&
+ permeability(1,n),permeability(2,n),permeability(3,n),Qs
+ write(iout,900) D_biot,H_biot,C_biot,M_biot
+ endif
+ endif
+
enddo
-!
-!---- formats
-!
- 100 format(//,' M a t e r i a l s e t s : ', &
- ' 2 D (p o r o) e l a s t i c i t y', &
- /1x,54('='),//5x,'Number of material sets . . . . . . (numat) =',i6)
+ !
+ !---- formats
+ !
+100 format(//,' M a t e r i a l s e t s : ', &
+ ' 2 D (p o r o) e l a s t i c i t y', &
+ /1x,54('='),//5x,'Number of material sets . . . . . . (numat) =',i6)
- 200 format(//5x,'----------------------------------------',/5x, &
- '-- Elastic (solid) isotropic material --',/5x, &
- '----------------------------------------',/5x, &
- 'Material set number. . . . . . . . (jmat) =',i6,/5x, &
- 'P-wave velocity. . . . . . . . . . . (cp) =',1pe15.8,/5x, &
- 'S-wave velocity. . . . . . . . . . . (cs) =',1pe15.8,/5x, &
- 'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
- 'Poisson''s ratio. . . . . . . . .(poisson) =',1pe15.8,/5x, &
- '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,/5x, &
- 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
- 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
+200 format(//5x,'----------------------------------------',/5x, &
+ '-- Elastic (solid) isotropic material --',/5x, &
+ '----------------------------------------',/5x, &
+ 'Material set number. . . . . . . . (jmat) =',i6,/5x, &
+ 'P-wave velocity. . . . . . . . . . . (cp) =',1pe15.8,/5x, &
+ 'S-wave velocity. . . . . . . . . . . (cs) =',1pe15.8,/5x, &
+ 'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
+ 'Poisson''s ratio. . . . . . . . .(poisson) =',1pe15.8,/5x, &
+ '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,/5x, &
+ 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
+ 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
- 300 format(//5x,'-------------------------------',/5x, &
- '-- Acoustic (fluid) material --',/5x, &
- '-------------------------------',/5x, &
- '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,/5x, &
- 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
- 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
+300 format(//5x,'-------------------------------',/5x, &
+ '-- Acoustic (fluid) material --',/5x, &
+ '-------------------------------',/5x, &
+ '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,/5x, &
+ 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
+ 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
- 400 format(//5x,'-------------------------------------',/5x, &
- '-- Transverse anisotropic material --',/5x, &
- '-------------------------------------',/5x, &
- 'Material set number. . . . . . . . (jmat) =',i6,/5x, &
- 'c11 coefficient (Pascal). . . . . . (c11) =',1pe15.8,/5x, &
- 'c13 coefficient (Pascal). . . . . . (c13) =',1pe15.8,/5x, &
- 'c33 coefficient (Pascal). . . . . . (c33) =',1pe15.8,/5x, &
- 'c44 coefficient (Pascal). . . . . . (c44) =',1pe15.8,/5x, &
- 'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
- '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,/5x, &
- 'Qp_attenuation. . . . . . . . . . . .(Qp) =',1pe15.8,/5x, &
- 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
+400 format(//5x,'-------------------------------------',/5x, &
+ '-- Transverse anisotropic material --',/5x, &
+ '-------------------------------------',/5x, &
+ 'Material set number. . . . . . . . (jmat) =',i6,/5x, &
+ 'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
+ 'c11 coefficient (Pascal). . . . . . (c11) =',1pe15.8,/5x, &
+ 'c13 coefficient (Pascal). . . . . . (c13) =',1pe15.8,/5x, &
+ 'c15 coefficient (Pascal). . . . . . (c15) =',1pe15.8,/5x, &
+ 'c33 coefficient (Pascal). . . . . . (c33) =',1pe15.8,/5x, &
+ 'c35 coefficient (Pascal). . . . . . (c35) =',1pe15.8,/5x, &
+ 'c55 coefficient (Pascal). . . . . . (c55) =',1pe15.8,/5x)
- 500 format(//5x,'----------------------------------------',/5x, &
- '-- Poroelastic isotropic material --',/5x, &
- '----------------------------------------',/5x, &
- 'Material set number. . . . . . . . (jmat) =',i6,/5x, &
- 'First P-wave velocity. . . . . . . . . . . (cpI) =',1pe15.8,/5x, &
- 'Second P-wave velocity. . . . . . . . . . . (cpII) =',1pe15.8,/5x, &
- 'S-wave velocity. . . . . . . . . . . (cs) =',1pe15.8)
+500 format(//5x,'----------------------------------------',/5x, &
+ '-- Poroelastic isotropic material --',/5x, &
+ '----------------------------------------',/5x, &
+ 'Material set number. . . . . . . . (jmat) =',i6,/5x, &
+ 'First P-wave velocity. . . . . . . . . . . (cpI) =',1pe15.8,/5x, &
+ 'Second P-wave velocity. . . . . . . . . . . (cpII) =',1pe15.8,/5x, &
+ 'S-wave velocity. . . . . . . . . . . (cs) =',1pe15.8)
- 600 format(//5x,'-------------------------------',/5x, &
- '-- Solid phase properties --',/5x, &
- 'Mass density. . . . . . . . . . (density_s) =',1pe15.8,/5x, &
- 'Poisson''s ratio. . . . . . . . .(poisson_s) =',1pe15.8,/5x, &
- 'First Lame parameter Lambda. . . (lambda_s) =',1pe15.8,/5x, &
- 'Second Lame parameter Mu. . . . . . .(mu_s) =',1pe15.8,/5x, &
- 'Solid bulk modulus Kappa . . . . . . . .(kappa_s) =',1pe15.8,/5x, &
- 'Young''s modulus E. . . . . . . . .(young_s) =',1pe15.8)
+600 format(//5x,'-------------------------------',/5x, &
+ '-- Solid phase properties --',/5x, &
+ 'Mass density. . . . . . . . . . (density_s) =',1pe15.8,/5x, &
+ 'Poisson''s ratio. . . . . . . . .(poisson_s) =',1pe15.8,/5x, &
+ 'First Lame parameter Lambda. . . (lambda_s) =',1pe15.8,/5x, &
+ 'Second Lame parameter Mu. . . . . . .(mu_s) =',1pe15.8,/5x, &
+ 'Solid bulk modulus Kappa . . . . . . . .(kappa_s) =',1pe15.8,/5x, &
+ 'Young''s modulus E. . . . . . . . .(young_s) =',1pe15.8)
- 700 format(//5x,'-------------------------------',/5x, &
- '-- Fluid phase properties --',/5x, &
- 'Mass density. . . . . . . . . . (density_f) =',1pe15.8,/5x, &
- 'Fluid bulk modulus Kappa . . . . . . . .(kappa_f) =',1pe15.8,/5x, &
- 'Fluid viscosity Eta . . . . . . . .(eta_f) =',1pe15.8)
+700 format(//5x,'-------------------------------',/5x, &
+ '-- Fluid phase properties --',/5x, &
+ 'Mass density. . . . . . . . . . (density_f) =',1pe15.8,/5x, &
+ 'Fluid bulk modulus Kappa . . . . . . . .(kappa_f) =',1pe15.8,/5x, &
+ 'Fluid viscosity Eta . . . . . . . .(eta_f) =',1pe15.8)
- 800 format(//5x,'-------------------------------',/5x, &
- '-- Frame properties --',/5x, &
- 'First Lame parameter Lambda. . . (lambda_fr) =',1pe15.8,/5x, &
- 'Second Lame parameter Mu. . . . . . .(mu_fr) =',1pe15.8,/5x, &
- 'Frame bulk modulus Kappa . . . . . . . .(kappa_fr) =',1pe15.8,/5x, &
- 'Porosity. . . . . . . . . . . . . . . . .(phi) =',1pe15.8,/5x,&
- 'Tortuosity. . . . . . . . . . . . . . . . .(c) =',1pe15.8,/5x,&
- 'Permeability xx component. . . . . . . . . . =',1pe15.8,/5x,&
- 'Permeability zx component. . . . . . . . . . =',1pe15.8,/5x,&
- 'Permeability zz component. . . . . . . . . . =',1pe15.8,/5x,&
- 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
+800 format(//5x,'-------------------------------',/5x, &
+ '-- Frame properties --',/5x, &
+ 'First Lame parameter Lambda. . . (lambda_fr) =',1pe15.8,/5x, &
+ 'Second Lame parameter Mu. . . . . . .(mu_fr) =',1pe15.8,/5x, &
+ 'Frame bulk modulus Kappa . . . . . . . .(kappa_fr) =',1pe15.8,/5x, &
+ 'Porosity. . . . . . . . . . . . . . . . .(phi) =',1pe15.8,/5x,&
+ 'Tortuosity. . . . . . . . . . . . . . . . .(c) =',1pe15.8,/5x,&
+ 'Permeability xx component. . . . . . . . . . =',1pe15.8,/5x,&
+ 'Permeability zx component. . . . . . . . . . =',1pe15.8,/5x,&
+ 'Permeability zz component. . . . . . . . . . =',1pe15.8,/5x,&
+ 'Qs_attenuation. . . . . . . . . . . .(Qs) =',1pe15.8)
- 900 format(//5x,'-------------------------------',/5x, &
- '-- Biot coefficients --',/5x, &
- '-------------------------------',/5x, &
- 'D. . . . . . . . =',1pe15.8,/5x, &
- 'H. . . . . . . . =',1pe15.8,/5x, &
- 'C. . . . . . . . =',1pe15.8,/5x, &
- 'M. . . . . . . . =',1pe15.8)
+900 format(//5x,'-------------------------------',/5x, &
+ '-- Biot coefficients --',/5x, &
+ '-------------------------------',/5x, &
+ 'D. . . . . . . . =',1pe15.8,/5x, &
+ 'H. . . . . . . . =',1pe15.8,/5x, &
+ 'C. . . . . . . . =',1pe15.8,/5x, &
+ 'M. . . . . . . . =',1pe15.8)
- end subroutine gmat01
+end subroutine gmat01
Added: seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/invert_mass_matrix.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -0,0 +1,74 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.0
+! ------------------------------
+!
+! Copyright Universite de Pau et des Pays de l'Adour, CNRS and INRIA, France,
+! and Princeton University, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+! Pieyre Le Loher, pieyre DOT le-loher aT inria.fr
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+
+subroutine invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic,npoin,rmass_inverse_elastic,&
+ rmass_inverse_acoustic,rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic)
+
+ implicit none
+ include 'constants.h'
+
+ logical any_elastic,any_acoustic,any_poroelastic
+
+ integer npoin
+
+! inverse mass matrices
+ real(kind=CUSTOM_REAL), dimension(npoin) :: rmass_inverse_elastic,rmass_inverse_acoustic
+ real(kind=CUSTOM_REAL), dimension(npoin) :: rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic
+
+
+! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
+ if(any_elastic) where(rmass_inverse_elastic <= 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
+ if(any_poroelastic) where(rmass_s_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_s_inverse_poroelastic = 1._CUSTOM_REAL
+ if(any_poroelastic) where(rmass_w_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_w_inverse_poroelastic = 1._CUSTOM_REAL
+ if(any_acoustic) where(rmass_inverse_acoustic <= 0._CUSTOM_REAL) rmass_inverse_acoustic = 1._CUSTOM_REAL
+
+! compute the inverse of the mass matrix
+ if(any_elastic) rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
+ if(any_poroelastic) rmass_s_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_s_inverse_poroelastic(:)
+ if(any_poroelastic) rmass_w_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_w_inverse_poroelastic(:)
+ if(any_acoustic) rmass_inverse_acoustic(:) = 1._CUSTOM_REAL / rmass_inverse_acoustic(:)
+
+end subroutine invert_mass_matrix
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -214,37 +214,37 @@
include "constants.h"
-! coordinates of the grid points of the mesh
+ ! coordinates of the grid points of the mesh
double precision, dimension(:,:), allocatable :: x,z
-! to compute the coordinate transformation
+ ! to compute the coordinate transformation
integer :: ioffset
double precision :: gamma,absx,a00,a01,bot0,top0
-! to store density and velocity model
- double precision, dimension(:), allocatable :: rho_s,cp,cs,aniso3,aniso4,Qp,Qs
+ ! to store density and velocity model
+ double precision, dimension(:), allocatable :: rho_s,cp,cs,aniso3,aniso4,aniso5,aniso6,aniso7,aniso8,Qp,Qs
double precision, dimension(:), allocatable :: rho_f,phi,tortuosity,permxx,permxz,&
- permzz,kappa_s,kappa_f,kappa_fr,eta_f,mu_fr
+ permzz,kappa_s,kappa_f,kappa_fr,eta_f,mu_fr
integer, dimension(:), allocatable :: icodemat
integer, dimension(:), allocatable :: num_material
-! interface data
+ ! interface data
integer interface_current,ipoint_current,number_of_interfaces,npoints_interface_bottom,npoints_interface_top
integer ilayer,number_of_layers,max_npoints_interface
double precision xinterface_dummy,zinterface_dummy,xinterface_dummy_previous
integer, dimension(:), allocatable :: nz_layer
double precision, dimension(:), allocatable :: &
- xinterface_bottom,zinterface_bottom,coefs_interface_bottom, &
- xinterface_top,zinterface_top,coefs_interface_top
+ xinterface_bottom,zinterface_bottom,coefs_interface_bottom, &
+ xinterface_top,zinterface_top,coefs_interface_top
-! for the source and receivers
+ ! for the source and receivers
integer, dimension(:), allocatable :: source_type,time_function_type
integer nrec_total,irec_global_number
double precision, dimension(:),allocatable :: xs,zs,f0,t0,angleforce,Mxx,Mzz,Mxz,factor
integer NSOURCE, NSOURCES, i_source, icounter, ios
logical, dimension(:),allocatable :: source_surf
double precision xrec,zrec
-! file number for source file
+ ! file number for source file
integer, parameter :: IIN_SOURCE = 22
character(len=150) dummystring
@@ -267,14 +267,14 @@
double precision tang1,tangN,vpregion,vsregion,poisson_ratio
double precision cutsnaps,sizemax_arrows,anglerec,xmin,xmax,deltat
double precision val0read,val1read,val2read,val3read,val4read,val5read,val6read,val7read,&
- val8read,val9read,val10read,val11read,val12read
+ val8read,val9read,val10read,val11read,val12read
double precision, dimension(:), allocatable :: xdeb,zdeb,xfin,zfin
logical interpol,gnuplot,assign_external_model,outputgrid
logical abstop,absbottom,absleft,absright,any_abs
logical meshvect,initialfield,modelvect,boundvect,add_Bielak_conditions
- logical TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON
+ logical TURN_ATTENUATION_ON,TURN_VISCATTENUATION_ON
double precision :: Q0,freq0
@@ -285,27 +285,30 @@
integer, external :: num_4, num_9
double precision, external :: value_spline
-! flag to save the last frame for kernels calculation purpose and type of simulation
+ ! flag to save the last frame for kernels calculation purpose and type of simulation
logical :: save_forward
integer :: isolver
-! flag to indicate an anisotropic material
+ ! flag to indicate an isotropic elastic/acoustic material
+ integer, parameter :: ISOTROPIC_MATERIAL = 1
+
+ ! flag to indicate an anisotropic material
integer, parameter :: ANISOTROPIC_MATERIAL = 2
-! flag to indicate a poroelastic material
+ ! flag to indicate a poroelastic material
integer, parameter :: POROELASTIC_MATERIAL = 3
-! file number for interface file
+ ! file number for interface file
integer, parameter :: IIN_INTERFACES = 15
-! ignore variable name field (junk) at the beginning of each input line
+ ! ignore variable name field (junk) at the beginning of each input line
logical, parameter :: IGNORE_JUNK = .true.,DONT_IGNORE_JUNK = .false.
-! parameters for external mesh
+ ! parameters for external mesh
logical :: read_external_mesh
character(len=256) :: mesh_file, nodes_coords_file, materials_file, free_surface_file, absorbing_surface_file
-! variables used for storing info about the mesh and partitions
+ ! variables used for storing info about the mesh and partitions
integer, dimension(:), pointer :: elmnts
integer, dimension(:), pointer :: elmnts_bis
double precision, dimension(:,:), pointer :: nodes_coords
@@ -329,13 +332,13 @@
integer, dimension(:), allocatable :: my_interfaces
integer, dimension(:), allocatable :: my_nb_interfaces
-! for acoustic/elastic coupled elements
+ ! for acoustic/elastic coupled elements
integer :: nedges_coupled, nedges_coupled_loc
integer, dimension(:,:), pointer :: edges_coupled
-! for acoustic/poroelastic coupled elements
+ ! for acoustic/poroelastic coupled elements
integer :: nedges_acporo_coupled, nedges_acporo_coupled_loc
integer, dimension(:,:), pointer :: edges_acporo_coupled
-! for poroelastic/elastic coupled elements
+ ! for poroelastic/elastic coupled elements
integer :: nedges_elporo_coupled, nedges_elporo_coupled_loc
integer, dimension(:,:), pointer :: edges_elporo_coupled
@@ -354,7 +357,7 @@
integer, dimension(:), pointer :: ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right
-! variables used for partitioning
+ ! variables used for partitioning
integer :: nproc
integer :: partitioning_method
character(len=256) :: partitioning_strategy
@@ -362,11 +365,11 @@
integer, dimension(0:4) :: metis_options
character(len=256) :: prname
-! variables used for attenuation
+ ! variables used for attenuation
integer :: N_SLS
double precision :: f0_attenuation
-! variables used for tangential detection
+ ! variables used for tangential detection
logical :: force_normal_to_surface,rec_normal_to_surface
character(len=256) :: tangential_detection_curve_file
integer :: nnodes_tangential_curve
@@ -380,16 +383,16 @@
-! ***
-! *** read the parameter file
-! ***
+ ! ***
+ ! *** read the parameter file
+ ! ***
print *,'Reading the parameter file ... '
print *
open(unit=IIN,file='DATA/Par_file',status='old')
-! read file names and path for output
+ ! read file names and path for output
call read_value_string(IIN,IGNORE_JUNK,title)
call read_value_string(IIN,IGNORE_JUNK,interfacesfile)
@@ -397,7 +400,7 @@
write(*,*) title
print *
-! read info about external mesh
+ ! read info about external mesh
call read_value_logical(IIN,IGNORE_JUNK,read_external_mesh)
call read_value_string(IIN,IGNORE_JUNK,mesh_file)
call read_value_string(IIN,IGNORE_JUNK,nodes_coords_file)
@@ -406,7 +409,7 @@
call read_value_string(IIN,IGNORE_JUNK,absorbing_surface_file)
call read_value_string(IIN,IGNORE_JUNK,tangential_detection_curve_file)
-! read info about partitioning
+ ! read info about partitioning
call read_value_integer(IIN,IGNORE_JUNK,nproc)
if ( nproc <= 0 ) then
print *, 'Number of processes (nproc) must be greater than or equal to one.'
@@ -446,7 +449,7 @@
end select
-! read grid parameters
+ ! read grid parameters
call read_value_double_precision(IIN,IGNORE_JUNK,xmin)
call read_value_double_precision(IIN,IGNORE_JUNK,xmax)
call read_value_integer(IIN,IGNORE_JUNK,nx)
@@ -460,15 +463,14 @@
call read_value_logical(IIN,IGNORE_JUNK,initialfield)
call read_value_logical(IIN,IGNORE_JUNK,add_Bielak_conditions)
call read_value_logical(IIN,IGNORE_JUNK,assign_external_model)
- call read_value_logical(IIN,IGNORE_JUNK,TURN_ANISOTROPY_ON)
call read_value_logical(IIN,IGNORE_JUNK,TURN_ATTENUATION_ON)
-! read viscous attenuation parameters (poroelastic media)
+ ! read viscous attenuation parameters (poroelastic media)
call read_value_logical(IIN,IGNORE_JUNK,TURN_VISCATTENUATION_ON)
call read_value_double_precision(IIN,IGNORE_JUNK,Q0)
call read_value_double_precision(IIN,IGNORE_JUNK,freq0)
-! determine if body or surface (membrane) waves calculation
+ ! determine if body or surface (membrane) waves calculation
call read_value_logical(IIN,IGNORE_JUNK,p_sv)
if ( read_external_mesh ) then
@@ -551,7 +553,7 @@
num_elmnt = num_elmnt + 1
enddo
enddo
- else
+ else
num_elmnt = 0
do j = 1, nzread
do i = 1, nxread
@@ -571,7 +573,7 @@
endif
endif
-! read absorbing boundaries parameters
+ ! read absorbing boundaries parameters
call read_value_logical(IIN,IGNORE_JUNK,any_abs)
call read_value_logical(IIN,IGNORE_JUNK,absbottom)
call read_value_logical(IIN,IGNORE_JUNK,absright)
@@ -584,12 +586,12 @@
absleft = .false.
endif
-! read time step parameters
+ ! read time step parameters
call read_value_integer(IIN,IGNORE_JUNK,nt)
call read_value_double_precision(IIN,IGNORE_JUNK,deltat)
call read_value_integer(IIN,IGNORE_JUNK,isolver)
-! read source parameters
+ ! read source parameters
call read_value_integer(IIN,IGNORE_JUNK,NSOURCE)
allocate(source_surf(NSOURCE))
allocate(xs(NSOURCE))
@@ -608,8 +610,8 @@
if(ios /= 0) stop 'error opening SOURCE file'
icounter = 0
do while(ios == 0)
- read(IIN_SOURCE,"(a)",iostat=ios) dummystring
- if(ios == 0) icounter = icounter + 1
+ read(IIN_SOURCE,"(a)",iostat=ios) dummystring
+ if(ios == 0) icounter = icounter + 1
enddo
close(IIN_SOURCE)
if(mod(icounter,NLINES_PER_SOURCE) /= 0) &
@@ -617,7 +619,7 @@
NSOURCES = icounter / NLINES_PER_SOURCE
if(NSOURCES < 1) stop 'need at least one source in SOURCE file'
if(NSOURCES /= NSOURCE) &
- stop 'total number of sources read is different than declared in Par_file'
+ stop 'total number of sources read is different than declared in Par_file'
open(unit=IIN_SOURCE,file='DATA/SOURCE',status='old',action='read')
do i_source=1,NSOURCE
@@ -634,15 +636,15 @@
call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,Mxz(i_source))
call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,factor(i_source))
-! if Dirac source time function, use a very thin Gaussian instead
-! if Heaviside source time function, use a very thin error function instead
- if(time_function_type(i_source) == 4 .or. time_function_type(i_source) == 5) f0(i_source) = 1.d0 / (10.d0 * deltat)
+ ! if Dirac source time function, use a very thin Gaussian instead
+ ! if Heaviside source time function, use a very thin error function instead
+ if(time_function_type(i_source) == 4 .or. time_function_type(i_source) == 5) f0(i_source) = 1.d0 / (10.d0 * deltat)
- ! time delay of the source in seconds, use a 20 % security margin (use 2 / f0 if error function)
+ ! time delay of the source in seconds, use a 20 % security margin (use 2 / f0 if error function)
if(time_function_type(i_source)== 5) then
- t0(i_source) = 2.0d0 / f0(i_source)+t0(i_source)
+ t0(i_source) = 2.0d0 / f0(i_source)+t0(i_source)
else
- t0(i_source) = 1.20d0 / f0(i_source)+t0(i_source)
+ t0(i_source) = 1.20d0 / f0(i_source)+t0(i_source)
endif
print *
@@ -661,16 +663,16 @@
call read_value_logical(IIN,IGNORE_JUNK,force_normal_to_surface)
-! read constants for attenuation
+ ! read constants for attenuation
call read_value_integer(IIN,IGNORE_JUNK,N_SLS)
call read_value_double_precision(IIN,IGNORE_JUNK,f0_attenuation)
-! if source is not a Dirac or Heavyside then f0_attenuation is f0 of the first source
+ ! 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
f0_attenuation = f0(1)
endif
-! read receiver line parameters
+ ! read receiver line parameters
call read_value_integer(IIN,IGNORE_JUNK,seismotype)
call read_value_logical(IIN,IGNORE_JUNK,save_forward)
call read_value_logical(IIN,IGNORE_JUNK,generate_STATIONS)
@@ -680,7 +682,7 @@
if(nreceiverlines < 1) stop 'number of receiver lines must be greater than 1'
-! allocate receiver line arrays
+ ! allocate receiver line arrays
allocate(nrec(nreceiverlines))
allocate(xdeb(nreceiverlines))
allocate(zdeb(nreceiverlines))
@@ -688,20 +690,20 @@
allocate(zfin(nreceiverlines))
allocate(enreg_surf(nreceiverlines))
-! loop on all the receiver lines
+ ! loop on all the receiver lines
do ireceiverlines = 1,nreceiverlines
- call read_value_integer(IIN,IGNORE_JUNK,nrec(ireceiverlines))
- call read_value_double_precision(IIN,IGNORE_JUNK,xdeb(ireceiverlines))
- call read_value_double_precision(IIN,IGNORE_JUNK,zdeb(ireceiverlines))
- call read_value_double_precision(IIN,IGNORE_JUNK,xfin(ireceiverlines))
- 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!'
- endif
+ call read_value_integer(IIN,IGNORE_JUNK,nrec(ireceiverlines))
+ call read_value_double_precision(IIN,IGNORE_JUNK,xdeb(ireceiverlines))
+ call read_value_double_precision(IIN,IGNORE_JUNK,zdeb(ireceiverlines))
+ call read_value_double_precision(IIN,IGNORE_JUNK,xfin(ireceiverlines))
+ 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!'
+ endif
enddo
-! read display parameters
+ ! read display parameters
call read_value_integer(IIN,IGNORE_JUNK,NTSTEP_BETWEEN_OUTPUT_INFO)
call read_value_logical(IIN,IGNORE_JUNK,output_postscript_snapshot)
call read_value_logical(IIN,IGNORE_JUNK,output_color_image)
@@ -717,15 +719,15 @@
call read_value_logical(IIN,IGNORE_JUNK,gnuplot)
call read_value_logical(IIN,IGNORE_JUNK,outputgrid)
-! can use only one point to display lower-left corner only for interpolated snapshot
+ ! can use only one point to display lower-left corner only for interpolated snapshot
if(pointsdisp < 3) then
- pointsdisp = 3
- plot_lowerleft_corner_only = .true.
+ pointsdisp = 3
+ plot_lowerleft_corner_only = .true.
else
- plot_lowerleft_corner_only = .false.
+ plot_lowerleft_corner_only = .false.
endif
-! read the different material materials
+ ! read the different material materials
call read_value_integer(IIN,IGNORE_JUNK,nb_materials)
if(nb_materials <= 0) stop 'Negative number of materials not allowed!'
@@ -734,6 +736,10 @@
allocate(cs(nb_materials))
allocate(aniso3(nb_materials))
allocate(aniso4(nb_materials))
+ allocate(aniso5(nb_materials))
+ allocate(aniso6(nb_materials))
+ allocate(aniso7(nb_materials))
+ allocate(aniso8(nb_materials))
allocate(Qp(nb_materials))
allocate(Qs(nb_materials))
allocate(rho_s(nb_materials))
@@ -755,6 +761,10 @@
cs(:) = 0.d0
aniso3(:) = 0.d0
aniso4(:) = 0.d0
+ aniso5(:) = 0.d0
+ aniso6(:) = 0.d0
+ aniso7(:) = 0.d0
+ aniso8(:) = 0.d0
Qp(:) = 0.d0
Qs(:) = 0.d0
rho_s(:) = 0.d0
@@ -772,88 +782,104 @@
num_material(:) = 0
do imaterial=1,nb_materials
- call read_material_parameters(IIN,DONT_IGNORE_JUNK,i,icodematread,val0read,val1read,val2read,val3read, &
- val4read,val5read,val6read,val7read,val8read,val9read,val10read,val11read,val12read)
- if(i < 1 .or. i > nb_materials) stop 'Wrong material number!'
- icodemat(i) = icodematread
+ call read_material_parameters(IIN,DONT_IGNORE_JUNK,i,icodematread,val0read,val1read,val2read,val3read, &
+ val4read,val5read,val6read,val7read,val8read,val9read,val10read,val11read,val12read)
+ if(i < 1 .or. i > nb_materials) stop 'Wrong material number!'
+ icodemat(i) = icodematread
- if(icodemat(i) /= POROELASTIC_MATERIAL) then
- rho_s(i) = val0read
- cp(i) = val1read
- cs(i) = val2read
- Qp(i) = val5read
- Qs(i) = val6read
+ if(icodemat(i) == ISOTROPIC_MATERIAL) then
+ rho_s(i) = val0read
+ cp(i) = val1read
+ cs(i) = val2read
+ Qp(i) = val5read
+ Qs(i) = val6read
- if(rho_s(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'
+ if(rho_s(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) = val3read
- aniso4(i) = val4read
- if(cs(i) /= 0.d0) then
- phi(i) = 0.d0 ! elastic
- else
- phi(i) = 1.d0 ! acoustic
- endif
- else ! poroelastic
- rho_s(i) = val0read
- rho_f(i) = val1read
- phi(i) = val2read
- tortuosity(i) = val3read
- permxx(i) = val4read
- permxz(i) = val5read
- permzz(i) = val6read
- kappa_s(i) = val7read
- kappa_f(i) = val8read
- kappa_fr(i) = val9read
- eta_f(i) = val10read
- mu_fr(i) = val11read
- Qs(i) = val12read
+ aniso3(i) = val3read
+ aniso4(i) = val4read
+ if(cs(i) /= 0.d0) then
+ phi(i) = 0.d0 ! elastic
+ else
+ phi(i) = 1.d0 ! acoustic
+ endif
+ elseif (icodemat(i) == ANISOTROPIC_MATERIAL) then
+ rho_s(i) = val0read
+ cp(i) = val1read
+ cs(i) = val2read
+ aniso3(i) = val3read
+ aniso4(i) = val4read
+ aniso5(i) = val5read
+ aniso6(i) = val6read
+ aniso7(i) = val7read
+ aniso8(i) = val8read
+ Qp(i) = val9read
+ Qs(i) = val10read
+ else ! poroelastic
+ rho_s(i) = val0read
+ rho_f(i) = val1read
+ phi(i) = val2read
+ tortuosity(i) = val3read
+ permxx(i) = val4read
+ permxz(i) = val5read
+ permzz(i) = val6read
+ kappa_s(i) = val7read
+ kappa_f(i) = val8read
+ kappa_fr(i) = val9read
+ eta_f(i) = val10read
+ mu_fr(i) = val11read
+ Qs(i) = val12read
- if(rho_s(i) <= 0.d0 .or. rho_f(i) <= 0.d0) stop 'negative value of density'
- if(phi(i) <= 0.d0 .or. tortuosity(i) <= 0.d0) stop 'negative value of porosity or tortuosity'
- if(kappa_s(i) <= 0.d0 .or. kappa_f(i) <= 0.d0 .or. kappa_fr(i) <= 0.d0 .or. mu_fr(i) <= 0.d0) stop 'negative value of modulus'
- if(Qs(i) <= 0.d0) stop 'negative value of Qs'
- endif
+ if(rho_s(i) <= 0.d0 .or. rho_f(i) <= 0.d0) stop 'negative value of density'
+ if(phi(i) <= 0.d0 .or. tortuosity(i) <= 0.d0) stop 'negative value of porosity or tortuosity'
+ if(kappa_s(i) <= 0.d0 .or. kappa_f(i) <= 0.d0 .or. kappa_fr(i) <= 0.d0 .or. mu_fr(i) <= 0.d0) then
+ stop 'negative value of modulus'
+ end if
+ if(Qs(i) <= 0.d0) stop 'negative value of Qs'
+ endif
enddo
print *
print *, 'Nb of solid, fluid or porous materials = ',nb_materials
print *
do i=1,nb_materials
- if(icodemat(i) /= ANISOTROPIC_MATERIAL .and. icodemat(i) /= POROELASTIC_MATERIAL) then
- print *,'Material #',i,' isotropic'
- print *,'rho,cp,cs = ',rho_s(i),cp(i),cs(i),Qp(i),Qs(i)
- if(cs(i) < TINYVAL) then
- print *,'Material is fluid'
- else
- print *,'Material is solid'
- endif
- elseif(icodemat(i) == POROELASTIC_MATERIAL) then
- print *,'Material #',i,' isotropic'
- print *,'rho_s, kappa_s= ',rho_s(i),kappa_s(i)
- print *,'rho_f, kappa_f, eta_f= ',rho_f(i),kappa_f(i),eta_f(i)
- print *,'phi, tortuosity, permxx, permxz, permzz= ',phi(i),tortuosity(i),permxx(i),permxz(i),permzz(i)
- print *,'kappa_fr, mu_fr, Qs= ',kappa_fr(i),mu_fr(i),Qs(i)
- print *,'Material is porous'
- else
- print *,'Material #',i,' anisotropic'
- print *,'rho,c11,c13,c33,c44 = ',rho_s(i),cp(i),cs(i),aniso3(i),aniso4(i),Qp(i),Qs(i)
- endif
- print *
+ if(icodemat(i) /= ANISOTROPIC_MATERIAL .and. icodemat(i) /= POROELASTIC_MATERIAL) then
+ print *,'Material #',i,' isotropic'
+ print *,'rho,cp,cs = ',rho_s(i),cp(i),cs(i),Qp(i),Qs(i)
+ if(cs(i) < TINYVAL) then
+ print *,'Material is fluid'
+ else
+ print *,'Material is solid'
+ endif
+ elseif(icodemat(i) == POROELASTIC_MATERIAL) then
+ print *,'Material #',i,' isotropic'
+ print *,'rho_s, kappa_s= ',rho_s(i),kappa_s(i)
+ print *,'rho_f, kappa_f, eta_f= ',rho_f(i),kappa_f(i),eta_f(i)
+ print *,'phi, tortuosity, permxx, permxz, permzz= ',phi(i),tortuosity(i),permxx(i),permxz(i),permzz(i)
+ print *,'kappa_fr, mu_fr, Qs= ',kappa_fr(i),mu_fr(i),Qs(i)
+ print *,'Material is porous'
+ else
+ print *,'Material #',i,' anisotropic'
+ print *,'rho,cp,cs = ',rho_s(i),cp(i),cs(i)
+ print*,'c11,c13,c15,c33,c35,c55 = ',aniso3(i),aniso4(i),aniso5(i),aniso6(i),aniso7(i),aniso8(i)
+ print *,'Qp,Qs = ',Qp(i),Qs(i)
+ endif
+ print *
enddo
-! tangential detection
+ ! tangential detection
if (force_normal_to_surface .or. rec_normal_to_surface) then
- open(unit=IIN,file=tangential_detection_curve_file,status='old',action='read')
- read(IIN,*) nnodes_tangential_curve
- allocate(nodes_tangential_curve(2,nnodes_tangential_curve))
- do i = 1, nnodes_tangential_curve
- read(IIN,*) nodes_tangential_curve(1,i), nodes_tangential_curve(2,i)
- enddo
- close(IIN)
+ open(unit=IIN,file=tangential_detection_curve_file,status='old',action='read')
+ read(IIN,*) nnodes_tangential_curve
+ allocate(nodes_tangential_curve(2,nnodes_tangential_curve))
+ do i = 1, nnodes_tangential_curve
+ read(IIN,*) nodes_tangential_curve(1,i), nodes_tangential_curve(2,i)
+ enddo
+ close(IIN)
else
- nnodes_tangential_curve = 0
- allocate(nodes_tangential_curve(2,1))
+ nnodes_tangential_curve = 0
+ allocate(nodes_tangential_curve(2,1))
endif
if ( read_external_mesh ) then
@@ -901,13 +927,17 @@
print *,'Qs = ',Qs(imaterial_number)
elseif(icodemat(imaterial_number) == POROELASTIC_MATERIAL) then
print *,'Material # ',imaterial_number,' isotropic'
- print *,'Material is poroelastic'
+ print *,'Material is poroelastic'
else
print *,'Material # ',imaterial_number,' anisotropic'
- print *,'c11 = ',cp(imaterial_number)
- print *,'c13 = ',cs(imaterial_number)
- print *,'c33 = ',aniso3(imaterial_number)
- print *,'c44 = ',aniso4(imaterial_number)
+ print *,'cp = ',cp(imaterial_number)
+ print *,'cs = ',cs(imaterial_number)
+ print *,'c11 = ',aniso3(imaterial_number)
+ print *,'c13 = ',aniso4(imaterial_number)
+ print *,'c15 = ',aniso5(imaterial_number)
+ print *,'c33 = ',aniso6(imaterial_number)
+ print *,'c35 = ',aniso7(imaterial_number)
+ print *,'c55 = ',aniso8(imaterial_number)
print *,'rho = ',rho_s(imaterial_number)
print *,'Qp = ',Qp(imaterial_number)
print *,'Qs = ',Qs(imaterial_number)
@@ -932,7 +962,7 @@
print *
print *,'Parameter file successfully read... '
-!---
+ !---
if(ngnod /= 4 .and. ngnod /= 9) stop 'ngnod different from 4 or 9!'
@@ -942,7 +972,7 @@
print *,'Control elements have ',ngnod,' nodes'
print *
-!---
+ !---
if ( .not. read_external_mesh ) then
! allocate arrays for the grid
@@ -1002,10 +1032,10 @@
! check if we are in the last layer, which contains topography,
! and modify the position of the source accordingly if it is located exactly at the surface
- do i_source=1,NSOURCE
- if(source_surf(i_source) .and. ilayer == number_of_layers) &
- zs(i_source) = value_spline(xs(i_source),xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
- enddo
+ do i_source=1,NSOURCE
+ if(source_surf(i_source) .and. ilayer == number_of_layers) &
+ zs(i_source) = value_spline(xs(i_source),xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
+ enddo
! compute the offset of this layer in terms of number of spectral elements below along Z
if(ilayer > 1) then
@@ -1090,42 +1120,42 @@
! count the number of acoustic free-surface elements
nelem_acoustic_surface = 0
-! if the surface is absorbing, it cannot be free at the same time
+ ! if the surface is absorbing, it cannot be free at the same time
if(.not. abstop) then
- j = nzread
- do i = 1,nxread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- endif
- enddo
+ j = nzread
+ do i = 1,nxread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ endif
+ enddo
endif
if(.not. absbottom) then
- j = 1
- do i = 1,nxread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- endif
- enddo
+ j = 1
+ do i = 1,nxread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ endif
+ enddo
endif
if(.not. absleft) then
- i = 1
- do j = 1,nzread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- endif
- enddo
+ i = 1
+ do j = 1,nzread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ endif
+ enddo
endif
if(.not. absright) then
- i = nxread
- do j = 1,nzread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- endif
- enddo
+ i = nxread
+ do j = 1,nzread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >= 1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ endif
+ enddo
endif
@@ -1134,56 +1164,56 @@
nelem_acoustic_surface = 0
if(.not. abstop) then
- j = nzread
- do i = 1,nxread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
- acoustic_surface(2,nelem_acoustic_surface) = 2
- acoustic_surface(3,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
- acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
- endif
- enddo
+ j = nzread
+ do i = 1,nxread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+ acoustic_surface(2,nelem_acoustic_surface) = 2
+ acoustic_surface(3,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
+ acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
+ endif
+ enddo
endif
if(.not. absbottom) then
- j = 1
- do i = 1,nxread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
- acoustic_surface(2,nelem_acoustic_surface) = 2
- acoustic_surface(3,nelem_acoustic_surface) = elmnts(0+ngnod*((j-1)*nxread+i-1))
- acoustic_surface(4,nelem_acoustic_surface) = elmnts(1+ngnod*((j-1)*nxread+i-1))
- endif
- enddo
+ j = 1
+ do i = 1,nxread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+ acoustic_surface(2,nelem_acoustic_surface) = 2
+ acoustic_surface(3,nelem_acoustic_surface) = elmnts(0+ngnod*((j-1)*nxread+i-1))
+ acoustic_surface(4,nelem_acoustic_surface) = elmnts(1+ngnod*((j-1)*nxread+i-1))
+ endif
+ enddo
endif
if(.not. absleft) then
- i = 1
- do j = 1,nzread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
- acoustic_surface(2,nelem_acoustic_surface) = 2
- acoustic_surface(3,nelem_acoustic_surface) = elmnts(0+ngnod*((j-1)*nxread+i-1))
- acoustic_surface(4,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
- endif
- enddo
+ i = 1
+ do j = 1,nzread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+ acoustic_surface(2,nelem_acoustic_surface) = 2
+ acoustic_surface(3,nelem_acoustic_surface) = elmnts(0+ngnod*((j-1)*nxread+i-1))
+ acoustic_surface(4,nelem_acoustic_surface) = elmnts(3+ngnod*((j-1)*nxread+i-1))
+ endif
+ enddo
endif
if(.not. absright) then
- i = nxread
- do j = 1,nzread
- imaterial_number = num_material((j-1)*nxread+i)
- if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
- nelem_acoustic_surface = nelem_acoustic_surface + 1
- acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
- acoustic_surface(2,nelem_acoustic_surface) = 2
- acoustic_surface(3,nelem_acoustic_surface) = elmnts(1+ngnod*((j-1)*nxread+i-1))
- acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
- endif
- enddo
+ i = nxread
+ do j = 1,nzread
+ imaterial_number = num_material((j-1)*nxread+i)
+ if(icodemat(imaterial_number) /= ANISOTROPIC_MATERIAL .and. phi(imaterial_number) >=1.d0 ) then
+ nelem_acoustic_surface = nelem_acoustic_surface + 1
+ acoustic_surface(1,nelem_acoustic_surface) = (j-1)*nxread + (i-1)
+ acoustic_surface(2,nelem_acoustic_surface) = 2
+ acoustic_surface(3,nelem_acoustic_surface) = elmnts(1+ngnod*((j-1)*nxread+i-1))
+ acoustic_surface(4,nelem_acoustic_surface) = elmnts(2+ngnod*((j-1)*nxread+i-1))
+ endif
+ enddo
endif
!
@@ -1238,84 +1268,84 @@
endif
-! compute min and max of X and Z in the grid
+ ! compute min and max of X and Z in the grid
print *
print *,'Min and max value of X in the grid = ',minval(nodes_coords(1,:)),maxval(nodes_coords(1,:))
print *,'Min and max value of Z in the grid = ',minval(nodes_coords(2,:)),maxval(nodes_coords(2,:))
print *
-! ***
-! *** create a Gnuplot file that displays the grid
-! ***
+ ! ***
+ ! *** create a Gnuplot file that displays the grid
+ ! ***
if ( .not. read_external_mesh ) then
- print *
- print *,'Saving the grid in Gnuplot format...'
+ print *
+ print *,'Saving the grid in Gnuplot format...'
- open(unit=20,file='OUTPUT_FILES/gridfile.gnu',status='unknown')
+ open(unit=20,file='OUTPUT_FILES/gridfile.gnu',status='unknown')
-! draw horizontal lines of the grid
- print *,'drawing horizontal lines of the grid'
- istepx = 1
- if(ngnod == 4) then
- istepz = 1
- else
- istepz = 2
- endif
- do ili=0,nz,istepz
- do icol=0,nx-istepx,istepx
- write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
- write(20,*) sngl(x(icol+istepx,ili)),sngl(z(icol+istepx,ili))
- write(20,10)
- enddo
- enddo
+ ! draw horizontal lines of the grid
+ print *,'drawing horizontal lines of the grid'
+ istepx = 1
+ if(ngnod == 4) then
+ istepz = 1
+ else
+ istepz = 2
+ endif
+ do ili=0,nz,istepz
+ do icol=0,nx-istepx,istepx
+ write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
+ write(20,*) sngl(x(icol+istepx,ili)),sngl(z(icol+istepx,ili))
+ write(20,10)
+ enddo
+ enddo
-! draw vertical lines of the grid
- print *,'drawing vertical lines of the grid'
- if(ngnod == 4) then
- istepx = 1
- else
- istepx = 2
- endif
- istepz = 1
- do icol=0,nx,istepx
- do ili=0,nz-istepz,istepz
- write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
- write(20,*) sngl(x(icol,ili+istepz)),sngl(z(icol,ili+istepz))
- write(20,10)
- enddo
- enddo
+ ! draw vertical lines of the grid
+ print *,'drawing vertical lines of the grid'
+ if(ngnod == 4) then
+ istepx = 1
+ else
+ istepx = 2
+ endif
+ istepz = 1
+ do icol=0,nx,istepx
+ do ili=0,nz-istepz,istepz
+ write(20,*) sngl(x(icol,ili)),sngl(z(icol,ili))
+ write(20,*) sngl(x(icol,ili+istepz)),sngl(z(icol,ili+istepz))
+ write(20,10)
+ enddo
+ enddo
- 10 format('')
+10 format('')
- close(20)
+ close(20)
-! create a Gnuplot script to display the grid
- open(unit=20,file='OUTPUT_FILES/plotgnu',status='unknown')
- write(20,*) '#set term X11'
- write(20,*) 'set term postscript landscape monochrome solid "Helvetica" 22'
- write(20,*) 'set output "grid.ps"'
- write(20,*) '#set xrange [',sngl(minval(x)),':',sngl(maxval(x)),']'
- write(20,*) '#set yrange [',sngl(minval(z)),':',sngl(maxval(z)),']'
-! use same unit length on both X and Y axes
- write(20,*) 'set size ratio -1'
- write(20,*) 'plot "gridfile.gnu" title "Macrobloc mesh" w l'
- write(20,*) 'pause -1 "Hit any key..."'
- close(20)
+ ! create a Gnuplot script to display the grid
+ open(unit=20,file='OUTPUT_FILES/plotgnu',status='unknown')
+ write(20,*) '#set term X11'
+ write(20,*) 'set term postscript landscape monochrome solid "Helvetica" 22'
+ write(20,*) 'set output "grid.ps"'
+ write(20,*) '#set xrange [',sngl(minval(x)),':',sngl(maxval(x)),']'
+ write(20,*) '#set yrange [',sngl(minval(z)),':',sngl(maxval(z)),']'
+ ! use same unit length on both X and Y axes
+ write(20,*) 'set size ratio -1'
+ write(20,*) 'plot "gridfile.gnu" title "Macrobloc mesh" w l'
+ write(20,*) 'pause -1 "Hit any key..."'
+ close(20)
- print *,'Grid saved in Gnuplot format...'
- print *
+ print *,'Grid saved in Gnuplot format...'
+ print *
endif
-!*****************************
-! partitioning
-!*****************************
+ !*****************************
+ ! partitioning
+ !*****************************
allocate(part(0:nelmnts-1))
-! if ngnod == 9, we work on a subarray of elmnts, which represents the elements with for nodes only
-! construction of the graph
+ ! if ngnod == 9, we work on a subarray of elmnts, which represents the elements with for nodes only
+ ! construction of the graph
if ( ngnod == 9 ) then
allocate(elmnts_bis(0:ESIZE*nelmnts-1))
do i = 0, nelmnts-1
@@ -1323,27 +1353,27 @@
enddo
if ( nproc > 1 ) then
- call mesh2dual_ncommonnodes(nelmnts, (nxread+1)*(nzread+1), elmnts_bis, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
+ call mesh2dual_ncommonnodes(nelmnts, (nxread+1)*(nzread+1), elmnts_bis, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
endif
else
if ( nproc > 1 ) then
- call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
+ call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,1)
endif
endif
if ( nproc == 1 ) then
- part(:) = 0
+ part(:) = 0
else
- nb_edges = xadj(nelmnts)
+ nb_edges = xadj(nelmnts)
-! giving weight to edges and vertices. Currently not used.
- call read_weights(nelmnts, vwgt, nb_edges, adjwgt)
+ ! giving weight to edges and vertices. Currently not used.
+ call read_weights(nelmnts, vwgt, nb_edges, adjwgt)
-! partitioning
+ ! partitioning
select case (partitioning_method)
case(1)
@@ -1377,7 +1407,7 @@
endif
-! beware of fluid solid edges : coupled elements are transfered to the same partition
+ ! beware of fluid solid edges : coupled elements are transfered to the same partition
if ( ngnod == 9 ) then
call acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts_bis, nb_materials, phi, num_material, &
nproc, part, nedges_coupled, edges_coupled)
@@ -1385,7 +1415,7 @@
call acoustic_elastic_repartitioning (nelmnts, nnodes, elmnts, nb_materials, phi, num_material, &
nproc, part, nedges_coupled, edges_coupled)
endif
-! beware of fluid porous edges : coupled elements are transfered to the same partition
+ ! beware of fluid porous edges : coupled elements are transfered to the same partition
if ( ngnod == 9 ) then
call acoustic_poro_repartitioning (nelmnts, nnodes, elmnts_bis, nb_materials, phi, num_material, &
nproc, part, nedges_acporo_coupled, edges_acporo_coupled)
@@ -1393,7 +1423,7 @@
call acoustic_poro_repartitioning (nelmnts, nnodes, elmnts, nb_materials, phi, num_material, &
nproc, part, nedges_acporo_coupled, edges_acporo_coupled)
endif
-! beware of porous solid edges : coupled elements are transfered to the same partition
+ ! beware of porous solid edges : coupled elements are transfered to the same partition
if ( ngnod == 9 ) then
call poro_elastic_repartitioning (nelmnts, nnodes, elmnts_bis, nb_materials, phi, num_material, &
nproc, part, nedges_elporo_coupled, edges_elporo_coupled)
@@ -1402,13 +1432,13 @@
nproc, part, nedges_elporo_coupled, edges_elporo_coupled)
endif
-! local number of each element for each partition
+ ! local number of each element for each partition
call Construct_glob2loc_elmnts(nelmnts, part, nproc, glob2loc_elmnts)
if ( ngnod == 9 ) then
if ( nproc > 1 ) then
- deallocate(nnodes_elmnts)
- deallocate(nodes_elmnts)
+ deallocate(nnodes_elmnts)
+ deallocate(nodes_elmnts)
endif
allocate(nnodes_elmnts(0:nnodes-1))
allocate(nodes_elmnts(0:nsize*nnodes-1))
@@ -1421,25 +1451,25 @@
enddo
else
if ( nproc < 2 ) then
- allocate(nnodes_elmnts(0:nnodes-1))
- allocate(nodes_elmnts(0:nsize*nnodes-1))
- nnodes_elmnts(:) = 0
- nodes_elmnts(:) = 0
- do i = 0, ngnod*nelmnts-1
- nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
- nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
+ allocate(nnodes_elmnts(0:nnodes-1))
+ allocate(nodes_elmnts(0:nsize*nnodes-1))
+ nnodes_elmnts(:) = 0
+ nodes_elmnts(:) = 0
+ do i = 0, ngnod*nelmnts-1
+ nodes_elmnts(elmnts(i)*nsize+nnodes_elmnts(elmnts(i))) = i/ngnod
+ nnodes_elmnts(elmnts(i)) = nnodes_elmnts(elmnts(i)) + 1
- enddo
+ enddo
endif
endif
-! local number of each node for each partition
+ ! local number of each node for each partition
call Construct_glob2loc_nodes(nelmnts, nnodes, nnodes_elmnts, nodes_elmnts, part, nproc, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
-! construct the interfaces between partitions (used for MPI assembly)
+ ! construct the interfaces between partitions (used for MPI assembly)
if ( nproc /= 1 ) then
if ( ngnod == 9 ) then
call Construct_interfaces(nelmnts, nproc, part, elmnts_bis, xadj, adjncy, tab_interfaces, &
@@ -1452,7 +1482,7 @@
allocate(my_nb_interfaces(0:ninterfaces-1))
endif
-! setting absorbing boundaries by elements instead of edges
+ ! setting absorbing boundaries by elements instead of edges
if ( any_abs ) then
call merge_abs_boundaries(nelemabs, nelemabs_merge, abs_surface, abs_surface_char, abs_surface_merge, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
@@ -1462,7 +1492,7 @@
elmnts, ngnod)
endif
-! *** generate the databases for the solver
+ ! *** generate the databases for the solver
do iproc = 0, nproc-1
@@ -1510,8 +1540,8 @@
write(15,*) 'seismotype imagetype save_forward'
write(15,*) seismotype,imagetype,save_forward
- write(15,*) 'assign_external_model outputgrid TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
- write(15,*) assign_external_model,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ write(15,*) 'assign_external_model outputgrid TURN_ATTENUATION_ON'
+ write(15,*) assign_external_model,outputgrid,TURN_ATTENUATION_ON
write(15,*) 'TURN_VISCATTENUATION_ON Q0 freq0'
write(15,*) TURN_VISCATTENUATION_ON,Q0,freq0
@@ -1525,9 +1555,9 @@
write(15,*) NSOURCE
do i_source=1,NSOURCE
- write(15,*) 'source', i_source
- write(15,*) source_type(i_source),time_function_type(i_source),xs(i_source),zs(i_source),f0(i_source),t0(i_source), &
- factor(i_source),angleforce(i_source),Mxx(i_source),Mzz(i_source),Mxz(i_source)
+ write(15,*) 'source', i_source
+ write(15,*) source_type(i_source),time_function_type(i_source),xs(i_source),zs(i_source),f0(i_source),t0(i_source), &
+ factor(i_source),angleforce(i_source),Mxx(i_source),Mzz(i_source),Mxz(i_source)
enddo
write(15,*) 'attenuation'
@@ -1565,18 +1595,21 @@
write(15,*) 'nelemabs nelem_acoustic_surface num_fluid_solid_edges num_fluid_poro_edges'
write(15,*) 'num_solid_poro_edges nnodes_tangential_curve'
write(15,*) nelemabs_loc,nelem_acoustic_surface_loc,nedges_coupled_loc,nedges_acporo_coupled_loc,&
- nedges_elporo_coupled_loc,nnodes_tangential_curve
+ nedges_elporo_coupled_loc,nnodes_tangential_curve
write(15,*) 'Material sets (num 1 rho vp vs 0 0 Qp Qs 0 0 0 0 0 0) or '
write(15,*) '(num 2 rho c11 c13 c33 c44 Qp Qs 0 0 0 0 0 0) or '
write(15,*) '(num 3 rhos rhof phi c k_xx k_xz k_zz Ks Kf Kfr etaf mufr Qs)'
do i=1,nb_materials
- if (icodemat(i) /= POROELASTIC_MATERIAL)then
- write(15,*) i,icodemat(i),rho_s(i),cp(i),cs(i),aniso3(i),aniso4(i),Qp(i),Qs(i),0,0,0,0,0,0
- else
- write(15,*) i,icodemat(i),rho_s(i),rho_f(i),phi(i),tortuosity(i),permxx(i),permxz(i),permzz(i),kappa_s(i),&
- kappa_f(i),kappa_fr(i),eta_f(i),mu_fr(i),Qs(i)
- endif
+ if (icodemat(i) == ISOTROPIC_MATERIAL) then
+ write(15,*) i,icodemat(i),rho_s(i),cp(i),cs(i),0,0,Qp(i),Qs(i),0,0,0,0,0,0
+ elseif(icodemat(i) == POROELASTIC_MATERIAL) then
+ write(15,*) i,icodemat(i),rho_s(i),rho_f(i),phi(i),tortuosity(i),permxx(i),permxz(i),permzz(i),kappa_s(i),&
+ kappa_f(i),kappa_fr(i),eta_f(i),mu_fr(i),Qs(i)
+ else
+ write(15,*) i,icodemat(i),rho_s(i),cp(i),cs(i),aniso3(i),aniso4(i),aniso5(i),aniso6(i),&
+ aniso7(i),aniso8(i),Qp(i),Qs(i),0,0
+ endif
enddo
write(15,*) 'Arrays kmato and knods for each bloc:'
@@ -1633,82 +1666,82 @@
!write(15,*) nnodes_tangential_curve
write(15,*) force_normal_to_surface,rec_normal_to_surface
do i = 1, nnodes_tangential_curve
- write(15,*) nodes_tangential_curve(1,i),nodes_tangential_curve(2,i)
+ write(15,*) nodes_tangential_curve(1,i),nodes_tangential_curve(2,i)
enddo
enddo
-! print position of the source
+ ! print position of the source
do i_source=1,NSOURCE
print *
print *,'Position (x,z) of the source = ',xs(i_source),zs(i_source)
print *
enddo
-!--- compute position of the receivers and write the STATIONS file
+ !--- compute position of the receivers and write the STATIONS file
if (generate_STATIONS) then
- print *
- print *,'writing the DATA/STATIONS file'
- print *
+ print *
+ print *,'writing the DATA/STATIONS file'
+ print *
-! total number of receivers in all the receiver lines
- nrec_total = sum(nrec)
+ ! total number of receivers in all the receiver lines
+ nrec_total = sum(nrec)
- print *
- print *,'There are ',nrec_total,' receivers'
+ print *
+ print *,'There are ',nrec_total,' receivers'
- print *
- print *,'Position (x,z) of the ',nrec_total,' receivers'
- print *
+ print *
+ print *,'Position (x,z) of the ',nrec_total,' receivers'
+ print *
- open(unit=15,file='DATA/STATIONS',status='unknown')
+ open(unit=15,file='DATA/STATIONS',status='unknown')
- irec_global_number = 0
+ irec_global_number = 0
-! loop on all the receiver lines
- do ireceiverlines = 1,nreceiverlines
+ ! loop on all the receiver lines
+ do ireceiverlines = 1,nreceiverlines
-! loop on all the receivers of this receiver line
- do irec = 1,nrec(ireceiverlines)
+ ! loop on all the receivers of this receiver line
+ do irec = 1,nrec(ireceiverlines)
-! compute global receiver number
- irec_global_number = irec_global_number + 1
+ ! compute global receiver number
+ irec_global_number = irec_global_number + 1
-! compute coordinates of the receiver
- if(nrec(ireceiverlines) > 1) then
- xrec = xdeb(ireceiverlines) + dble(irec-1)*(xfin(ireceiverlines)-xdeb(ireceiverlines))/dble(nrec(ireceiverlines)-1)
- zrec = zdeb(ireceiverlines) + dble(irec-1)*(zfin(ireceiverlines)-zdeb(ireceiverlines))/dble(nrec(ireceiverlines)-1)
- else
- xrec = xdeb(ireceiverlines)
- zrec = zdeb(ireceiverlines)
- endif
+ ! compute coordinates of the receiver
+ if(nrec(ireceiverlines) > 1) then
+ xrec = xdeb(ireceiverlines) + dble(irec-1)*(xfin(ireceiverlines)-xdeb(ireceiverlines))/dble(nrec(ireceiverlines)-1)
+ zrec = zdeb(ireceiverlines) + dble(irec-1)*(zfin(ireceiverlines)-zdeb(ireceiverlines))/dble(nrec(ireceiverlines)-1)
+ else
+ xrec = xdeb(ireceiverlines)
+ zrec = zdeb(ireceiverlines)
+ endif
-! modify position of receiver if we must record exactly at the surface
- if(enreg_surf(ireceiverlines)) &
- zrec = value_spline(xrec,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
+ ! modify position of receiver if we must record exactly at the surface
+ if(enreg_surf(ireceiverlines)) &
+ zrec = value_spline(xrec,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
-! display position of the receiver
- print *,'Receiver ',irec_global_number,' = ',xrec,zrec
+ ! display position of the receiver
+ print *,'Receiver ',irec_global_number,' = ',xrec,zrec
- write(15,"('S',i4.4,' AA ',f20.7,1x,f20.7,' 0.0 0.0')") irec_global_number,xrec,zrec
+ write(15,"('S',i4.4,' AA ',f20.7,1x,f20.7,' 0.0 0.0')") irec_global_number,xrec,zrec
- enddo
- enddo
+ enddo
+ enddo
- close(15)
+ close(15)
endif
print *
if (nproc == 1) then
- print *,'This will be a serial simulation'
+ print *,'This will be a serial simulation'
else
- print *,'This will be a parallel simulation on ',nproc,' processors'
+ print *,'This will be a parallel simulation on ',nproc,' processors'
endif
print *
- end program meshfem2D
+end program meshfem2D
! *******************
! meshing subroutines
@@ -1717,54 +1750,54 @@
!--- global node number
- integer function num(i,j,nx)
+integer function num(i,j,nx)
- implicit none
+ implicit none
- integer i,j,nx
+ integer i,j,nx
- num = j*(nx+1) + i + 1
+ num = j*(nx+1) + i + 1
- end function num
+end function num
- !--- global node number (when ngnod==4).
- integer function num_4(i,j,nx)
+!--- global node number (when ngnod==4).
+integer function num_4(i,j,nx)
- implicit none
+ implicit none
- integer i,j,nx
+ integer i,j,nx
- num_4 = j*(nx+1) + i + 1
+ num_4 = j*(nx+1) + i + 1
- end function num_4
+end function num_4
- !--- global node number (when ngnod==9).
- integer function num_9(i,j,nx,nz)
+!--- global node number (when ngnod==9).
+integer function num_9(i,j,nx,nz)
- implicit none
+ implicit none
- integer i,j,nx,nz
+ integer i,j,nx,nz
- if ( (mod(i,2) == 0) .and. (mod(j,2) == 0) ) then
- num_9 = j/2 * (nx+1) + i/2 + 1
- else
- if ( mod(j,2) == 0 ) then
- num_9 = (nx+1)*(nz+1) + j/2 * nx + ceiling(real(i)/real(2))
- else
- num_9 = (nx+1)*(nz+1) + nx*(nz+1) + floor(real(j)/real(2))*(nx*2+1) + i + 1
+ if ( (mod(i,2) == 0) .and. (mod(j,2) == 0) ) then
+ num_9 = j/2 * (nx+1) + i/2 + 1
+ else
+ if ( mod(j,2) == 0 ) then
+ num_9 = (nx+1)*(nz+1) + j/2 * nx + ceiling(real(i)/real(2))
+ else
+ num_9 = (nx+1)*(nz+1) + nx*(nz+1) + floor(real(j)/real(2))*(nx*2+1) + i + 1
- endif
- endif
+ endif
+ endif
- end function num_9
+end function num_9
!--- spline to describe the interfaces
- double precision function value_spline(x,xinterface,zinterface,coefs_interface,npoints_interface)
+double precision function value_spline(x,xinterface,zinterface,coefs_interface,npoints_interface)
implicit none
@@ -1776,11 +1809,11 @@
xp = x
-! assign the value on the edge if point is outside the model
+ ! assign the value on the edge if point is outside the model
if(xp < xinterface(1)) xp = xinterface(1)
if(xp > xinterface(npoints_interface)) xp = xinterface(npoints_interface)
call spline_evaluation(xinterface,zinterface,coefs_interface,npoints_interface,xp,value_spline)
- end function value_spline
+end function value_spline
Modified: seismo/2D/SPECFEM2D/trunk/part_unstruct.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/part_unstruct.F90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -1032,6 +1032,7 @@
jend_left(:) = NGLLZ
is_acoustic(:) = .false.
+
do i = 1, nb_materials
if (phi_material(i) >= 1.d0) then
is_acoustic(i) = .true.
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.F90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.F90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -3066,8 +3066,8 @@
620 format(f6.3,' neg CM 0 MR (Cut =',f5.2,' \%) show')
640 format(f6.3,' neg CM 0 MR (Max norm =',1pe12.3,') show')
- 499 format(f6.2,1x,f6.2,' L')
- 500 format(f6.2,1x,f6.2,' M')
+ 499 format(f8.3,1x,f8.3,' L')
+ 500 format(f8.3,1x,f8.3,' M')
502 format('fN (',i4,') Cshow')
679 format(f12.6,1x,f12.6,1x,f12.6,' RG fill stroke')
680 format(f12.6,1x,f12.6,1x,f12.6,' RG GF')
Added: seismo/2D/SPECFEM2D/trunk/read_external_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_external_model.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/read_external_model.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -0,0 +1,170 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.0
+! ------------------------------
+!
+! Copyright Universite de Pau et des Pays de l'Adour, CNRS and INRIA, France,
+! and Princeton University, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+! Pieyre Le Loher, pieyre DOT le-loher aT inria.fr
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+
+subroutine read_external_model(any_acoustic,any_elastic,any_poroelastic, &
+ elastic,poroelastic,anisotropic,nspec,npoin,N_SLS,ibool, &
+ f0_attenuation,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,myrank,rhoext,vpext,vsext, &
+ Qp_attenuationext,Qs_attenuationext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext)
+
+ implicit none
+ include "constants.h"
+
+integer :: nspec,myrank,npoin
+double precision :: f0_attenuation
+
+! Mesh
+integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+double precision, dimension(NDIM,npoin) :: coord
+
+! Material properties
+logical :: any_acoustic,any_elastic,any_poroelastic
+integer, dimension(nspec) :: kmato
+logical, dimension(nspec) :: elastic,poroelastic
+double precision, dimension(NGLLX,NGLLZ,nspec) :: rhoext,vpext,vsext
+
+! for attenuation
+integer :: N_SLS
+double precision :: Mu_nu1_sent,Mu_nu2_sent
+double precision, dimension(N_SLS) :: inv_tau_sigma_nu1_sent,phi_nu1_sent,inv_tau_sigma_nu2_sent,phi_nu2_sent
+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
+double precision, dimension(NGLLX,NGLLZ,nspec) :: Qp_attenuationext,Qs_attenuationext
+
+! for anisotropy
+logical, dimension(nspec) :: anisotropic
+double precision, dimension(NGLLX,NGLLZ,nspec) :: c11ext,c13ext,c15ext,c33ext,c35ext,c55ext
+
+! Local variables
+integer :: i,j,ispec,iglob
+double precision :: previous_vsext
+double precision :: tmp1, tmp2,tmp3
+
+if(READ_EXTERNAL_MODEL_FILE) then
+ write(IOUT,*)
+ write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+ write(IOUT,*) 'Assigning external velocity and density model (elastic and/or acoustic)...'
+ write(IOUT,*) 'Read outside SEG model...'
+ write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
+
+ open(unit=1001,file='DATA/model_velocity.dat_input',status='unknown')
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec)
+ read(1001,*) tmp1,tmp2,tmp3,rhoext(i,j,ispec),vpext(i,j,ispec),vsext(i,j,ispec)
+! vsext(i,j,ispec)=0.0
+ end do
+ end do
+ end do
+ close(1001)
+
+ else
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+ iglob = ibool(i,j,ispec)
+ call define_external_model(coord(1,iglob),coord(2,iglob),kmato(ispec),myrank,&
+ rhoext(i,j,ispec),vpext(i,j,ispec),vsext(i,j,ispec), &
+ Qp_attenuationext(i,j,ispec),Qs_attenuationext(i,j,ispec),&
+ c11ext(i,j,ispec),c13ext(i,j,ispec),c15ext(i,j,ispec),c33ext(i,j,ispec),c35ext(i,j,ispec),c55ext(i,j,ispec))
+ if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
+ (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
+ ! vp, vs : dummy values, trick to avoid floating point errors
+ vpext(i,j,ispec) = 20.d0
+ vsext(i,j,ispec) = 10.d0
+ end if
+ end do
+ end do
+ end do
+ end if
+
+ any_acoustic = .false.
+ any_elastic = .false.
+ any_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))) &
+ call exit_MPI('external velocity model cannot be both fluid and solid inside the same spectral element')
+ if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
+ (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
+ anisotropic(ispec) = .true.
+ poroelastic(ispec) = .false.
+ elastic(ispec) = .true.
+ any_elastic = .true.
+ Qp_attenuationext(i,j,ispec) = 10.d0
+ Qs_attenuationext(i,j,ispec) = 10.d0
+ elseif(vsext(i,j,ispec) < TINYVAL) then
+ elastic(ispec) = .false.
+ poroelastic(ispec) = .false.
+ any_acoustic = .true.
+ else
+ poroelastic(ispec) = .false.
+ elastic(ispec) = .true.
+ any_elastic = .true.
+ endif
+ call attenuation_model(N_SLS,Qp_attenuationext(i,j,ispec),Qs_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)
+ 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
+ previous_vsext = vsext(i,j,ispec)
+ enddo
+ enddo
+ enddo
+
+end subroutine read_external_model
Modified: seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -156,6 +156,7 @@
character(len=100) string_read
call read_next_line(iin,ignore_junk,string_read)
+ print*,string_read
read(string_read,*) i,icodematread,val0read,val1read,val2read,val3read,val4read,val5read,&
val6read,val7read,val8read,val9read,val10read,val11read,val12read
Added: seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/setup_sources_receivers.f90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -0,0 +1,164 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 6.0
+! ------------------------------
+!
+! Copyright Universite de Pau et des Pays de l'Adour, CNRS and INRIA, France,
+! and Princeton University, USA.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+! Christina Morency, cmorency aT princeton DOT edu
+! Pieyre Le Loher, pieyre DOT le-loher aT inria.fr
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+subroutine setup_sources_receivers(NSOURCE,assign_external_model,initialfield,numat,source_type,&
+ coord,ibool,kmato,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
+ x_source,z_source,ix_source,iz_source,ispec_selected_source,ispec_selected_rec,iglob_source, &
+ is_proc_source,nb_proc_source,rho_at_source_location,ipass,&
+ sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,density,rhoext,&
+ nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod, &
+ nrec,nrecloc,recloc,which_proc_receiver,st_xval,st_zval, &
+ xi_receiver,gamma_receiver,station_name,network_name,x_final_receiver,z_final_receiver)
+
+ implicit none
+
+ include "constants.h"
+
+ logical :: assign_external_model,initialfield
+ integer :: NSOURCE, numat
+ integer :: npgeo,ngnod,myrank,ipass,nproc
+ integer :: npoin,nspec,nelem_acoustic_surface
+
+ ! Gauss-Lobatto-Legendre points
+ double precision, dimension(NGLLX) :: xigll
+ double precision, dimension(NGLLZ) :: zigll
+
+ ! for receivers
+ integer :: nrec,nrecloc
+ integer, dimension(nrec) :: recloc, which_proc_receiver
+ integer, dimension(nrec) :: ispec_selected_rec
+ double precision, dimension(nrec) :: xi_receiver,gamma_receiver,st_xval,st_zval
+ double precision, dimension(nrec) :: x_final_receiver, z_final_receiver
+
+ ! timing information for the stations
+ character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
+ character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
+
+ ! for sources
+ double precision :: rho_at_source_location
+ integer, dimension(NSOURCE) :: source_type
+ integer, dimension(NSOURCE) :: ispec_selected_source,iglob_source,ix_source,iz_source,is_proc_source,nb_proc_source
+ real(kind=CUSTOM_REAL), dimension(NSOURCE,NDIM,NGLLX,NGLLZ) :: sourcearray
+ double precision, dimension(NSOURCE) :: x_source,z_source,xi_source,gamma_source,Mxx,Mzz,Mxz
+
+ logical, dimension(nspec) :: elastic,poroelastic
+ integer, dimension(nspec) :: kmato
+ integer, dimension(ngnod,nspec) :: knods
+ integer, dimension(5,nelem_acoustic_surface) :: acoustic_surface
+ integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
+ double precision, dimension(NDIM,npgeo) :: coorg
+ double precision, dimension(2,numat) :: density
+ double precision, dimension(NGLLX,NGLLZ,nspec) :: rhoext
+ double precision, dimension(NDIM,npoin) :: coord
+
+ ! Local variables
+ integer i_source,ispec,ispec_acoustic_surface,i,j,iglob
+
+ do i_source=1,NSOURCE
+
+ if(source_type(i_source) == 1) then
+
+ ! collocated force source
+ call locate_source_force(coord,ibool,npoin,nspec,x_source(i_source),z_source(i_source), &
+ ix_source(i_source),iz_source(i_source),ispec_selected_source(i_source),iglob_source(i_source), &
+ is_proc_source(i_source),nb_proc_source(i_source),ipass)
+
+ ! get density at the source in order to implement collocated force with the right
+ ! amplitude later
+ if(is_proc_source(i_source) == 1) then
+ rho_at_source_location = density(1,kmato(ispec_selected_source(i_source)))
+ ! external velocity model
+ if(assign_external_model) rho_at_source_location = &
+ rhoext(ix_source(i_source),iz_source(i_source),ispec_selected_source(i_source))
+ endif
+
+ ! check that acoustic source is not exactly on the free surface because pressure is zero there
+ if(is_proc_source(i_source) == 1) then
+ do ispec_acoustic_surface = 1,nelem_acoustic_surface
+ ispec = acoustic_surface(1,ispec_acoustic_surface)
+ if( .not. elastic(ispec) .and. .not. poroelastic(ispec) .and. ispec == ispec_selected_source(i_source) ) then
+ do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
+ do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
+ iglob = ibool(i,j,ispec)
+ if ( iglob_source(i_source) == iglob ) then
+
+ call exit_MPI('an acoustic source cannot be located exactly on the free surface because pressure is zero there')
+
+ endif
+ enddo
+ enddo
+ endif
+ enddo
+ endif
+
+ else if(source_type(i_source) == 2) then
+ ! moment-tensor source
+ call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
+ ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source),&
+ nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo,ipass)
+
+ ! compute source array for moment-tensor source
+ call compute_arrays_source(ispec_selected_source(i_source),xi_source(i_source),gamma_source(i_source),&
+ sourcearray(i_source,:,:,:), &
+ Mxx(i_source),Mzz(i_source),Mxz(i_source),xix,xiz,gammax,gammaz,xigll,zigll,nspec)
+
+ else if(.not.initialfield) then
+ call exit_MPI('incorrect source type')
+ endif
+
+
+ ! locate receivers in the mesh
+ call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
+ st_xval,st_zval,ispec_selected_rec, &
+ xi_receiver,gamma_receiver,station_name,network_name,x_source(i_source),z_source(i_source),&
+ coorg,knods,ngnod,npgeo,ipass, &
+ x_final_receiver, z_final_receiver)
+
+ enddo ! do i_source=1,NSOURCE
+
+
+end subroutine setup_sources_receivers
+
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -210,6 +210,7 @@
! 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
@@ -289,6 +290,12 @@
! poroelastic and elastic coefficients
double precision, dimension(:,:,:), allocatable :: poroelastcoef
+! anisotropy parameters
+ logical :: all_anisotropic
+ double precision :: c11,c13,c15,c33,c35,c55
+ logical, dimension(:), allocatable :: anisotropic
+ double precision, dimension(:,:), allocatable :: anisotropy
+
! for acoustic medium
real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
@@ -308,8 +315,10 @@
double precision, dimension(:), allocatable :: vp_display
- double precision, dimension(:,:,:), allocatable :: vpext,vsext,rhoext
- double precision :: previous_vsext,rho_at_source_location
+ double precision, dimension(:,:,:), allocatable :: vpext,vsext,rhoext
+ double precision :: rho_at_source_location
+ double precision, dimension(:,:,:), allocatable :: Qp_attenuationext,Qs_attenuationext
+ double precision, dimension(:,:,:), allocatable :: c11ext,c13ext,c15ext,c33ext,c35ext,c55ext
double precision, dimension(:,:,:), allocatable :: shape2D,shape2D_display
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: xix,xiz,gammax,gammaz,jacobian
@@ -334,7 +343,7 @@
integer :: numat,ngnod,nspec,pointsdisp,nelemabs,nelem_acoustic_surface,ispecabs,UPPER_LIMIT_DISPLAY
logical interpol,meshvect,modelvect,boundvect,assign_external_model,initialfield, &
- outputgrid,gnuplot,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,output_postscript_snapshot,output_color_image, &
+ outputgrid,gnuplot,TURN_ATTENUATION_ON,output_postscript_snapshot,output_color_image, &
plot_lowerleft_corner_only,add_Bielak_conditions
double precision :: cutsnaps,sizemax_arrows,anglerec,xirec,gammarec
@@ -629,7 +638,6 @@
double precision :: x_final_receiver_dummy, z_final_receiver_dummy
!!!!!!!!!!
double precision, dimension(:,:,:),allocatable:: rho_local,vp_local,vs_local
- double precision :: tmp1, tmp2,tmp3
!!!! hessian
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhorho_el_hessian_final1, rhorho_el_hessian_temp1
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rhorho_el_hessian_final2, rhorho_el_hessian_temp2
@@ -759,7 +767,7 @@
endif
read(IIN,"(a80)") datlin
- read(IIN,*) assign_external_model,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ read(IIN,*) assign_external_model,outputgrid,TURN_ATTENUATION_ON
read(IIN,"(a80)") datlin
read(IIN,*) TURN_VISCATTENUATION_ON,Q0,freq0
@@ -772,7 +780,7 @@
write(IOUT,200) npgeo,NDIM
write(IOUT,600) NTSTEP_BETWEEN_OUTPUT_INFO,colors,numbers
write(IOUT,700) seismotype,anglerec
- write(IOUT,750) initialfield,add_Bielak_conditions,assign_external_model,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
+ write(IOUT,750) initialfield,add_Bielak_conditions,assign_external_model,TURN_ATTENUATION_ON,outputgrid
write(IOUT,800) imagetype,100.d0*cutsnaps,subsamp
endif
@@ -781,7 +789,7 @@
read(IIN,*) NSTEP,deltat,isolver
if (myrank == 0 .and. ipass == 1) write(IOUT,703) NSTEP,deltat,NSTEP*deltat
- if(isolver == 1 .and. save_forward .and. (TURN_ANISOTROPY_ON .or. TURN_ATTENUATION_ON .or. TURN_VISCATTENUATION_ON)) then
+ if(isolver == 1 .and. save_forward .and. (TURN_ATTENUATION_ON .or. TURN_VISCATTENUATION_ON)) then
print*, '*************** WARNING ***************'
print*, 'Anisotropy & Attenuation & Viscous damping are not presently implemented for adjoint calculations'
stop
@@ -881,7 +889,7 @@
ipoin = 0
read(IIN,"(a80)") datlin
- allocate(coorgread(NDIM))
+ allocate(coorgread(NDIM))
do ip = 1,npgeo
read(IIN,*) ipoin,(coorgread(id),id =1,NDIM)
if(ipoin<1 .or. ipoin>npgeo) call exit_MPI('Wrong control point number')
@@ -918,6 +926,7 @@
allocate(Uxinterp(pointsdisp,pointsdisp))
allocate(Uzinterp(pointsdisp,pointsdisp))
allocate(density(2,numat))
+ allocate(anisotropy(6,numat))
allocate(porosity(numat))
allocate(tortuosity(numat))
allocate(permeability(3,numat))
@@ -929,6 +938,7 @@
allocate(ibool(NGLLX,NGLLZ,nspec))
allocate(elastic(nspec))
allocate(poroelastic(nspec))
+ allocate(anisotropic(nspec))
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))
@@ -940,12 +950,14 @@
endif
! --- allocate arrays for absorbing boundary conditions
+
if(nelemabs <= 0) then
nelemabs = 1
anyabs = .false.
- else
+ else
anyabs = .true.
endif
+
if(ipass == 1) then
allocate(numabs(nelemabs))
allocate(codeabs(4,nelemabs))
@@ -985,7 +997,7 @@
!
!---- read the material properties
!
- call gmat01(density,porosity,tortuosity,permeability,poroelastcoef,numat,&
+ call gmat01(density,porosity,tortuosity,anisotropy,permeability,poroelastcoef,numat,&
myrank,ipass,Qp_attenuation,Qs_attenuation)
!
!---- read spectral macrobloc data
@@ -1013,6 +1025,7 @@
any_acoustic = .false.
any_elastic = .false.
any_poroelastic = .false.
+ anisotropic(:) = .false.
do ispec = 1,nspec
if(porosity(kmato(ispec)) == 1.d0) then ! acoustic domain
@@ -1023,6 +1036,9 @@
elastic(ispec) = .true.
poroelastic(ispec) = .false.
any_elastic = .true.
+ if(any(anisotropy(:,kmato(ispec)) /= 0)) then
+ anisotropic(ispec) = .true.
+ end if
else ! poroelastic domain
elastic(ispec) = .false.
poroelastic(ispec) = .true.
@@ -1038,12 +1054,13 @@
print*, '*************** WARNING ***************'
stop
endif
- if(.not. p_sv .and. (TURN_ATTENUATION_ON .or. TURN_ANISOTROPY_ON)) then
+ if(.not. p_sv .and. (TURN_ATTENUATION_ON)) then
print*, '*************** WARNING ***************'
print*, 'Attenuation and anisotropy are not implemented for surface (membrane) waves calculation'
print*, '*************** WARNING ***************'
stop
endif
+
if(TURN_ATTENUATION_ON) then
nspec_allocate = nspec
@@ -1088,7 +1105,7 @@
Mu_nu2(i,j,ispec) = Mu_nu2_sent
enddo
enddo
- enddo
+ enddo
! allocate memory variables for viscous attenuation (poroelastic media)
if(ipass == 1) then
@@ -1169,7 +1186,7 @@
if (myrank == 0 .and. ipass == 1) then
write(IOUT,*)
write(IOUT,*) 'Number of absorbing elements: ',nelemabs
- endif
+ endif
nspec_xmin = ZERO
nspec_xmax = ZERO
@@ -1250,7 +1267,29 @@
write(IOUT,*) 'nspec_zmin = ',nspec_zmin
write(IOUT,*) 'nspec_zmax = ',nspec_zmax
- endif
+ elseif(ipass == 1) then
+ allocate(ib_xmin(1))
+ allocate(ib_xmax(1))
+ allocate(ib_zmin(1))
+ allocate(ib_zmax(1))
+ allocate(b_absorb_elastic_left(1,1,1,1))
+ allocate(b_absorb_elastic_right(1,1,1,1))
+ allocate(b_absorb_elastic_bottom(1,1,1,1))
+ allocate(b_absorb_elastic_top(1,1,1,1))
+ allocate(b_absorb_poro_s_left(1,1,1,1))
+ allocate(b_absorb_poro_s_right(1,1,1,1))
+ allocate(b_absorb_poro_s_bottom(1,1,1,1))
+ allocate(b_absorb_poro_s_top(1,1,1,1))
+ allocate(b_absorb_poro_w_left(1,1,1,1))
+ allocate(b_absorb_poro_w_right(1,1,1,1))
+ allocate(b_absorb_poro_w_bottom(1,1,1,1))
+ allocate(b_absorb_poro_w_top(1,1,1,1))
+ allocate(b_absorb_acoustic_left(1,1,1))
+ allocate(b_absorb_acoustic_right(1,1,1))
+ allocate(b_absorb_acoustic_bottom(1,1,1))
+ allocate(b_absorb_acoustic_top(1,1,1))
+
+ endif
!
!---- read acoustic free surface data
!
@@ -1577,10 +1616,24 @@
allocate(vpext(NGLLX,NGLLZ,nspec))
allocate(vsext(NGLLX,NGLLZ,nspec))
allocate(rhoext(NGLLX,NGLLZ,nspec))
+ allocate(Qp_attenuationext(NGLLX,NGLLZ,nspec))
+ allocate(Qs_attenuationext(NGLLX,NGLLZ,nspec))
+ allocate(c11ext(NGLLX,NGLLZ,nspec))
+ allocate(c13ext(NGLLX,NGLLZ,nspec))
+ allocate(c15ext(NGLLX,NGLLZ,nspec))
+ allocate(c33ext(NGLLX,NGLLZ,nspec))
+ allocate(c35ext(NGLLX,NGLLZ,nspec))
+ allocate(c55ext(NGLLX,NGLLZ,nspec))
else
allocate(vpext(1,1,1))
allocate(vsext(1,1,1))
allocate(rhoext(1,1,1))
+ allocate(c11ext(1,1,1))
+ allocate(c13ext(1,1,1))
+ allocate(c15ext(1,1,1))
+ allocate(c33ext(1,1,1))
+ allocate(c35ext(1,1,1))
+ allocate(c55ext(1,1,1))
endif
endif
@@ -1676,53 +1729,24 @@
!
if(gnuplot .and. myrank == 0 .and. ipass == 1) call plotgll(knods,ibool,coorg,coord,npoin,npgeo,ngnod,nspec)
- write(IOUT,*) 'assign_external_model = ', assign_external_model
-if ( assign_external_model .and. ipass == 1) then
- write(IOUT,*)
- write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
- write(IOUT,*) 'Assigning external velocity and density model (elastic and/or acoustic)...'
- write(IOUT,*) 'Read outside SEG model...'
- write(IOUT,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
- if(TURN_ANISOTROPY_ON .or. TURN_ATTENUATION_ON) &
- call exit_MPI('cannot have anisotropy nor attenuation if external model in current version')
- any_acoustic = .false.
- any_elastic = .false.
- any_poroelastic = .false.
- open(unit=1001,file='DATA/model_velocity.dat_input',status='unknown')
- do ispec = 1,nspec
- do j = 1,NGLLZ
- do i = 1,NGLLX
- iglob = ibool(i,j,ispec)
- read(1001,*) tmp1,tmp2,tmp3,rhoext(i,j,ispec),vpext(i,j,ispec),vsext(i,j,ispec)
-! vsext(i,j,ispec)=0.0
- end do
- end do
- end do
- close(1001)
- 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))) &
- call exit_MPI('external velocity model cannot be both fluid and solid inside the same spectral element')
- if(vsext(i,j,ispec) < TINYVAL) then
- elastic(ispec) = .false.
- poroelastic(ispec) = .false.
- any_acoustic = .true.
- else
- poroelastic(ispec) = .false.
- elastic(ispec) = .true.
- any_elastic = .true.
- endif
- previous_vsext = vsext(i,j,ispec)
- enddo
- enddo
- enddo
-endif
+ if(myrank == 0 .and. ipass == 1) write(IOUT,*) 'assign_external_model = ', assign_external_model
+!if ( assign_external_model .and. ipass == 1) then
+ if ( assign_external_model) then
+ call read_external_model(any_acoustic,any_elastic,any_poroelastic, &
+ elastic,poroelastic,anisotropic,nspec,npoin,N_SLS,ibool, &
+ f0_attenuation,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,myrank,rhoext,vpext,vsext, &
+ Qp_attenuationext,Qs_attenuationext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext)
+ end if
+ if(count(anisotropic(:) .eqv. .true.) == nspec) all_anisotropic = .true.
+ if(all_anisotropic .and. anyabs) stop 'Cannot put absorbing boundaries if anisotropic materials along edges'
+ if(TURN_ATTENUATION_ON .and. all_anisotropic) then
+ stop 'Cannot turn attenuation on in anisotropic materials'
+ end if
+
!
!---- perform basic checks on parameters read
!
@@ -1742,16 +1766,9 @@
#endif
! for acoustic
- if(TURN_ANISOTROPY_ON .and. .not. any_elastic_glob) &
- call exit_MPI('cannot have anisotropy if acoustic/poroelastic simulation only')
-
if(TURN_ATTENUATION_ON .and. .not. any_elastic_glob) &
call exit_MPI('currently cannot have attenuation if acoustic/poroelastic simulation only')
-! for attenuation
- if(TURN_ANISOTROPY_ON .and. TURN_ATTENUATION_ON) &
- call exit_MPI('cannot have anisotropy and attenuation both turned on in current version')
-
!
!---- define coefficients of the Newmark time scheme
!
@@ -1766,65 +1783,16 @@
endif
!---- define actual location of source and receivers
- do i_source=1,NSOURCE
- if(source_type(i_source) == 1) then
+call setup_sources_receivers(NSOURCE,assign_external_model,initialfield,numat,source_type,&
+ coord,ibool,kmato,npoin,nspec,nelem_acoustic_surface,acoustic_surface,elastic,poroelastic, &
+ x_source,z_source,ix_source,iz_source,ispec_selected_source,ispec_selected_rec,iglob_source, &
+ is_proc_source,nb_proc_source,rho_at_source_location,ipass,&
+ sourcearray,Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,npgeo,density,rhoext,&
+ nproc,myrank,xi_source,gamma_source,coorg,knods,ngnod, &
+ nrec,nrecloc,recloc,which_proc_receiver,st_xval,st_zval, &
+ xi_receiver,gamma_receiver,station_name,network_name,x_final_receiver,z_final_receiver)
-! collocated force source
- call locate_source_force(coord,ibool,npoin,nspec,x_source(i_source),z_source(i_source), &
- ix_source(i_source),iz_source(i_source),ispec_selected_source(i_source),iglob_source(i_source), &
- is_proc_source(i_source),nb_proc_source(i_source),ipass)
-
-! get density at the source in order to implement collocated force with the right
-! amplitude later
- if(is_proc_source(i_source) == 1) then
- rho_at_source_location = density(1,kmato(ispec_selected_source(i_source)))
-! external velocity model
- if(assign_external_model) rho_at_source_location = &
- rhoext(ix_source(i_source),iz_source(i_source),ispec_selected_source(i_source))
- endif
-
-! check that acoustic source is not exactly on the free surface because pressure is zero there
- if(is_proc_source(i_source) == 1) then
- do ispec_acoustic_surface = 1,nelem_acoustic_surface
- ispec = acoustic_surface(1,ispec_acoustic_surface)
- if( .not. elastic(ispec) .and. .not. poroelastic(ispec) .and. ispec == ispec_selected_source(i_source) ) then
- do j = acoustic_surface(4,ispec_acoustic_surface), acoustic_surface(5,ispec_acoustic_surface)
- do i = acoustic_surface(2,ispec_acoustic_surface), acoustic_surface(3,ispec_acoustic_surface)
- iglob = ibool(i,j,ispec)
- if ( iglob_source(i_source) == iglob ) then
- call exit_MPI('an acoustic source cannot be located exactly on the free surface because pressure is zero there')
- endif
- enddo
- enddo
- endif
- enddo
- endif
-
- else if(source_type(i_source) == 2) then
-! moment-tensor source
- call locate_source_moment_tensor(ibool,coord,nspec,npoin,xigll,zigll,x_source(i_source),z_source(i_source), &
- ispec_selected_source(i_source),is_proc_source(i_source),nb_proc_source(i_source),&
- nproc,myrank,xi_source(i_source),gamma_source(i_source),coorg,knods,ngnod,npgeo,ipass)
-
-! compute source array for moment-tensor source
- call compute_arrays_source(ispec_selected_source(i_source),xi_source(i_source),gamma_source(i_source),&
- sourcearray(i_source,:,:,:), &
- Mxx(i_source),Mzz(i_source),Mxz(i_source),xix,xiz,gammax,gammaz,xigll,zigll,nspec)
-
- else if(.not.initialfield) then
- call exit_MPI('incorrect source type')
- endif
-
-
-! locate receivers in the mesh
- call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
- st_xval,st_zval,ispec_selected_rec, &
- xi_receiver,gamma_receiver,station_name,network_name,x_source(i_source),z_source(i_source),coorg,knods,ngnod,npgeo,ipass, &
- x_final_receiver, z_final_receiver)
-
- enddo ! do i_source=1,NSOURCE
-
! compute source array for adjoint source
if(isolver == 2) then ! adjoint calculation
nadj_rec_local = 0
@@ -2390,7 +2358,7 @@
wxgll(i)*wzgll(j)*jacobian(i,j,ispec)*(rhol_bar*rhol_f*tortl - &
phil*rhol_f*rhol_f)/(rhol_bar*phil)
elseif(elastic(ispec)) then ! for elastic medium
- rmass_inverse_elastic(iglob) = rmass_inverse_elastic(iglob) + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
+ rmass_inverse_elastic(iglob) = rmass_inverse_elastic(iglob) + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
else ! for acoustic medium
rmass_inverse_acoustic(iglob) = rmass_inverse_acoustic(iglob) + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / kappal
endif
@@ -2653,18 +2621,9 @@
!--- end of section performed in two passes
!---
-! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
- if(any_elastic) where(rmass_inverse_elastic <= 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
- if(any_poroelastic) where(rmass_s_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_s_inverse_poroelastic = 1._CUSTOM_REAL
- if(any_poroelastic) where(rmass_w_inverse_poroelastic <= 0._CUSTOM_REAL) rmass_w_inverse_poroelastic = 1._CUSTOM_REAL
- if(any_acoustic) where(rmass_inverse_acoustic <= 0._CUSTOM_REAL) rmass_inverse_acoustic = 1._CUSTOM_REAL
+call invert_mass_matrix(any_elastic,any_acoustic,any_poroelastic,npoin,rmass_inverse_elastic,&
+ rmass_inverse_acoustic,rmass_s_inverse_poroelastic,rmass_w_inverse_poroelastic)
-! compute the inverse of the mass matrix
- if(any_elastic) rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
- if(any_poroelastic) rmass_s_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_s_inverse_poroelastic(:)
- if(any_poroelastic) rmass_w_inverse_poroelastic(:) = 1._CUSTOM_REAL / rmass_w_inverse_poroelastic(:)
- if(any_acoustic) rmass_inverse_acoustic(:) = 1._CUSTOM_REAL / rmass_inverse_acoustic(:)
-
! check the mesh, stability and number of points per wavelength
if(DISPLAY_SUBSET_OPTION == 1) then
UPPER_LIMIT_DISPLAY = nspec
@@ -2677,11 +2636,12 @@
else
stop 'incorrect value of DISPLAY_SUBSET_OPTION'
endif
+
call checkgrid(vpext,vsext,rhoext,density,poroelastcoef,porosity,tortuosity,ibool,kmato, &
coord,npoin,vpImin,vpImax,vpIImin,vpIImax, &
assign_external_model,nspec,UPPER_LIMIT_DISPLAY,numat,deltat,f0,t0,initialfield, &
time_function_type,coorg,xinterp,zinterp,shape2D_display,knods,simulation_title, &
- npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,myrank,nproc,NSOURCE,poroelastic)
+ npgeo,pointsdisp,ngnod,any_elastic,any_poroelastic,all_anisotropic,myrank,nproc,NSOURCE,poroelastic)
! convert receiver angle to radians
anglerec = anglerec * pi / 180.d0
@@ -2753,6 +2713,7 @@
iglob_image_color(:,:) = -1
! cheking which pixels are inside each elements
+
nb_pixel_loc = 0
do ispec = 1, nspec
@@ -2865,7 +2826,7 @@
if (myrank == 0) write(IOUT,*) 'done locating all the pixels of color images'
- endif
+endif
!
!---- initialize seismograms
@@ -5243,11 +5204,12 @@
call compute_forces_elastic(p_sv,npoin,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
- initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
+ initialfield,TURN_ATTENUATION_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic,b_accel_elastic,b_displ_elastic, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
- jacobian,vpext,vsext,rhoext,source_time_function,sourcearray,adj_sourcearrays, &
+ jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
+ source_time_function,sourcearray,adj_sourcearrays, &
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,N_SLS, &
@@ -5634,11 +5596,27 @@
endif
! compute stress tensor
! full anisotropy
- if(TURN_ANISOTROPY_ON) then
+ if(kmato(ispec_elastic) == 2) then
! implement anisotropy in 2D
- sigma_xx = sigma_xx + c11val*dux_dxl + c15val*(duz_dxl + dux_dzl) + c13val*duz_dzl
- sigma_zz = sigma_zz + c13val*dux_dxl + c35val*(duz_dxl + dux_dzl) + c33val*duz_dzl
- sigma_xz = sigma_xz + c15val*dux_dxl + c55val*(duz_dxl + dux_dzl) + c35val*duz_dzl
+ if(assign_external_model) then
+ c11 = c11ext(ii2,jj2,ispec_elastic)
+ c13 = c13ext(ii2,jj2,ispec_elastic)
+ c15 = c15ext(ii2,jj2,ispec_elastic)
+ c33 = c33ext(ii2,jj2,ispec_elastic)
+ c35 = c35ext(ii2,jj2,ispec_elastic)
+ c55 = c55ext(ii2,jj2,ispec_elastic)
+ else
+ c11 = anisotropy(1,kmato(ispec_elastic))
+ c13 = anisotropy(2,kmato(ispec_elastic))
+ c15 = anisotropy(3,kmato(ispec_elastic))
+ c33 = anisotropy(4,kmato(ispec_elastic))
+ c35 = anisotropy(5,kmato(ispec_elastic))
+ c55 = anisotropy(6,kmato(ispec_elastic))
+ end if
+
+ sigma_xx = sigma_xx + c11*dux_dxl + c15*(duz_dxl + dux_dzl) + c13*duz_dzl
+ sigma_zz = sigma_zz + c13*dux_dxl + c35*(duz_dxl + dux_dzl) + c33*duz_dzl
+ sigma_xz = sigma_xz + c15*dux_dxl + c55*(duz_dxl + dux_dzl) + c35*duz_dzl
else
! no attenuation
sigma_xx = sigma_xx + lambdalplus2mul_relaxed*dux_dxl + lambdal_relaxed*duz_dzl
@@ -6103,11 +6081,26 @@
sigma_zz = lambdalplus2mul_relaxed*duz_dzl + lambdal_relaxed*dux_dxl
! full anisotropy
- if(TURN_ANISOTROPY_ON) then
+ if(anisotropic(ispec_elastic)) then
! implement anisotropy in 2D
- sigma_xx = c11val*dux_dxl + c15val*(duz_dxl + dux_dzl) + c13val*duz_dzl
- sigma_zz = c13val*dux_dxl + c35val*(duz_dxl + dux_dzl) + c33val*duz_dzl
- sigma_xz = c15val*dux_dxl + c55val*(duz_dxl + dux_dzl) + c35val*duz_dzl
+ if(assign_external_model) then
+ c11 = c11ext(i,j,ispec_elastic)
+ c13 = c13ext(i,j,ispec_elastic)
+ c15 = c15ext(i,j,ispec_elastic)
+ c33 = c33ext(i,j,ispec_elastic)
+ c35 = c35ext(i,j,ispec_elastic)
+ c55 = c55ext(i,j,ispec_elastic)
+ else
+ c11 = anisotropy(1,kmato(ispec_elastic))
+ c13 = anisotropy(2,kmato(ispec_elastic))
+ c15 = anisotropy(3,kmato(ispec_elastic))
+ c33 = anisotropy(4,kmato(ispec_elastic))
+ c35 = anisotropy(5,kmato(ispec_elastic))
+ c55 = anisotropy(6,kmato(ispec_elastic))
+ end if
+ sigma_xx = c11*dux_dxl + c15*(duz_dxl + dux_dzl) + c13*duz_dzl
+ sigma_zz = c13*dux_dxl + c35*(duz_dxl + dux_dzl) + c33*duz_dzl
+ sigma_xz = c15*dux_dxl + c55*(duz_dxl + dux_dzl) + c35*duz_dzl
endif
if(isolver == 2) then
@@ -6527,9 +6520,9 @@
xix,xiz,gammax,gammaz,jacobian,ibool,elastic,poroelastic,hprime_xx,hprime_zz, &
nspec,npoin,assign_external_model,it,deltat,t0(1),kmato,poroelastcoef,density, &
porosity,tortuosity, &
- vpext,vsext,rhoext,wxgll,wzgll,numat, &
+ vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,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,N_SLS)
+ potential_dot_acoustic,potential_dot_dot_acoustic,TURN_ATTENUATION_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 .or. it == NSTEP) then
@@ -6629,8 +6622,9 @@
call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,&
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
+ numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy,ispec,e1,e11, &
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
else if(.not. elastic(ispec) .and. .not. poroelastic(ispec)) then
@@ -6737,7 +6731,7 @@
valcurl = valcurl + dcurld*hlagrange
enddo
- enddo
+ enddo
! rotate seismogram components if needed, except if recording pressure, which is a scalar
if(seismotype /= 4 .and. seismotype /= 6) then
@@ -7235,8 +7229,9 @@
call compute_pressure_whole_medium(potential_dot_dot_acoustic,displ_elastic,&
displs_poroelastic,displw_poroelastic,elastic,poroelastic,vector_field_display, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
- numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext,e1,e11, &
- TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
+ numat,kmato,density,porosity,tortuosity,poroelastcoef,vpext,vsext,rhoext, &
+ c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,e1,e11, &
+ TURN_ATTENUATION_ON,Mu_nu1,Mu_nu2,N_SLS)
else if(imagetype == 4 .and. .not. p_sv) then
call exit_MPI('cannot draw pressure field for SH (membrane) waves')
@@ -7541,7 +7536,6 @@
'Read external initial field. . . . . .(initialfield) = ',l6/5x, &
'Add Bielak conditions . . . .(add_Bielak_conditions) = ',l6/5x, &
'Assign external model . . . .(assign_external_model) = ',l6/5x, &
- 'Turn anisotropy on or off. . . .(TURN_ANISOTROPY_ON) = ',l6/5x, &
'Turn attenuation on or off. . .(TURN_ATTENUATION_ON) = ',l6/5x, &
'Save grid in external file or not. . . .(outputgrid) = ',l6)
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2010-03-12 00:41:46 UTC (rev 16404)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2010-03-12 20:23:22 UTC (rev 16405)
@@ -130,27 +130,51 @@
open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown')
close(11,status='delete')
+! print*,'apres 1'
+! call flush(6)
+
open(unit=11,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown')
close(11,status='delete')
+! print*,'apres 2'
+! call flush(6)
+
open(unit=11,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown')
close(11,status='delete')
+! print*,'apres 3'
+! call flush(6)
+
open(unit=11,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown')
close(11,status='delete')
+ print*,'apres 4'
+ call flush(6)
+
open(unit=11,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown')
close(11,status='delete')
+! print*,'apres 5'
+! call flush(6)
+
open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown')
close(11,status='delete')
+! print*,'apres 6'
+! call flush(6)
+
open(unit=11,file='OUTPUT_FILES/Curl_file_single.bin',status='unknown')
close(11,status='delete')
+ ! print*,'apres 7'
+! call flush(6)
+
open(unit=11,file='OUTPUT_FILES/Curl_file_double.bin',status='unknown')
close(11,status='delete')
+! print*,'apres 8'
+! call flush(6)
+
endif
if ( myrank == 0 ) then
More information about the CIG-COMMITS
mailing list