[cig-commits] r8608 - in seismo/2D/SPECFEM2D/trunk: . DATA UTILS
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 16:00:48 PST 2007
Author: walter
Date: 2007-12-07 16:00:47 -0800 (Fri, 07 Dec 2007)
New Revision: 8608
Added:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30
seismo/2D/SPECFEM2D/trunk/DATA/interfaces_Ronan_SV_30.dat
seismo/2D/SPECFEM2D/trunk/UTILS/deriv_ricker_spatial.m
seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
added analytical plane wave and Bielak's conditions on the absorbing edges
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-08 00:00:47 UTC (rev 8608)
@@ -22,6 +22,7 @@
nx = 80 # number of elements along X
ngnod = 9 # number of control nodes per element (4 or 9)
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
@@ -100,4 +101,4 @@
1 80 1 20 1
1 80 21 40 2
1 80 41 60 3
-60 70 21 40 4
\ No newline at end of file
+60 70 21 40 4
Added: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Ronan_SV_30 2007-12-08 00:00:47 UTC (rev 8608)
@@ -0,0 +1,98 @@
+
+# title of job, and file that contains interface data
+title = Test for flat canyon Ronan SV 30 degrees
+interfacesfile = interfaces_Ronan_SV_30.dat
+
+# data concerning mesh, when generated using third-party app (more info in README)
+read_external_mesh = .false.
+mesh_file = ./DATA/unstructured_fluide_solide_test/mesh # file containing the mesh
+nodes_coords_file = ./DATA/unstructured_fluide_solide_test/nodes_coords # file containing the nodes coordinates
+materials_file = ./DATA/unstructured_fluide_solide_test/mat # file containing the material number for each element
+free_surface_file = ./DATA/unstructured_fluide_solide_test/surface_free # file containing the free surface
+absorbing_surface_file = ./DATA/unstructured_fluide_solide_test/surface_abs # file containing the absorbing surface
+
+# parameters concerning partitionning
+nproc = 1 # number of processes
+partionning_method = 1 # ascending order = 1, Metis = 2, Scotch = 3
+partitionning_strategy = 01110 #b{strat=m{asc=g,low=g,rat=1.0,type=h,vert=4}} #01110 # options concerning partitionning strategy.
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 19.d0 # abscissa of right side of the model
+nx = 70 # number of elements along X
+ngnod = 4 # number of control nodes per element (4 or 9)
+initialfield = .true. # use a plane wave as source or not
+add_Bielak_conditions = .true. # 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
+
+# absorbing boundaries parameters
+absorbing_conditions = .true. # absorbing boundary active or not
+absorbbottom = .true.
+absorbright = .true.
+absorbtop = .false.
+absorbleft = .true.
+
+# time step parameters
+nt = 20000 # total number of time steps
+deltat = 6.25e-4 # duration of a time step
+
+# source parameters
+source_surf = .true. # source inside the medium or at the surface
+xs = 0. # source location x in meters
+zs = 0. # source location z in meters
+source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
+time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0 = 8.0 # dominant source frequency (Hz) if not Dirac or Heaviside
+angleforce = 0. # angle of the source (for a force only)
+Mxx = 1. # Mxx component (for a moment tensor source only)
+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
+
+# constants for attenuation
+N_SLS = 2 # number of standard linear solids for attenuation
+Qp_attenuation = 136.4376068115 # quality factor P for attenuation
+Qs_attenuation = 136.4376068115 # quality factor S for attenuation
+f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver line parameters for seismograms
+seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS = .true. # creates a STATION file in ./DATA
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+
+# first receiver line
+nrec = 11 # number of receivers
+xdeb = 0. # first receiver x in meters
+zdeb = 0. # first receiver z in meters
+xfin = 10. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 0. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf = .true. # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO = 500 # display frequency in time steps
+output_postscript_snapshot = .true. # output Postscript snapshot of the results
+output_color_image = .true. # output color image of the results
+imagetype = 1 # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps = 1. # minimum amplitude in % for snapshots
+meshvect = .true. # display mesh on vector plots or not
+modelvect = .false. # display velocity model on vector plots
+boundvect = .true. # display boundary conditions on plots
+interpol = .true. # interpolation of the display or not
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp = 1 # subsampling of color snapshots
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
+gnuplot = .false. # generate a GNUPLOT file for the grid
+outputgrid = .false. # save the grid in a text file or not
+
+# velocity and density models
+nbmodels = 1 # nb of different models
+# define models as (model_number,1,rho,vp,vs,0,0) or (model_number,2,rho,c11,c13,c33,c44)
+# set vs to zero to make a given model acoustic
+# the mesh can contain both acoustic and elastic models simultaneously
+1 1 1.d0 2.d0 1.d0 0 0
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions = 1 # nb of regions and model number for each
+1 70 1 28 1
Added: seismo/2D/SPECFEM2D/trunk/DATA/interfaces_Ronan_SV_30.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/interfaces_Ronan_SV_30.dat 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/DATA/interfaces_Ronan_SV_30.dat 2007-12-08 00:00:47 UTC (rev 8608)
@@ -0,0 +1,26 @@
+#
+# number of interfaces
+#
+ 2
+#
+# for each interface below, we give the number of points and then x,y for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 5000 0
+#
+# interface number 2 (topography, top of the mesh)
+#
+ 2
+ 0 9
+ 5000 9
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 28
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-08 00:00:47 UTC (rev 8608)
@@ -46,24 +46,24 @@
# Portland
#F90 = /opt/openmpi-1.2.2/pgi64/bin/mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
-F90 = pgf90
-CC = pgcc
-FLAGS_NOCHECK=-fast -Mnobounds -Minline -Mneginfo -Mdclchk -Knoieee -Minform=warn -fastsse -tp amd64e -Msmart
-FLAGS_CHECK=-fast -Mbounds -Mneginfo -Mdclchk -Minform=warn
+#F90 = pgf90
+#CC = pgcc
+#FLAGS_NOCHECK=-fast -Mnobounds -Minline -Mneginfo -Mdclchk -Knoieee -Minform=warn -fastsse -tp amd64e -Msmart
+#FLAGS_CHECK=-fast -Mbounds -Mneginfo -Mdclchk -Minform=warn
# Intel
#F90 = ifort
#CC = gcc
-#FLAGS_NOCHECK=-O3 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -assume byterecl -check nobounds
-#FLAGS_CHECK = $(FLAGS_NOCHECK) -check bounds
+#FLAGS_NOCHECK=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -assume byterecl -check nobounds
+#FLAGS_CHECK = $(FLAGS_NOCHECK)
# GNU gfortran
#F90 = /opt/openmpi-1.2.1/gfortran64/bin/mpif90 -DUSE_MPI -DUSE_METIS -DUSE_SCOTCH
-#F90 = gfortran
-#CC = gcc
+F90 = gfortran
+CC = gcc
#FLAGS_NOCHECK = -O3 -march=opteron -m64 -mfpmath=sse,387
-#FLAGS_NOCHECK = -std=gnu -fimplicit-none -frange-check -O2 -Wunused-labels -Waliasing -Wampersand -Wsurprising -Wline-truncation -Wunderflow
-#FLAGS_CHECK = $(FLAGS_NOCHECK) -fbounds-check
+FLAGS_NOCHECK = -std=gnu -fimplicit-none -frange-check -O2 -Wunused-labels -Waliasing -Wampersand -Wsurprising -Wline-truncation -Wunderflow
+FLAGS_CHECK = $(FLAGS_NOCHECK) -fbounds-check
LINK = $(F90)
@@ -80,7 +80,7 @@
$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/attenuation_compute_param.o
+ $O/attenuation_compute_param.o $O/compute_Bielak_conditions.o
default: clean meshfem2D specfem2D convolve_source_timefunction
@@ -185,6 +185,9 @@
$O/compute_pressure.o: compute_pressure.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/compute_pressure.o compute_pressure.f90
+$O/compute_Bielak_conditions.o: compute_Bielak_conditions.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/compute_Bielak_conditions.o compute_Bielak_conditions.f90
+
$O/compute_arrays_source.o: compute_arrays_source.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/compute_arrays_source.o compute_arrays_source.f90
Added: seismo/2D/SPECFEM2D/trunk/UTILS/deriv_ricker_spatial.m
===================================================================
--- seismo/2D/SPECFEM2D/trunk/UTILS/deriv_ricker_spatial.m 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/UTILS/deriv_ricker_spatial.m 2007-12-08 00:00:47 UTC (rev 8608)
@@ -0,0 +1,50 @@
+
+(* Dimitri Komatitsch, University of Pau, Dec 2007 *)
+
+(* 'Mathematica' script to compute the spatial derivative of a Ricker in order
+to calculate the analytical traction for a plane wave for Bielak's technique
+on the absorbing edges *)
+
+ rac3 = Sqrt[3];
+ rac3sur2 = rac3 / 2;
+ HALF = 1 / 2;
+
+ rickertest[t_] := - (1 - 2 a t^2) Exp[-a t^2]
+
+ Ux[t_,x_,z_] := rac3sur2 rickertest[t - x/2 + (9 - z) rac3sur2] + rac3sur2 rickertest[t - x/2 - (9 - z) rac3sur2] + rac3 rickertest[t - x/2]
+
+ Uz[t_,x_,z_] := - HALF rickertest[t - x/2 + (9 - z) rac3sur2] + HALF rickertest[t - x/2 - (9 - z) rac3sur2]
+
+(* check displacements Ux and Uz *)
+(* Print[Simplify[Ux[t,x,z]]] *)
+
+(* Print[Simplify[Uz[t,x,z]]] *)
+
+(* output results to a file *)
+
+" " >> result.txt
+
+(* compute spatial derivatives of displacement *)
+
+" " >>> result.txt
+" dxUx = " >>> result.txt
+" " >>> result.txt
+ FortranForm[Simplify[D[Ux[t,x,z],x]]] >>> result.txt
+
+" " >>> result.txt
+" dzUx = " >>> result.txt
+" " >>> result.txt
+ FortranForm[Simplify[D[Ux[t,x,z],z]]] >>> result.txt
+
+" " >>> result.txt
+" dxUz = " >>> result.txt
+" " >>> result.txt
+ FortranForm[Simplify[D[Uz[t,x,z],x]]] >>> result.txt
+
+" " >>> result.txt
+" dzUz = " >>> result.txt
+" " >>> result.txt
+ FortranForm[Simplify[D[Uz[t,x,z],z]]] >>> result.txt
+
+" " >>> result.txt
+
Added: seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/compute_Bielak_conditions.f90 2007-12-08 00:00:47 UTC (rev 8608)
@@ -0,0 +1,164 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Copyright Universite de Pau et des Pays de l'Adour and CNRS, France.
+! 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
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic 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".
+!
+!========================================================================
+
+! compute analytical initial plane wave for Bielak's conditions
+
+ subroutine compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert)
+
+ implicit none
+
+ include "constants.h"
+
+ integer, intent(in) :: iglob,npoin,it
+
+ double precision, intent(in) :: deltat
+
+ double precision, intent(out) :: dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert
+
+ double precision, dimension(NDIM,npoin), intent(in) :: coord
+
+ double precision :: time_veloc,time_traction,t,x,z
+
+ double precision, external :: ricker_Bielak_veloc
+
+! get the coordinates of the mesh point
+ x = coord(1,iglob)
+ z = coord(2,iglob)
+
+! times for velocity and traction are staggered i.e. separated by deltat/2.d0
+!! DK DK YYYYYYYYYYYYYYYY UGLY il faudra peut-etre ajouter ou enlever deltat ici
+ time_veloc = (it-1)*deltat + deltat/2.d0 + time_offset
+ time_traction = time_veloc + deltat/2.d0
+
+ t = time_traction
+
+! derivatives of analytical expression of horizontal and vertical displacements,
+! computed using the "Mathematica" script in UTILS/deriv_ricker_spatial.m
+ dxUx = (sqrt(3.d0)*a*((-8*t + 4*x)*exp(-a*(t - x/2.d0)**2) + &
+ ((2*t - x)*(-2 + a*(-2*t + x)**2))*exp(-a*(t - x/2.d0)**2) + &
+ (2*(-2*t + x - sqrt(3.d0)*(-9 + z)))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ ((1 - (a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/2.d0)* &
+ (-2*t + x - sqrt(3.d0)*(-9 + z)))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ (2*(-2*t + x + sqrt(3.d0)*(-9 + z)))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ ((1 - (a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/2.d0)* &
+ (-2*t + x + sqrt(3.d0)*(-9 + z)))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0)))/4.d0
+
+ dzUx = (3*a*(((t + (-x + sqrt(3.d0)*(-9 + z))/2.d0)* &
+ (1 - (a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/2.d0))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0) - &
+ ((1 - (a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/2.d0)* &
+ (t - x/2.d0 - (sqrt(3.d0)*(-9 + z))/2.d0))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ (2*t - x + sqrt(3.d0)*(-9 + z))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ (-2*t + x + sqrt(3.d0)*(-9 + z))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0)))/2.d0
+
+ dxUz = (a*((2*t - x - sqrt(3.d0)*(-9 + z))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ (-2*t + x - sqrt(3.d0)*(-9 + z))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ ((1 - (a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/2.d0)* &
+ (-2*t + x - sqrt(3.d0)*(-9 + z)))/2.d0*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0) - &
+ ((1 - (a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/2.d0)* &
+ (-2*t + x + sqrt(3.d0)*(-9 + z)))/2.d0*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0)))/2.d0
+
+ dzUz = (sqrt(3.d0)*a*(((t + (-x + sqrt(3.d0)*(-9 + z))/2.d0)* &
+ (1 - (a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/2.d0))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ (2*t - x - sqrt(3.d0)*(-9 + z))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ ((1 - (a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/2.d0)* &
+ (t - x/2.d0 - (sqrt(3.d0)*(-9 + z))/2.d0))*exp(-(a*(-2*t + x + sqrt(3.d0)*(-9 + z))**2)/4.d0) + &
+ (2*t - x + sqrt(3.d0)*(-9 + z))*exp(-(a*(2*t - x + sqrt(3.d0)*(-9 + z))**2)/4.d0)))/2.d0
+
+ t = time_veloc
+
+! analytical expression of the two components of the velocity vector
+ veloc_horiz = rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2) &
+ + rac3 * ricker_Bielak_veloc(t - x/2.d0)
+ veloc_vert = - HALF * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + HALF * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2)
+
+ end subroutine compute_Bielak_conditions
+
+! ********
+
+! compute time variation of the source for analytical initial plane wave
+ double precision function ricker_Bielak_displ(t)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision :: t
+
+! Ricker
+ ricker_Bielak_displ = - (1 - 2*a*t**2)*exp(-a*t**2)
+
+ end function ricker_Bielak_displ
+
+! *******
+
+! compute time variation of the source for analytical initial plane wave
+ double precision function ricker_Bielak_veloc(t)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision :: t
+
+! first time derivative of a Ricker
+ ricker_Bielak_veloc = 2*a*t*(3 - 2*a*t**2)*exp(-a*t**2)
+
+ end function ricker_Bielak_veloc
+
+! *******
+
+! compute time variation of the source for analytical initial plane wave
+ double precision function ricker_Bielak_accel(t)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision :: t
+
+! second time derivative of a Ricker
+ ricker_Bielak_accel = 2*a*(3 - 12*a*t**2 + 4*a**2*t**4)* exp(-a*t**2)
+
+ end function ricker_Bielak_accel
+
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2007-12-08 00:00:47 UTC (rev 8608)
@@ -49,8 +49,7 @@
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, &
- nspec_inner_outer, ispec_inner_outer_to_glob, num_phase_inner_outer &
- )
+ nspec_inner_outer,ispec_inner_outer_to_glob,num_phase_inner_outer,deltat,coord,add_Bielak_conditions)
! compute forces for the elastic elements
@@ -60,7 +59,7 @@
integer :: npoin,nspec,nelemabs,numat,iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP
- logical :: anyabs,assign_external_model,initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON
+ logical :: anyabs,assign_external_model,initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,add_Bielak_conditions
double precision :: angleforce,deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
@@ -126,6 +125,11 @@
! 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
+ 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, external :: rickertest,rickertestvel
+
! only for the first call to compute_forces_elastic (during computation on outer elements)
if ( num_phase_inner_outer ) then
! compute Grad(displ_elastic) at time step n for attenuation
@@ -311,6 +315,19 @@
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) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert)
+ traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
+ traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ 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)
@@ -331,16 +348,16 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob)
- vz = veloc_elastic(2,iglob)
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ vz = veloc_elastic(2,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - tx*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - tz*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz + traction_z_t0)*weight
endif
enddo
@@ -356,6 +373,19 @@
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) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert)
+ traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
+ traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ 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)
@@ -376,16 +406,16 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob)
- vz = veloc_elastic(2,iglob)
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ vz = veloc_elastic(2,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - tx*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - tz*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz - traction_z_t0)*weight
endif
enddo
@@ -407,6 +437,19 @@
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) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert)
+ 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)
@@ -427,16 +470,16 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob)
- vz = veloc_elastic(2,iglob)
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ vz = veloc_elastic(2,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - tx*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - tz*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx + traction_x_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz + traction_z_t0)*weight
endif
enddo
@@ -458,6 +501,19 @@
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) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert)
+ 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)
@@ -478,16 +534,16 @@
! Clayton-Engquist condition if elastic
if(elastic(ispec)) then
- vx = veloc_elastic(1,iglob)
- vz = veloc_elastic(2,iglob)
+ vx = veloc_elastic(1,iglob) - veloc_horiz
+ vz = veloc_elastic(2,iglob) - veloc_vert
vn = nx*vx+nz*vz
tx = rho_vp*vn*nx+rho_vs*(vx-vn*nx)
tz = rho_vp*vn*nz+rho_vs*(vz-vn*nz)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) - tx*weight
- accel_elastic(2,iglob) = accel_elastic(2,iglob) - tz*weight
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) - (tx - traction_x_t0)*weight
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - (tz - traction_z_t0)*weight
endif
enddo
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-08 00:00:47 UTC (rev 8608)
@@ -1,3 +1,5 @@
+
+! size of real and double precision numbers in bytes
integer, parameter :: SIZE_REAL = 4
integer, parameter :: SIZE_DOUBLE = 8
@@ -81,6 +83,13 @@
! was found by trial and error
double precision, parameter :: SOURCE_DECAY_MIMIC_TRIANGLE = 1.628d0
+! for the analytical initial plane wave source for Bielak conditions
+ double precision, parameter :: time_offset = 3.575d0
+ double precision, parameter :: f0_ricker_Bielak = 1.d0
+ double precision, parameter :: a = PI*PI*f0_ricker_Bielak*f0_ricker_Bielak
+ double precision, parameter :: rac3 = 1.7320508075688772935d0
+ double precision, parameter :: rac3sur2 = rac3 / 2.d0
+
! non linear display to enhance small amplitudes in color images
double precision, parameter :: POWER_DISPLAY_COLOR = 0.30d0
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-08 00:00:47 UTC (rev 8608)
@@ -123,7 +123,7 @@
logical interpol,gnuplot,assign_external_model,outputgrid
logical abstop,absbottom,absleft,absright,any_abs
- logical source_surf,meshvect,initialfield,modelvect,boundvect
+ logical source_surf,meshvect,initialfield,modelvect,boundvect,add_Bielak_conditions
logical TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
logical, dimension(:), allocatable :: enreg_surf
@@ -285,6 +285,7 @@
endif
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)
@@ -1112,8 +1113,8 @@
write(15,*) 'anglerec'
write(15,*) anglerec
- write(15,*) 'initialfield'
- write(15,*) initialfield
+ write(15,*) 'initialfield add_Bielak_conditions'
+ write(15,*) initialfield,add_Bielak_conditions
write(15,*) 'seismotype imagetype'
write(15,*) seismotype,imagetype
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-03 15:01:18 UTC (rev 8607)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2007-12-08 00:00:47 UTC (rev 8608)
@@ -198,7 +198,7 @@
ibegin_bottom,iend_bottom,ibegin_top,iend_top,jbegin_left,jend_left,jbegin_right,jend_right
integer ispec_selected_source,iglob_source,ix_source,iz_source,is_proc_source,nb_proc_source
- double precision a,displnorm_all,displnorm_all_glob
+ double precision aval,displnorm_all,displnorm_all_glob
real(kind=CUSTOM_REAL), dimension(:), allocatable :: source_time_function
double precision, external :: netlib_specfun_erf
@@ -209,7 +209,7 @@
logical interpol,meshvect,modelvect,boundvect,assign_external_model,initialfield, &
outputgrid,gnuplot,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,output_postscript_snapshot,output_color_image, &
- plot_lowerleft_corner_only
+ plot_lowerleft_corner_only,add_Bielak_conditions
double precision :: cutsnaps,sizemax_arrows,anglerec,xirec,gammarec
@@ -342,6 +342,10 @@
integer, dimension(:,:,:), allocatable :: copy_ibool_ori
integer :: inumber
+! to compute analytical initial plane wave field
+ double precision :: t
+ double precision, external :: ricker_Bielak_displ,ricker_Bielak_veloc,ricker_Bielak_accel
+
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
@@ -423,7 +427,8 @@
read(IIN,*) anglerec
read(IIN,"(a80)") datlin
- read(IIN,*) initialfield
+ read(IIN,*) initialfield,add_Bielak_conditions
+ if(add_Bielak_conditions .and. .not. initialfield) stop 'need to have an initial field to add Bielak plane wave conditions'
read(IIN,"(a80)") datlin
read(IIN,*) seismotype,imagetype
@@ -437,7 +442,7 @@
write(IOUT,200) npgeo,NDIM
write(IOUT,600) NTSTEP_BETWEEN_OUTPUT_INFO,colors,numbers
write(IOUT,700) seismotype,anglerec
- write(IOUT,750) initialfield,assign_external_model,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
+ write(IOUT,750) initialfield,add_Bielak_conditions,assign_external_model,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
write(IOUT,800) imagetype,100.d0*cutsnaps,subsamp
!---- read time step
@@ -479,7 +484,7 @@
endif
! for the source time function
- a = pi*pi*f0*f0
+ aval = pi*pi*f0*f0
!----- convert angle from degrees to radians
angleforce = angleforce * pi / 180.d0
@@ -1442,24 +1447,65 @@
write(IOUT,*) 'Reading initial fields from external file...'
write(IOUT,*)
if(any_acoustic) call exit_MPI('initial field currently implemented for purely elastic simulation only')
- open(unit=55,file='OUTPUT_FILES/wavefields.txt',status='unknown')
- read(55,*) nbpoin
- if(nbpoin /= npoin) call exit_MPI('Wrong number of points in input file')
- allocate(displread(NDIM))
- allocate(velocread(NDIM))
- allocate(accelread(NDIM))
- do n = 1,npoin
- read(55,*) inump, (displread(i), i=1,NDIM), &
- (velocread(i), i=1,NDIM), (accelread(i), i=1,NDIM)
- if(inump<1 .or. inump>npoin) call exit_MPI('Wrong point number')
- displ_elastic(:,inump) = displread
- veloc_elastic(:,inump) = velocread
- accel_elastic(:,inump) = accelread
- enddo
- deallocate(displread)
- deallocate(velocread)
- deallocate(accelread)
- close(55)
+
+ if(.not. add_Bielak_conditions) then
+
+ open(unit=55,file='OUTPUT_FILES/wavefields.txt',status='unknown')
+ read(55,*) nbpoin
+ if(nbpoin /= npoin) call exit_MPI('Wrong number of points in input file')
+ allocate(displread(NDIM))
+ allocate(velocread(NDIM))
+ allocate(accelread(NDIM))
+ do n = 1,npoin
+ read(55,*) inump, (displread(i), i=1,NDIM), (velocread(i), i=1,NDIM), (accelread(i), i=1,NDIM)
+ if(inump<1 .or. inump>npoin) call exit_MPI('Wrong point number')
+ displ_elastic(:,inump) = displread
+ veloc_elastic(:,inump) = velocread
+ accel_elastic(:,inump) = accelread
+ enddo
+ deallocate(displread)
+ deallocate(velocread)
+ deallocate(accelread)
+ close(55)
+
+ else
+
+! compute analytical initial plane wave field
+ print *,'computing analytical initial plane wave field for SV wave at 30 degrees and Poisson = 0.25'
+
+ do i = 1,npoin
+
+ x = coord(1,i)
+ z = coord(2,i)
+
+! add a time offset in order for the initial field to be inside the medium
+ t = 0.d0 + time_offset
+
+! initial analytical displacement
+ displ_elastic(1,i) = rac3sur2 * ricker_Bielak_displ(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + rac3sur2 * ricker_Bielak_displ(t - x/2.d0 - (9 - z) * rac3sur2) &
+ + rac3 * ricker_Bielak_displ(t - x/2.d0)
+ displ_elastic(2,i) = - HALF * ricker_Bielak_displ(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + HALF * ricker_Bielak_displ(t - x/2.d0 - (9 - z) * rac3sur2)
+
+! initial analytical velocity
+ veloc_elastic(1,i) = rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + rac3sur2 * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2) &
+ + rac3 * ricker_Bielak_veloc(t - x/2.d0)
+ veloc_elastic(2,i) = - HALF * ricker_Bielak_veloc(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + HALF * ricker_Bielak_veloc(t - x/2.d0 - (9 - z) * rac3sur2)
+
+! initial analytical acceleration
+ accel_elastic(1,i) = rac3sur2 * ricker_Bielak_accel(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + rac3sur2 * ricker_Bielak_accel(t - x/2.d0 - (9 - z) * rac3sur2) &
+ + rac3 * ricker_Bielak_accel(t - x/2.d0)
+ accel_elastic(2,i) = - HALF * ricker_Bielak_accel(t - x/2.d0 + (9 - z) * rac3sur2) &
+ + HALF * ricker_Bielak_accel(t - x/2.d0 - (9 - z) * rac3sur2)
+
+ enddo
+
+ endif
+
write(IOUT,*) 'Max norm of initial elastic displacement = ',maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(2,:)**2))
endif
@@ -1490,15 +1536,15 @@
! Ricker (second derivative of a Gaussian) source time function
if(time_function_type == 1) then
- source_time_function(it) = - factor * (ONE-TWO*a*(time-t0)**2) * exp(-a*(time-t0)**2)
+ source_time_function(it) = - factor * (ONE-TWO*aval*(time-t0)**2) * exp(-aval*(time-t0)**2)
! first derivative of a Gaussian source time function
else if(time_function_type == 2) then
- source_time_function(it) = - factor * TWO*a*(time-t0) * exp(-a*(time-t0)**2)
+ source_time_function(it) = - factor * TWO*aval*(time-t0) * exp(-aval*(time-t0)**2)
! Gaussian or Dirac (we use a very thin Gaussian instead) source time function
else if(time_function_type == 3 .or. time_function_type == 4) then
- source_time_function(it) = factor * exp(-a*(time-t0)**2)
+ source_time_function(it) = factor * exp(-aval*(time-t0)**2)
! Heaviside source time function (we use a very thin error function instead)
else if(time_function_type == 5) then
@@ -2009,8 +2055,7 @@
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, &
- nspec_outer, ispec_outer_to_glob, .true. &
- )
+ nspec_outer, ispec_outer_to_glob,.true.,deltat,coord,add_Bielak_conditions)
! *********************************************************
! ************* add coupling with the acoustic side
@@ -2108,8 +2153,7 @@
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, &
- nspec_inner, ispec_inner_to_glob, .false. &
- )
+ nspec_inner, ispec_inner_to_glob,.false.,deltat,coord,add_Bielak_conditions)
! assembling accel_elastic for elastic elements (receive)
#ifdef USE_MPI
@@ -2567,6 +2611,7 @@
750 format(//1x,'C o n t r o l',/1x,13('='),//5x, &
'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, &
More information about the cig-commits
mailing list