[cig-commits] [commit] devel, master: (SPECFEM2D) BENCHMARKED (725bdb8)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Jun 18 15:19:37 PDT 2014
Repository : https://github.com/geodynamics/specfem2d
On branches: devel,master
Link : https://github.com/geodynamics/specfem2d/compare/fc67e6fd7ad890705b2b72b4b3c509accb22249e...e9ca46c40131588d89d7b0883250bc6584ce6b4c
>---------------------------------------------------------------
commit 725bdb8d99fb8cd0b6ded0380fdf200229300060
Author: Yang Luo <yangl at princeton.edu>
Date: Tue Sep 13 16:06:31 2011 +0000
(SPECFEM2D) BENCHMARKED
>---------------------------------------------------------------
725bdb8d99fb8cd0b6ded0380fdf200229300060
.../Par_file_M2_UPPA => INDUSTRIAL_FORMAT/Par_file | 54 +++--
INDUSTRIAL_FORMAT/README | 14 ++
INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.H | 12 ++
INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.bin | Bin 0 -> 387000 bytes
M2_UPPA/SOURCE_M2_UPPA => INDUSTRIAL_FORMAT/SOURCE | 4 +-
INDUSTRIAL_FORMAT/constants.h | 186 +++++++++++++++++
.../interfaces_M2_UPPA_curved.dat | 11 +-
INDUSTRIAL_FORMAT/interpolate.bash | 15 ++
INDUSTRIAL_FORMAT/interpolate.f90 | 228 +++++++++++++++++++++
INDUSTRIAL_FORMAT/run.bash | 11 +
10 files changed, 498 insertions(+), 37 deletions(-)
diff --git a/M2_UPPA/Par_file_M2_UPPA b/INDUSTRIAL_FORMAT/Par_file
similarity index 79%
copy from M2_UPPA/Par_file_M2_UPPA
copy to INDUSTRIAL_FORMAT/Par_file
index 2ff1f51..fcef543 100644
--- a/M2_UPPA/Par_file_M2_UPPA
+++ b/INDUSTRIAL_FORMAT/Par_file
@@ -1,9 +1,10 @@
# title of job
-title = Test for M2 UPPA
+title = Test Adjoint Pure Acoustic
# forward or adjoint simulation
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
-SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
+NOISE_TOMOGRAPHY = 0 # 0 = earthquake simulation, 1/2/3 = noise simulation
+SAVE_FORWARD = .true. # save the last frame, needed for adjoint simulation
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -13,7 +14,7 @@ ngnod = 9 # number of control nodes per e
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
+READ_EXTERNAL_SEP_FILE = .false. # Read external model 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
@@ -21,8 +22,8 @@ freq0 = 10 # frequency for viscous attenua
p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 1600 # total number of time steps
-deltat = 1.d-3 # duration of a time step
+nt = 10000 # total number of time steps
+deltat = 5.d-4 # duration of a time step
# source parameters
NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
@@ -33,28 +34,28 @@ N_SLS = 2 # number of standard li
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
+seismotype = 6 # record 1=displ 2=veloc 3=accel 4=pressure 6=potential
generate_STATIONS = .true. # 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 = 11 # number of receivers
-xdeb = 300. # first receiver x in meters
-zdeb = 2200. # first receiver z in meters
-xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
-zfin = 2200. # last receiver z in meters (ignored if onlyone receiver)
-enreg_surf_same_vertical = .true. # receivers inside the medium or at the surface
+nrec = 351 # number of receivers
+xdeb = 15520.0 # first receiver x in meters
+zdeb = 80.0 # first receiver z in meters
+xfin = 29520.0 # last receiver x in meters (ignored if onlyone receiver)
+zfin = 80.0 # 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_postscript_snapshot = .false. # 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
+modelvect = .true. # 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)
@@ -73,10 +74,10 @@ nbmodels = 4 # nb of different models
# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
# 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 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
-2 1 2500.d0 2700.d0 0 0 0 9999 9999 0 0 0 0 0 0
-3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
-4 1 2200.d0 2200.d0 1343.375d0 0 0 9999 9999 0 0 0 0 0 0
+1 1 4000.d0 2500.d0 2000.0d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 1000.d0 1500.d0 0 0 0 9999 9999 0 0 0 0 0 0
+3 1 4200.d0 4500.d0 2443.375d0 0 0 9999 9999 0 0 0 0 0 0
+4 1 4200.d0 4200.d0 2343.375d0 0 0 9999 9999 0 0 0 0 0 0
# external mesh or not
read_external_mesh = .false.
@@ -100,22 +101,19 @@ tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the
# PARAMETERS FOR INTERNAL MESHING
# file containing interfaces for internal mesh
-interfacesfile = ../interfaces_M2_UPPA_curved.dat
+interfacesfile = ./interfaces_M2_UPPA_curved.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
+xmax = 59040.d0 # abscissa of right side of the model
+nx = 369 # number of elements along X
# absorbing boundary parameters (see absorbing_conditions above)
-absorbbottom = .true.
+absorbbottom = .false.
absorbright = .true.
-absorbtop = .false.
+absorbtop = .true.
absorbleft = .true.
# define the different regions of the model in the (nx,nz) spectral element mesh
-nbregions = 4 # nb of regions and model number for each
-1 80 1 20 1
-1 80 21 40 2
-1 80 41 60 3
-60 70 21 40 4
+nbregions = 1 # nb of regions and model number for each
+1 369 1 76 2
diff --git a/INDUSTRIAL_FORMAT/README b/INDUSTRIAL_FORMAT/README
new file mode 100644
index 0000000..f7f76f6
--- /dev/null
+++ b/INDUSTRIAL_FORMAT/README
@@ -0,0 +1,14 @@
+1. modify interpolate.f90
+ mainly to change the location/names of your SEP/SU models
+
+2. run interpolation:
+ ./interpolate.bash
+
+3. then the model will be saved in "DATA/model_velocity.dat_input"
+ you can then run simulations in that model by setting in "DATA/Par_file"
+ assign_external_model = .true.
+ READ_EXTERNAL_SEP_FILE = .true.
+
+ a simple example can be done by:
+ ./run.bash
+
diff --git a/INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.H b/INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.H
new file mode 100755
index 0000000..362c7d4
--- /dev/null
+++ b/INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.H
@@ -0,0 +1,12 @@
+in=./vp.bin
+n1=645
+n2=1
+n3=150
+o1=0
+o2=0
+o3=0
+d1=80
+d2=80
+d3=80
+esize=4
+data_format="native_float"
diff --git a/INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.bin b/INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.bin
new file mode 100755
index 0000000..33f246d
Binary files /dev/null and b/INDUSTRIAL_FORMAT/SEG_2D_SALT/vp.bin differ
diff --git a/M2_UPPA/SOURCE_M2_UPPA b/INDUSTRIAL_FORMAT/SOURCE
similarity index 87%
copy from M2_UPPA/SOURCE_M2_UPPA
copy to INDUSTRIAL_FORMAT/SOURCE
index 73cd2c0..4047208 100644
--- a/M2_UPPA/SOURCE_M2_UPPA
+++ b/INDUSTRIAL_FORMAT/SOURCE
@@ -1,7 +1,7 @@
#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 = 2500. # source location x in meters
-zs = 2500. # source location z in meters
+xs = 29520. # source location x in meters
+zs = 80. # 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 = 10.0 # dominant source frequency (Hz) if not Dirac or Heaviside
diff --git a/INDUSTRIAL_FORMAT/constants.h b/INDUSTRIAL_FORMAT/constants.h
new file mode 100644
index 0000000..bb7b60c
--- /dev/null
+++ b/INDUSTRIAL_FORMAT/constants.h
@@ -0,0 +1,186 @@
+!=====================================================================
+!
+! S P E C F E M 2 D Version 6 . 2
+!
+!=====================================================================
+
+! setup/constants.h. Generated from constants.h.in by configure.
+
+!
+! solver in single or double precision depending on the machine (4 or 8 bytes)
+!
+! ALSO CHANGE FILE precision.h ACCORDINGLY
+!
+ integer, parameter :: SIZE_REAL = 4
+ integer, parameter :: SIZE_DOUBLE = 8
+
+
+! set to SIZE_REAL to run in single precision
+! set to SIZE_DOUBLE to run in double precision (increases memory size by 2)
+!
+! DO CHANGE precision.h accordingly
+!
+ integer, parameter :: CUSTOM_REAL = SIZE_DOUBLE
+
+!----------- parameters that can be changed by the user -----------
+
+! number of Gauss-Lobatto-Legendre (GLL) points (i.e., polynomial degree + 1)
+ integer, parameter :: NGLLX = 5
+! the code does NOT work if NGLLZ /= NGLLX because it then cannot handle a non-structured mesh
+! due to non matching polynomial degrees along common edges
+ integer, parameter :: NGLLZ = NGLLX
+
+! further reduce cache misses inner/outer in two passes in the case of an MPI simulation
+! this flag is ignored in the case of a serial simulation
+ logical, parameter :: FURTHER_REDUCE_CACHE_MISSES = .true.
+
+! output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+ logical,parameter :: SU_FORMAT=.true.
+
+! for inverse Cuthill-McKee (1969) permutation
+ logical, parameter :: PERFORM_CUTHILL_MCKEE = .true.
+ logical, parameter :: INVERSE = .true.
+ logical, parameter :: FACE = .false.
+ integer, parameter :: NGNOD_QUADRANGLE = 4
+! perform classical or multi-level Cuthill-McKee ordering
+ logical, parameter :: CMcK_MULTI = .false.
+! maximum size if multi-level Cuthill-McKee ordering
+ integer, parameter :: LIMIT_MULTI_CUTHILL = 50
+
+! implement Cuthill-McKee or replace with identity permutation
+ logical, parameter :: ACTUALLY_IMPLEMENT_PERM_OUT = .false.
+ logical, parameter :: ACTUALLY_IMPLEMENT_PERM_INN = .false.
+ logical, parameter :: ACTUALLY_IMPLEMENT_PERM_WHOLE = .true.
+
+! add MPI barriers and suppress seismograms if we generate traces of the run for analysis with "ParaVer"
+ logical, parameter :: GENERATE_PARAVER_TRACES = .false.
+
+! option to display only part of the mesh and not the whole mesh,
+! for instance to analyze Cuthill-McKee mesh partitioning etc.
+! Possible values are:
+! 1: display all the elements (i.e., the whole mesh)
+! 2: display inner elements only
+! 3: display outer elements only
+! 4: display a fixed number of elements (in each partition) only
+ integer, parameter :: DISPLAY_SUBSET_OPTION = 1
+! number of spectral elements to display in each subset when a fixed subset size is used (option 4 above)
+ integer, parameter :: NSPEC_DISPLAY_SUBSET = 2300
+
+! use this t0 as earliest starting time rather than the automatically calculated one
+! (must be positive and bigger than the automatically one to be effective;
+! simulation will start at t = - t0)
+ double precision, parameter :: USER_T0 = 0.0d0
+
+!--- beginning of Nicolas Le Goff's constants for an unstructured CUBIT/METIS/SCOTCH mesh
+
+! maximum number of neighbors per element
+ integer, parameter :: MAX_NEIGHBORS = 40
+
+! maximum number of elements that can contain the same node
+ integer, parameter :: nsize = 40
+
+!--- end of Nicolas Le Goff's constants for an unstructured CUBIT/METIS/SCOTCH mesh
+
+! output file for energy
+ integer, parameter :: IOUT_ENERGY = 77
+
+! select fast (Paul Fischer) or slow (topology only) global numbering algorithm
+ logical, parameter :: FAST_NUMBERING = .true.
+
+! mesh tolerance for fast global numbering
+ double precision, parameter :: SMALLVALTOL = 0.00001d0
+
+! displacement threshold above which we consider the code became unstable
+ double precision, parameter :: STABILITY_THRESHOLD = 1.d+25
+
+! input and output files
+ integer, parameter :: IIN = 40
+ integer, parameter :: ISTANDARD_OUTPUT = 6
+! uncomment this to write to standard output
+ integer, parameter :: IOUT = ISTANDARD_OUTPUT
+! uncomment this to write to file instead
+! integer, parameter :: IOUT = 41
+
+! number of lines per source in SOURCE file
+ integer, parameter :: NLINES_PER_SOURCE = 13
+
+! flags for absorbing boundaries
+ integer, parameter :: IBOTTOM = 1
+ integer, parameter :: IRIGHT = 2
+ integer, parameter :: ITOP = 3
+ integer, parameter :: ILEFT = 4
+
+! number of edges and corners of each element
+ integer, parameter :: NEDGES = 4
+ integer, parameter :: NCORNERS = 4
+
+! a few useful constants
+ double precision, parameter :: ZERO = 0.d0,ONE = 1.d0
+ double precision, parameter :: HALF = 0.5d0,TWO = 2.d0,QUART = 0.25d0
+
+! pi
+ double precision, parameter :: PI = 3.141592653589793d0
+
+! 4/3
+ double precision, parameter :: FOUR_THIRDS = 4.d0/3.d0
+
+! 1/24
+ double precision, parameter :: ONE_OVER_24 = 1.d0 / 24.d0
+
+! parameters to define the Gauss-Lobatto-Legendre points
+ double precision, parameter :: GAUSSALPHA = ZERO,GAUSSBETA = ZERO
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.d+30,TINYVAL = 1.d-9
+
+! number of spatial dimensions
+ integer, parameter :: NDIM = 2
+
+! maximum length of station and network name for receivers
+ integer, parameter :: MAX_LENGTH_STATION_NAME = 32
+ integer, parameter :: MAX_LENGTH_NETWORK_NAME = 8
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 4
+
+! we mimic a triangle of half duration equal to half_duration_triangle
+! using a Gaussian having a very close shape, as explained in Figure 4.2
+! of the manual. This source decay rate to mimic an equivalent triangle
+! was found by trial and error
+ double precision, parameter :: SOURCE_DECAY_MIMIC_TRIANGLE = 1.628d0
+
+! non linear display to enhance small amplitudes in color images
+ double precision, parameter :: POWER_DISPLAY_COLOR = 0.30d0
+
+! US letter paper or European A4
+ logical, parameter :: US_LETTER = .false.
+
+! X and Z axis origin of PostScript plot in centimeters
+ double precision, parameter :: ORIG_X = 2.4d0
+ double precision, parameter :: ORIG_Z = 2.9d0
+
+! dot to centimeter conversion for PostScript
+ double precision, parameter :: CENTIM = 28.5d0
+
+! parameters for arrows for PostScript snapshot
+ double precision, parameter :: ARROW_ANGLE = 20.d0
+ double precision, parameter :: ARROW_RATIO = 0.40d0
+
+! size of frame used for Postscript display in percentage of the size of the page
+ double precision, parameter :: RPERCENTX = 70.0d0,RPERCENTZ = 77.0d0
+
+! 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
+ integer, parameter :: POROELASTIC_MATERIAL = 3
+
+! file number for interface file
+ integer, parameter :: IIN_INTERFACES = 15
+
+! ignore variable name field (junk) at the beginning of each input line
+ logical, parameter :: IGNORE_JUNK = .true.,DONT_IGNORE_JUNK = .false.
+
diff --git a/Tromp2005/interfaces_Tromp2005_s100.dat b/INDUSTRIAL_FORMAT/interfaces_M2_UPPA_curved.dat
similarity index 72%
copy from Tromp2005/interfaces_Tromp2005_s100.dat
copy to INDUSTRIAL_FORMAT/interfaces_M2_UPPA_curved.dat
index e3e9061..7d281c8 100644
--- a/Tromp2005/interfaces_Tromp2005_s100.dat
+++ b/INDUSTRIAL_FORMAT/interfaces_M2_UPPA_curved.dat
@@ -6,13 +6,10 @@
# interface number 1 (bottom of the mesh)
2
0 0
- 2000 0
-# interface number 5 (topography, top of the mesh)
+ 59040 0
2
- 0 800
- 2000 800
+ 0 11840
+ 59040 11840
#
# for each layer, we give the number of spectral elements in the vertical direction
-#
-# layer number 1
- 32
\ No newline at end of file
+ 76
diff --git a/INDUSTRIAL_FORMAT/interpolate.bash b/INDUSTRIAL_FORMAT/interpolate.bash
new file mode 100755
index 0000000..8a97443
--- /dev/null
+++ b/INDUSTRIAL_FORMAT/interpolate.bash
@@ -0,0 +1,15 @@
+#!/bin/bash -eu
+
+DIR_SEM="../../"
+
+cp Par_file $DIR_SEM/DATA/
+cp interfaces_M2_UPPA_curved.dat $DIR_SEM/DATA/
+cp SOURCE $DIR_SEM/DATA/
+cp constants.h $DIR_SEM/setup/
+
+ifort interpolate.f90 -o $DIR_SEM/xinterpolate
+cd $DIR_SEM
+
+make
+cp bin/* ./
+./xinterpolate
diff --git a/INDUSTRIAL_FORMAT/interpolate.f90 b/INDUSTRIAL_FORMAT/interpolate.f90
new file mode 100644
index 0000000..b00f8b9
--- /dev/null
+++ b/INDUSTRIAL_FORMAT/interpolate.f90
@@ -0,0 +1,228 @@
+ program MAIN
+
+ implicit none
+ !include "constants.h"
+ character(len=512),parameter :: sep_directory='./EXAMPLES/INDUSTRIAL_FORMAT/SEG_2D_SALT/'
+ character(len=512),parameter :: sep_header_file_vp='vp.H'
+ character(len=512),parameter :: sep_header_file_vs='vs.H' ! not used
+ character(len=512),parameter :: sep_header_file_rho='rho.H' ! not used
+ integer :: NX,NY,NZ,esize
+ real :: OX,OY,OZ,DX,DY,DZ
+ real,dimension(:,:),allocatable :: vp_SEP,vs_SEP,rho_SEP
+ character(len=512) :: system_command,sep_file,data_format,mesh_file,wavespeed_file
+
+ ! creating mesh/grid
+ print*, 'generating mesh/grid...'
+ write(system_command,"('sed -e ""s#^SIMULATION_TYPE.*#SIMULATION_TYPE = 1 #g"" < ./DATA/Par_file > temp; mv temp ./DATA/Par_file')"); call system(system_command)
+ write(system_command,"('sed -e ""s#^SAVE_FORWARD.*#SAVE_FORWARD = .false. #g"" < ./DATA/Par_file > temp; mv temp ./DATA/Par_file')"); call system(system_command)
+ write(system_command,"('sed -e ""s#^assign_external_model.*#assign_external_model = .false. #g"" < ./DATA/Par_file > temp; mv temp ./DATA/Par_file')"); call system(system_command)
+ write(system_command,"('sed -e ""s#^READ_EXTERNAL_SEP_FILE.*#READ_EXTERNAL_SEP_FILE = .false. #g"" < ./DATA/Par_file > temp; mv temp ./DATA/Par_file')"); call system(system_command)
+ write(system_command,"('sed -e ""s#^nt.*#nt = 5 #g"" < ./DATA/Par_file > temp; mv temp ./DATA/Par_file')"); call system(system_command)
+ write(system_command,"('./xmeshfem2D')"); call system(system_command)
+ write(system_command,"('./xspecfem2D')"); call system(system_command)
+
+ ! reading SEP models
+ print*, 'reading SEP models...'
+ !!! VP model in SEP format
+ call READ_SEP_HEADER(sep_directory,sep_header_file_vp, &
+ sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ, &
+ esize,data_format)
+ allocate(vp_SEP(NX,NZ),vs_SEP(NX,NZ),rho_SEP(NX,NZ)) ! assuming that VP,VS,RHO models have the same dimension
+ call READ_SEP_MODEL_2D(sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ,esize,data_format,vp_SEP)
+
+ !!! VS model, scaled or in SEP format
+ vs_SEP=0.0
+ !vs_SEP=vp_SEP/1.732
+ !call READ_SEP_HEADER(sep_directory,sep_header_file_vs, &
+ ! sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ, &
+ ! esize,data_format)
+ !call READ_SEP_MODEL_2D(sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ,esize,data_format,vs_SEP)
+
+ !!! RHO model, scaled or in SEP format
+ rho_SEP = 1000.0
+ !rho_SEP = (1.6612 * (vp_SEP/14500*4500 / 1000.0) &
+ ! -0.4720 * (vp_SEP/14500*4500 / 1000.0)**2 &
+ ! +0.0671 * (vp_SEP/14500*4500 / 1000.0)**3 &
+ ! -0.0043 * (vp_SEP/14500*4500 / 1000.0)**4 &
+ ! +0.000106*(vp_SEP/14500*4500 / 1000.0)**5 )*1000.0
+ !call READ_SEP_HEADER(sep_directory,sep_header_file_rho, &
+ ! sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ, &
+ ! esize,data_format)
+ !call READ_SEP_MODEL_2D(sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ,esize,data_format,rho_SEP)
+
+ ! interpolating models
+ print*, 'interpolating models onto mesh/grid...'
+ mesh_file="./DATA/model_velocity.dat_output"
+ wavespeed_file="./DATA/model_velocity.dat_input"
+ call interpolate(vp_SEP,vs_SEP,rho_SEP,NX,NY,NZ,DX,DY,DZ,OX,OY,OZ,mesh_file,wavespeed_file)
+ end program MAIN
+
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+
+ subroutine READ_SEP_HEADER(sep_directory,sep_header_file, &
+ sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ, &
+ esize,data_format)
+
+ implicit none
+
+ !!! input
+ character(len=512),intent(in) :: sep_directory ! where the sep header file and sep file are saved
+ character(len=512),intent(in) :: sep_header_file ! file name of the sep header file
+
+ !!! output
+ character(len=512),intent(out) :: sep_file ! sep file specified in sep header file
+ integer,intent(out) :: NX,NY,NZ ! n1,n2,n3 in sep header file
+ real,intent(out) :: OX,OY,OZ ! o1,o2,o3 in sep header file
+ real,intent(out) :: DX,DY,DZ ! d1,d2,d3 in sep header file
+ integer,intent(out) :: esize ! esize in sep header file
+ character(len=512),intent(out) :: data_format ! data_format in sep header file
+
+ !!! local
+ integer :: ier
+ character(len=512) :: junk,sep_header_file_complete
+
+ sep_header_file_complete=trim(adjustl(sep_directory))//trim(adjustl(sep_header_file))
+
+ open(unit=13,file=trim(adjustl(sep_header_file_complete)),status='old',iostat=ier)
+ print*, ''
+ print*, '*******************************************************************************'
+ print*, 'reading sep header file: '
+ print*, trim(adjustl(sep_header_file_complete))
+ if (ier/=0) stop 'ERROR: cannot open sep header file'
+
+ read(13,'(a3a)') junk, sep_file
+ read(13,'(a3i10)') junk, NX
+ read(13,'(a3i10)') junk, NY
+ read(13,'(a3i10)') junk, NZ
+ read(13,'(a3f20.0)') junk, OX
+ read(13,'(a3f20.0)') junk, OY
+ read(13,'(a3f20.0)') junk, OZ
+ read(13,'(a3f20.0)') junk, DX
+ read(13,'(a3f20.0)') junk, DY
+ read(13,'(a3f20.0)') junk, DZ
+ read(13,'(a6i10)') junk, esize
+ read(13,'(a13a)') junk, data_format
+ close(13)
+ sep_file=trim(adjustl(sep_directory))//trim(adjustl(sep_file))
+ data_format=data_format(1:len_trim(adjustl(data_format))-1)
+
+ print*, ''
+ print*, 'sep file specified in the header file is: ', trim(adjustl(sep_file))
+ print*, 'NX,NY,NZ = ', NX,NY,NZ
+ print*, 'OX,OY,OZ = ', OX,OY,OZ
+ print*, 'DX,DY,DZ = ', DX,DY,DZ
+ print*, 'esize = ', esize
+ print*, 'data_format = ', trim(adjustl(data_format))
+ print*, '*******************************************************************************'
+ print*, ''
+
+ end subroutine READ_SEP_HEADER
+
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+
+ subroutine READ_SEP_MODEL_2D(sep_file,NX,NY,NZ,OX,OY,OZ,DX,DY,DZ,esize,data_format, &
+ model)
+
+ implicit none
+
+ !!! input
+ character(len=512),intent(in) :: sep_file ! sep file specified in sep header file
+ integer,intent(in) :: NX,NY,NZ ! n1,n2,n3 in sep header file
+ real,intent(in) :: OX,OY,OZ ! o1,o2,o3 in sep header file
+ real,intent(in) :: DX,DY,DZ ! d1,d2,d3 in sep header file
+ integer,intent(in) :: esize ! esize in sep header file
+ character(len=512),intent(in) :: data_format ! data_format in sep header file
+
+ !!! output
+ real,dimension(NX,NZ),intent(out) :: model
+
+ !!! local
+ integer :: ier
+
+ ! check whether the model is 2D (NY==1)
+ if (NY/=1) stop 'ERROR: this only works for 2D problems (NY/n2 must be 1)'
+
+ ! note that we keep NY as general in the following (for 3D problems in the future)
+ open(unit=14,file=trim(adjustl(sep_file)),access='direct',status='old',recl=4*NX*NY*NZ,iostat=ier)
+ print*, '*******************************************************************************'
+ print*, 'reading sep file: '
+ print*, trim(adjustl(sep_file))
+ if (ier/=0) stop 'ERROR: cannot open sep file'
+
+ read(14,rec=1,iostat=ier) model(:,:)
+ close(14)
+ if (ier/=0) stop 'ERROR: reading sep file'
+ print*, 'done reading sucessfully'
+ print*, '*******************************************************************************'
+ print*, ''
+
+ end subroutine READ_SEP_MODEL_2D
+
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+
+ subroutine interpolate(vp,vs,rho,NX,NY,NZ,DX,DY,DZ,OX,OY,OZ,mesh_file,wavespeed_file)
+
+ implicit none
+
+ !!! constants
+ integer, parameter :: NGLLX=5,NGLLZ=5
+ double precision, parameter :: vp_water=5500 ! upper bound
+
+ !!! input
+ integer,intent(in) :: NX,NY,NZ
+ real,intent(in) :: OX,OY,OZ
+ real,intent(in) :: DX,DY,DZ
+ real,dimension(NX,NZ),intent(in) :: vp,vs,rho
+ character(len=512),intent(in) :: mesh_file
+ character(len=512),intent(in) :: wavespeed_file
+
+ !!! output
+ ! currently no output, all information is saved in wavespeed_file
+
+ !!! local
+ integer :: ier,ix,iz,i
+ integer, dimension(NGLLX*NGLLZ) :: iglob
+ double precision :: rho_temp,vp_temp,vs_temp
+ double precision, dimension(NGLLX*NGLLZ) :: rho_new,vp_new,vs_new,x,z
+
+
+ open(unit=15,file=trim(adjustl(mesh_file)),status='old')
+ open(unit=16,file=trim(adjustl(wavespeed_file)),status='unknown')
+ ier=0
+ do while (ier==0)
+ do i=1,NGLLX*NGLLZ
+ read(15,'(I10, 5F13.4)',iostat=ier) iglob(i),x(i),z(i),rho_temp,vp_temp,vs_temp
+ ix=NINT((x(i)-OX)/DX)+1
+ iz=NINT((z(i)-OZ)/DZ)+1
+ if (ix>NX-2) ix=NX-2
+ if (ix<1) ix=1
+ if (iz>NZ) iz=NZ
+ if (iz<1) iz=1
+ rho_new(i)=rho(ix,iz)
+ vp_new(i)=vp(ix,iz)
+ vs_new(i)=vs(ix,iz)
+ if(vs_temp<1000.0) vs_new(i)=0.0
+ if(vs_temp<1000.0) rho_new(i)=1000.0
+ enddo
+ if (ier==0) then
+ do i=1,NGLLX*NGLLZ
+ write(16,'(I10, 5F13.4)') iglob(i),x(i),z(i),rho_new(i),vp_new(i),vs_new(i)
+ enddo
+ endif
+ enddo
+
+ close(15)
+ close(16)
+
+ end subroutine interpolate
+
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+!-----------------------------------------------------------------------------------
+
diff --git a/INDUSTRIAL_FORMAT/run.bash b/INDUSTRIAL_FORMAT/run.bash
new file mode 100755
index 0000000..5d425f2
--- /dev/null
+++ b/INDUSTRIAL_FORMAT/run.bash
@@ -0,0 +1,11 @@
+#!/bin/bash -eu
+
+DIR_SEM="../../"
+FILE="$DIR_SEM/DATA/Par_file"
+sed -e "s#^nt.*#nt = 10000 #g" < $FILE > temp; mv temp $FILE
+sed -e "s#^assign_external_model.*#assign_external_model = .true. #g" < $FILE > temp; mv temp $FILE
+sed -e "s#^READ_EXTERNAL_SEP_FILE.*#READ_EXTERNAL_SEP_FILE = .true. #g" < $FILE > temp; mv temp $FILE
+
+cd $DIR_SEM
+./xmeshfem2D
+./xspecfem2D
More information about the CIG-COMMITS
mailing list