[cig-commits] r18803 - in seismo/2D/SPECFEM2D/trunk: DATA src/specfem2D
rmodrak at geodynamics.org
rmodrak at geodynamics.org
Wed Jul 27 17:19:51 PDT 2011
Author: rmodrak
Date: 2011-07-27 17:19:50 -0700 (Wed, 27 Jul 2011)
New Revision: 18803
Added:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_1
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_2
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_3
seismo/2D/SPECFEM2D/trunk/DATA/SOURCE_noise
seismo/2D/SPECFEM2D/trunk/DATA/STATIONS_target_noise
seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
adds noise tomography routines
Added: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_1
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_1 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_1 2011-07-28 00:19:50 UTC (rev 18803)
@@ -0,0 +1,115 @@
+# title of job
+title = Noise_2D
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
+
+# parameters concerning partitioning
+nproc = 1 # number of processes
+partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
+
+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
+READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+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
+freq0 = 10 # frequency for viscous attenuation
+p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt = 6000 # total number of time steps
+deltat = 1.d-3 # duration of a time step
+
+# source parameters
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS = 2 # number of standard linear solids for attenuation
+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 = .false. # creates a STATION file in ./DATA
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
+
+# first receiver line (repeat these 6 lines and adjust nreceiverlines accordingly)
+nrec = 1 # number of receivers
+xdeb = 3000. # first receiver x in meters
+zdeb = 1500. # first receiver z in meters
+xfin = 3000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1500. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO = 100 # 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
+output_grid = .false. # save the grid in a text file or not
+output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels = 1 # nb of different models
+# define models as
+# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
+# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) 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, and poroelastic models simultaneously.
+1 1 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh = .false.
+
+# absorbing boundary active or not
+absorbing_conditions = .true.
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file = ./DATA/Mesh_canyon/canyon_mesh_file # file containing the mesh
+nodes_coords_file = ./DATA/Mesh_canyon/canyon_nodes_coords_file # file containing the nodes coordinates
+materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
+free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
+absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile = ../EXAMPLES/noise_example/uniform.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 4000.d0 # abscissa of right side of the model
+nx = 80 # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom = .true.
+absorbright = .true.
+absorbtop = .true.
+absorbleft = .true.
+
+# 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 80 1 60 1
Added: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_2
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_2 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_2 2011-07-28 00:19:50 UTC (rev 18803)
@@ -0,0 +1,115 @@
+# title of job
+title = Noise_2D
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+SAVE_FORWARD = .true. # save the last frame, needed for adjoint simulation
+
+# parameters concerning partitioning
+nproc = 1 # number of processes
+partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
+
+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
+READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+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
+freq0 = 10 # frequency for viscous attenuation
+p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt = 6000 # total number of time steps
+deltat = 1.d-3 # duration of a time step
+
+# source parameters
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS = 2 # number of standard linear solids for attenuation
+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 = .false. # creates a STATION file in ./DATA
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
+
+# first receiver line (repeat these 6 lines and adjust nreceiverlines accordingly)
+nrec = 1 # number of receivers
+xdeb = 3000. # first receiver x in meters
+zdeb = 1500. # first receiver z in meters
+xfin = 3000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1500. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO = 100 # 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
+output_grid = .false. # save the grid in a text file or not
+output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels = 1 # nb of different models
+# define models as
+# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
+# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) 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, and poroelastic models simultaneously.
+1 1 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh = .false.
+
+# absorbing boundary active or not
+absorbing_conditions = .true.
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file = ./DATA/Mesh_canyon/canyon_mesh_file # file containing the mesh
+nodes_coords_file = ./DATA/Mesh_canyon/canyon_nodes_coords_file # file containing the nodes coordinates
+materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
+free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
+absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile = ../EXAMPLES/noise_example/uniform.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 4000.d0 # abscissa of right side of the model
+nx = 80 # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom = .true.
+absorbright = .true.
+absorbtop = .true.
+absorbleft = .true.
+
+# 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 80 1 60 1
Added: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_3
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_3 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_noise_3 2011-07-28 00:19:50 UTC (rev 18803)
@@ -0,0 +1,115 @@
+# title of job
+title = Noise_2D
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 2 # 1 = forward, 2 = adjoint + kernels
+SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
+
+# parameters concerning partitioning
+nproc = 1 # number of processes
+partitioning_method = 3 # SCOTCH = 3, ascending order (very bad idea) = 1
+
+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
+READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+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
+freq0 = 10 # frequency for viscous attenuation
+p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt = 6000 # total number of time steps
+deltat = 1.d-3 # duration of a time step
+
+# source parameters
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS = 2 # number of standard linear solids for attenuation
+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 = .false. # creates a STATION file in ./DATA
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+rec_normal_to_surface = .false. # base anglerec normal to surface (external mesh and curve file needed)
+
+# first receiver line (repeat these 6 lines and adjust nreceiverlines accordingly)
+nrec = 1 # number of receivers
+xdeb = 3000. # first receiver x in meters
+zdeb = 1500. # first receiver z in meters
+xfin = 3000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1500. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO = 100 # 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
+output_grid = .false. # save the grid in a text file or not
+output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot = .false. # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels = 1 # nb of different models
+# define models as
+# I: (model_number 1 rho Vp Vs 0 0 Qp Qs 0 0 0 0 0 0) or
+# II: (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) 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, and poroelastic models simultaneously.
+1 1 2300.d0 2450.d0 1425.d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh = .false.
+
+# absorbing boundary active or not
+absorbing_conditions = .true.
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file = ./DATA/Mesh_canyon/canyon_mesh_file # file containing the mesh
+nodes_coords_file = ./DATA/Mesh_canyon/canyon_nodes_coords_file # file containing the nodes coordinates
+materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
+free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
+absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile = ../EXAMPLES/noise_examples/uniform.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 4000.d0 # abscissa of right side of the model
+nx = 80 # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom = .true.
+absorbright = .true.
+absorbtop = .true.
+absorbleft = .true.
+
+# 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 80 1 60 1
Added: seismo/2D/SPECFEM2D/trunk/DATA/SOURCE_noise
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/SOURCE_noise (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/DATA/SOURCE_noise 2011-07-28 00:19:50 UTC (rev 18803)
@@ -0,0 +1,13 @@
+#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 = 0. # source location z in meters
+source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
+time_function_type = 3 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0 = 10.0 # dominant source frequency (Hz) if not Dirac or Heaviside
+tshift = 0.0 # time shift when multi sources (if one source, must be zero)
+angleforce = 0. # angle of the source (for a force only)
+Mxx = 0.d0 # Mxx component (for a moment tensor source only)
+Mzz = 0.d0 # Mzz component (for a moment tensor source only)
+Mxz = 0.d0 # Mxz component (for a moment tensor source only)
+factor = 0.d0 # amplification factor
Added: seismo/2D/SPECFEM2D/trunk/DATA/STATIONS_target_noise
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/STATIONS_target_noise (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/DATA/STATIONS_target_noise 2011-07-28 00:19:50 UTC (rev 18803)
@@ -0,0 +1,3 @@
+S0001 AA 1000.0000000 1500.0000000 0.0 0.0
+S0002 AA 2000.0000000 1500.0000000 0.0 0.0
+S0003 AA 3000.0000000 1500.0000000 0.0 0.0
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in 2011-07-24 03:30:44 UTC (rev 18802)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/Makefile.in 2011-07-28 00:19:50 UTC (rev 18803)
@@ -140,6 +140,7 @@
$O/locate_source_force.o \
$O/locate_source_moment_tensor.o \
$O/netlib_specfun_erf.o \
+ $O/noise_tomography.o \
$O/paco_beyond_critical.o \
$O/paco_convolve_fft.o \
$O/plotgll.o \
@@ -332,6 +333,9 @@
$O/netlib_specfun_erf.o: ${S}/netlib_specfun_erf.f90
${F90} $(FLAGS_CHECK) -c -o $O/netlib_specfun_erf.o ${S}/netlib_specfun_erf.f90
+$O/noise_tomography.o: ${S}/noise_tomography.f90 ${SETUP}/constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/noise_tomography.o ${S}/noise_tomography.f90
+
$O/paco_beyond_critical.o: ${S}/paco_beyond_critical.f90 ${SETUP}/constants.h
${F90} $(FLAGS_CHECK) -c -o $O/paco_beyond_critical.o ${S}/paco_beyond_critical.f90
Added: seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-07-28 00:19:50 UTC (rev 18803)
@@ -0,0 +1,425 @@
+!========================================================================
+!
+! S P E C F E M 2 D Version 6 . 2
+! ------------------------------
+!
+! Copyright Universite de Pau, CNRS and INRIA, France,
+! and Princeton University / California Institute of Technology, 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
+!
+! 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".
+!
+!========================================================================
+
+!NOISE TOMOGRAPHY TO DO LIST
+
+!1. Replace ad hoc scaling in add_surface_movie_noise by expression involving Jacobian
+!2. Move NOISE_TOMOGRAPHY flag from source code to Par_file
+!3. Use separate STATIONS_ADJOINT file
+!4. Add test case (uniform model) under EXAMPLES
+!5. Update manual
+
+! =============================================================================================================
+! specify spatial distribution of microseismic noise sources
+! USERS need to modify this subroutine to suit their own needs
+ subroutine create_mask_noise(p_sv,nglob,coord,angle_noise,mask_noise)
+
+ implicit none
+ include "constants.h"
+
+ !input
+ logical :: p_sv
+ integer :: nglob
+ real(kind=CUSTOM_REAL) :: angle_noise
+ real(kind=CUSTOM_REAL), dimension(2,nglob) :: coord
+
+ !output
+ real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
+
+ !local
+ integer :: iglob, ios
+ real(kind=CUSTOM_REAL) :: xx,zz,mask_noise_out
+
+ !specify distribution of noise sources as a function of xx, zz
+ do iglob = 1,nglob
+ xx = coord(1,iglob)
+ zz = coord(2,iglob)
+
+ !below, the noise is assumed to be uniform; users are free to
+ !to change this expression to one involving xx, zz
+ mask_noise(iglob) = 1.0
+
+ enddo
+
+ end subroutine create_mask_noise
+
+
+! =============================================================================================================
+! read noise parameters and check for consitency
+ subroutine read_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
+ any_acoustic,any_poroelastic,p_sv,irec_master, &
+ Mxx,Mxz,Mzz,factor,NSOURCES, &
+ xi_receiver,gamma_receiver,ispec_selected_rec,nrec, &
+ xi_noise,gamma_noise,ispec_noise,angle_noise)
+
+ implicit none
+ include "constants.h"
+
+ !input
+ integer :: NOISE_TOMOGRAPHY, SIMULATION_TYPE
+ logical :: SAVE_FORWARD
+ logical :: any_acoustic, any_poroelastic
+ logical :: p_sv
+ integer :: NSOURCES, nrec
+ double precision, dimension(NSOURCES) :: Mxx, Mxz, Mzz, factor
+ double precision, dimension(nrec) :: xi_receiver, gamma_receiver
+ integer, dimension(nrec) :: ispec_selected_rec
+ double precision :: xi_noise, gamma_noise, angle_noise
+ integer :: ispec_noise
+
+
+ !output
+ integer :: irec_master
+
+ !local
+ integer :: i,ios
+
+ !define master receiver
+ open(unit=509,file='NOISE_TOMOGRAPHY/irec_master',status='old',action='read',iostat=ios)
+ if (ios == 0) then
+ read(509,*) irec_master
+ else
+ irec_master=1
+ write(*,*) ''
+ write(*,*) 'Error opening NOISE_TOMOGRAPHY/irec_master.'
+ write(*,*) 'Using irec_master=1. Continuing.'
+ write(*,*) ''
+ endif
+ close(509)
+
+ if ( (NOISE_TOMOGRAPHY == 1) .and. (irec_master > nrec .or. irec_master < 1) ) &
+ call exit_mpi('irec_master out of range of given number of receivers. Exiting.')
+
+ xi_noise = xi_receiver(irec_master)
+ gamma_noise = gamma_receiver(irec_master)
+ ispec_noise = ispec_selected_rec(irec_master)
+ angle_noise = 0.d0
+
+
+ !check simulation parameters
+ if (NOISE_TOMOGRAPHY==1) then
+ if (SIMULATION_TYPE/=1) call exit_mpi('NOISE_TOMOGRAPHY=1 requires SIMULATION_TYPE=1 -> check DATA/Par_file')
+
+
+
+ else if (NOISE_TOMOGRAPHY==2) then
+ if (SIMULATION_TYPE/=1) call exit_mpi('NOISE_TOMOGRAPHY=2 requires SIMULATION_TYPE=1 -> check DATA/Par_file')
+
+ if (.not. SAVE_FORWARD) call exit_mpi('NOISE_TOMOGRAPHY=2 requires SAVE_FORWARD=.true. -> check DATA/Par_file')
+
+
+
+ else if (NOISE_TOMOGRAPHY==3) then
+ if (SIMULATION_TYPE/=2) call exit_mpi('NOISE_TOMOGRAPHY=3 requires SIMULATION_TYPE=2 -> check DATA/Par_file')
+
+ if (SAVE_FORWARD) call exit_mpi('NOISE_TOMOGRAPHY=3 requires SAVE_FORWARD=.false. -> check DATA/Par_file')
+
+ endif
+
+
+! check model parameters
+ if (any_acoustic) call exit_mpi('Acoustic models not yet implemented for noise simulations. Exiting.')
+ if (any_poroelastic) call exit_mpi('Poroelastic models not yet implemented for noise simulations. Exiting.')
+
+
+! moment tensor elements must be zero!
+ do i=1,NSOURCES
+ if ( (Mxx(i) /= 0.d0) .or. (Mxz(i) /= 0.d0) .or. (Mzz(i) /= 0.d0) .or. &
+ (factor(i) /= 0.d0)) then
+ call exit_mpi('For noise simulations, all moment tensor elements must be zero. Exiting.')
+ endif
+ enddo
+
+
+ end subroutine read_parameters_noise
+
+
+! =============================================================================================================
+! read in time series based on noise spectrum and construct noise "source" array
+ subroutine compute_source_array_noise(p_sv,NSTEP,deltat,nglob,ibool,ispec_noise, &
+ xi_noise,gamma_noise,angle_noise,xigll,zigll, &
+ time_function_noise,source_array_noise)
+
+ implicit none
+ include "constants.h"
+
+ !input parameters
+ logical :: p_sv
+ integer NSTEP, ispec_noise,nglob
+ integer, dimension(NGLLX,NGLLZ,nglob) :: ibool
+ real(kind=CUSTOM_REAL) :: deltat, xi_noise, gamma_noise, angle_noise
+
+ !output paramters
+ real(kind=CUSTOM_REAL), dimension(NSTEP) :: time_function_noise
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLZ,NSTEP) :: source_array_noise
+
+ !local parameters
+ integer, parameter :: METHOD = 1
+ integer :: it
+ real(kind=CUSTOM_REAL) :: t
+
+ integer :: i,j,iglob
+ real(kind=CUSTOM_REAL) :: factor_noise
+ real(kind=CUSTOM_REAL) :: xigll, zigll
+ real(kind=CUSTOM_REAL), dimension(NGLLX) :: hxi, hpxi
+ real(kind=CUSTOM_REAL), dimension(NGLLZ) :: hgamma, hpgamma
+
+
+ !local paramters, for method 1
+ real(kind=CUSTOM_REAL) ::f0_noise, aval_noise, t0_noise
+
+ !local paramters, for method 2
+ character(len=60) :: file_in_noise
+
+
+ if(METHOD == 1) then
+
+ ! METHOD 1: use Ricker function
+ time_function_noise(:) = 0._CUSTOM_REAL
+ f0_noise = 15.0d0
+ aval_noise = PI*PI*f0_noise*f0_noise
+ t0_noise = ((NSTEP-1)/2.)*deltat + 0.05d0
+ factor_noise = 1.d3
+
+ do it = 1,NSTEP
+ t = it*deltat
+ time_function_noise(it) = factor_noise * &
+ (1. - 2.*aval_noise*((t-t0_noise)**2.)) * &
+ exp(-aval_noise*(t-t0_noise)**2.)
+ enddo
+
+ elseif(METHOD == 2) then
+
+ ! METHOD 2: read in custom noise function
+ open(unit=55,file='NOISE_TOMOGRAPHY/S_squared',status='old')
+ do it = 1,NSTEP
+ READ(55,*) time_function_noise(it)
+ enddo
+ close(55)
+
+ endif
+
+ !interpolate over GLL points
+ source_array_noise(:,:,:,:) = 0._CUSTOM_REAL
+ call lagrange_any(xi_noise,NGLLX,xigll,hxi,hpxi)
+ call lagrange_any(gamma_noise,NGLLZ,zigll,hgamma,hpgamma)
+
+ if(p_sv) then ! P-SV simulation
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_noise)
+ source_array_noise(1,i,j,:) = time_function_noise(:) * hxi(i) * hgamma(j)
+ source_array_noise(3,i,j,:) = time_function_noise(:) * hxi(i) * hgamma(j)
+ enddo
+ enddo
+ else ! SH (membrane) simulation
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_noise)
+ source_array_noise(2,i,j,:) = time_function_noise(:) * hxi(i) * hgamma(j)
+ enddo
+ enddo
+ endif
+
+ end subroutine compute_source_array_noise
+
+
+! =============================================================================================================
+! inject the "source" that drives the "generating wavefield"
+ subroutine add_point_source_noise(p_sv,it,NSTEP,nglob,ibool,ispec_noise, &
+ accel_elastic,angle_noise,time_function_noise, &
+ source_array_noise)
+ implicit none
+ include "constants.h"
+
+ !input parameters
+ logical :: p_sv
+ integer :: it, NSTEP
+ integer :: ispec_noise, nglob
+ integer, dimension(NGLLX,NGLLZ,nglob) :: ibool
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: accel_elastic
+ real(kind=CUSTOM_REAL) :: angle_noise
+ real(kind=CUSTOM_REAL), dimension(NSTEP) :: time_function_noise
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLZ,NSTEP) :: source_array_noise
+
+ !local variables
+ integer :: i,j,iglob
+
+ if(p_sv) then ! P-SV calculation
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_noise)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) &
+ + sin(angle_noise)*source_array_noise(1,i,j,it)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) &
+ - cos(angle_noise)*source_array_noise(3,i,j,it)
+ enddo
+ enddo
+ else ! SH (membrane) calculation
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec_noise)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) - source_array_noise(2,i,j,it)
+ enddo
+ enddo
+ endif
+
+ end subroutine add_point_source_noise
+
+
+! =============================================================================================================
+! read in and inject the "source" that drives the "enemble forward wavefield"
+! (recall that the ensemble forward wavefield has a spatially distributed source)
+ subroutine add_surface_movie_noise(p_sv,it,NSTEP,nglob,ibool, &
+ accel_elastic,surface_movie_x_noise,surface_movie_y_noise, &
+ surface_movie_z_noise,mask_noise)
+
+ implicit none
+ include "constants.h"
+
+ !input parameters
+ logical :: p_sv
+ integer :: it, NSTEP
+ integer :: nglob
+ integer, dimension(NGLLX,NGLLZ,nglob) :: ibool
+
+ real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
+ real(kind=CUSTOM_REAL), dimension(nglob) :: &
+ surface_movie_x_noise, surface_movie_y_noise, surface_movie_z_noise
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: accel_elastic
+
+ !local variables
+ integer :: i,j,iglob,ios
+ real(kind=CUSTOM_REAL) :: factor_noise, hlagrange
+ character(len=60) :: file_in_noise
+
+
+ write(file_in_noise,"('movie_noise_',i6.6)") NSTEP-it+1
+ open(unit=500,file='OUTPUT_FILES/'//trim(file_in_noise),status='old',form='unformatted',action='read',iostat=ios)
+ if( ios /= 0) write(*,*) 'Error retrieving generating wavefield.'
+ if(p_sv) then
+ read(500) surface_movie_x_noise
+ read(500) surface_movie_z_noise
+ else
+ read(500) surface_movie_y_noise
+ endif
+ close(500)
+
+ factor_noise = 1.d6
+
+ if(p_sv) then ! P-SV calculation
+ accel_elastic(1,:) = accel_elastic(1,:) + factor_noise * mask_noise(:) * surface_movie_x_noise(:)
+ accel_elastic(3,:) = accel_elastic(3,:) + factor_noise * mask_noise(:) * surface_movie_z_noise(:)
+! note: above implementation does not yet include non-vertical noise
+! excitation
+
+ else ! SH (membrane) calculation
+ accel_elastic(2,:) = accel_elastic(2,:) + factor_noise * mask_noise(:) * surface_movie_y_noise(:)
+
+! note: above implementation includes ad hoc scaling by a number
+! "factor_noise", meant to approximate the correct jacobian and gll weighting
+! correct implementation should look something like:
+! do iglob = 1,nglob
+! hlagrange = weight(iglob) * jacobian(iglob)
+! accel_elastic(2,iglob) = accel_elastic(2,iglob) &
+! + hlagrange * surface_movie_y_noise(iglob)
+! enddo
+! enddo
+! enddo
+
+ endif
+
+ end subroutine add_surface_movie_noise
+
+! =============================================================================================================
+! save a snapshot of the "generating wavefield" that will be used to drive
+! the "ensemble forward wavefield"
+ subroutine save_surface_movie_noise(p_sv,it,nglob,displ_elastic)
+
+ implicit none
+ include "constants.h"
+
+ !input paramters
+ logical :: p_sv
+ integer :: it, nglob
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: displ_elastic
+
+ !local parameters
+ character(len=60) file_out_noise
+
+ write(file_out_noise,"('movie_noise_',i6.6)") it
+ open(unit=500,file='OUTPUT_FILES/'//trim(file_out_noise),status='unknown',form='unformatted',action='write')
+
+ if(p_sv) then
+ write(500) displ_elastic(1,:)
+ write(500) displ_elastic(3,:)
+ else
+ write(500) displ_elastic(2,:)
+ endif
+
+ end subroutine save_surface_movie_noise
+
+! =============================================================================================================
+! auxillary routine for writing out the array "mask_noise"; uses a format similar to
+! the one currently used for writing kernels
+ subroutine write_mask_noise(nglob,coord,mask_noise)
+
+ implicit none
+ include "constants.h"
+
+ !input paramters
+ integer :: nglob
+ real(kind=CUSTOM_REAL), dimension(2,nglob) :: coord
+ real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
+
+ !local parameters
+ integer :: iglob
+
+ open(unit=504,file='OUTPUT_FILES/mask_noise',status='unknown',action='write')
+ do iglob = 1, nglob
+ write(504,*) coord(1,iglob), coord(2,iglob), mask_noise(iglob)
+ enddo
+ close(504)
+
+
+ end subroutine write_mask_noise
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-07-24 03:30:44 UTC (rev 18802)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-07-28 00:19:50 UTC (rev 18803)
@@ -791,6 +791,46 @@
integer :: ispecperio, ispecperio2, ispec2, i2, j2
integer :: iglob_target_to_replace, ispec3, i3, j3
+!<RMODRAK
+ ! NOISE_TOMOGRAPHY = 0 - turn noise tomography subroutines off, resulting in
+ ! an earthquake simulation rather than a noise simulation
+
+ ! NOISE_TOMOGRAPHY = 1 - compute "generating wavefield" and store the result
+
+ ! NOISE_TOMOGRAPHY = 2 - compute "ensemble forward wavefield"; if an adjoint
+ ! simulation is planned, users need to manually extract cross-correlograms
+
+ ! NOISE_TOMOGRAPHY = 3 - carry out adjoint simulation; for noise tomography
+ ! applications, users need to supply adjoint source(s) based on cross-
+ ! -correlograms from previous simulation
+
+ ! For an explanation of terms and concepts in noise tomography, see "Tromp et
+ ! al., 2011, Noise Cross-Correlation Sensitivity Kernels, Geophysical Journal
+ ! International"
+
+
+ integer, parameter :: NOISE_TOMOGRAPHY = 0
+ !setting this flag to a value other than 0 will result in a noise simulation
+
+ !for NOISE_TOMOGRAPHY = 1
+ integer :: irec_master, ispec_noise
+ real(kind=CUSTOM_REAL) :: xi_noise, gamma_noise, angle_noise
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: time_function_noise
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: source_array_noise
+
+ !for NOISE_TOMOGRAPHY = 2
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: mask_noise
+ !to avoid empty arrays depending on SH/P_SV, use separate arrays for x,y,z
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
+ surface_movie_x_noise, surface_movie_y_noise, surface_movie_z_noise
+
+ !for NOISE_TOMOGRAPHY = 3
+ !no additional declarations required for NOISE_TOMOGRAPHY = 3
+
+
+!>RMODRAK
+
+
!! DK DK Feb 2010 for periodic conditions: detect common points between left and right edges
!***********************************************************************
@@ -3527,6 +3567,43 @@
! open the file in which we will store the energy curve
if(output_energy .and. myrank == 0) open(unit=IOUT_ENERGY,file='energy.dat',status='unknown')
+!<RMODRAK
+
+ if (NOISE_TOMOGRAPHY /= 0) then
+
+ !allocate arrays for noise tomography
+ allocate(time_function_noise(NSTEP))
+ allocate(source_array_noise(3,NGLLX,NGLLZ,NSTEP))
+ allocate(mask_noise(nglob))
+ allocate(surface_movie_x_noise(nglob))
+ allocate(surface_movie_y_noise(nglob))
+ allocate(surface_movie_z_noise(nglob))
+
+ !read in parameters for noise tomography
+ call read_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
+ any_acoustic,any_poroelastic,p_sv,irec_master, &
+ Mxx,Mxz,Mzz,factor,NSOURCES, &
+ xi_receiver,gamma_receiver,ispec_selected_rec,nrec, &
+ xi_noise,gamma_noise,ispec_noise,angle_noise)
+
+ endif ! NOISE_TOMOGRAPHY /= 0
+
+
+ if (NOISE_TOMOGRAPHY == 1) then
+ call compute_source_array_noise(p_sv,NSTEP,deltat,nglob,ibool,ispec_noise, &
+ xi_noise,gamma_noise,angle_noise,xigll,zigll, &
+ time_function_noise,source_array_noise)
+
+ elseif (NOISE_TOMOGRAPHY == 2) then
+ call create_mask_noise(p_sv,nglob,coord,angle_noise,mask_noise)
+
+ elseif (NOISE_TOMOGRAPHY == 3) then
+ call create_mask_noise(p_sv,nglob,coord,angle_noise,mask_noise)
+
+ endif
+
+!>RMODRAK
+
!
!---- s t a r t t i m e i t e r a t i o n s
!
@@ -4862,6 +4939,30 @@
endif ! if this processor carries the source and the source element is elastic
enddo ! do i_source=1,NSOURCES
+!<RMODRAK
+
+ ! inject wavefield sources for noise simulations
+
+ if (NOISE_TOMOGRAPHY == 1) then
+ call add_point_source_noise(p_sv,it,NSTEP,nglob,ibool,ispec_noise, &
+ accel_elastic,angle_noise,time_function_noise, &
+ source_array_noise)
+
+ elseif (NOISE_TOMOGRAPHY == 2) then
+ call add_surface_movie_noise(p_sv,it,NSTEP,nglob,ibool,accel_elastic, &
+ surface_movie_x_noise,surface_movie_y_noise, &
+ surface_movie_z_noise,mask_noise)
+
+ elseif (NOISE_TOMOGRAPHY == 3) then
+ call add_surface_movie_noise(p_sv,it,NSTEP,nglob,ibool,b_accel_elastic, &
+ surface_movie_x_noise,surface_movie_y_noise, &
+ surface_movie_z_noise,mask_noise)
+
+ endif
+
+!>RMODRAK
+
+
endif ! if not using an initial field
endif !if(any_elastic)
@@ -6219,6 +6320,13 @@
endif ! if(SIMULATION_TYPE == 2)
+!<RMODRAK
+ if ( NOISE_TOMOGRAPHY == 1 ) then
+ call save_surface_movie_noise(p_sv,it,nglob,displ_elastic)
+ endif
+!>RMODRAK
+
+
!
!---- display results at given time steps
!
More information about the CIG-COMMITS
mailing list