Skip to content
Snippets Groups Projects
Commit 8133dee1 authored by joeri.brackenhoff's avatar joeri.brackenhoff
Browse files

3D

parent 5125b6c2
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 3105 additions and 0 deletions
fdelmodc3D/3dfd.c 0 → 100644
+ 1221
− 0
View file @ 8133dee1
#include<mpi.h>
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<time.h>
#include<assert.h>
#include<sys/time.h>
#include"par.h"
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include"fdelmodc3D.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
#define STRIPE_COUNT "4" /* must be an ascii string */
#define STRIPE_SIZE "1048576" /* 1 MB must be an ascii string */
//#define STRIPE_SIZE "268435456" /* 256 MB must be an ascii string */
#define C1 (9.0/8.0)
#define C2 (1.0/24.0)
#define Dx(f,ix,iy,iz,nz) C1*(f[iy+ix*nz+iz] - f[iy+(ix-1)*nz+iz]) - C2*(f[iy+(ix+1)*nz+iz] - f[iy+(ix-2)*nz+iz])
#define Dy(f,ix,iy,iz,nxz) C1*(f[iy*nxz+ix+iz] - f[(iy-1)*nxz+ix+iz]) - C2*(f[(iy+1)*nxz+ix+iz] - f[(iy-2)*nxz+ix+iz])
#define Dz(f,ix,iy,iz) C1*(f[iy+ix+iz] - f[iy+ix+iz-1]) - C2*(f[iy+ix+iz+1] - f[iy+ix+iz-2])
#define Dv(vx,vz,ix,iy,iz,nz,nxz) C1*((vx[iy*nxz+(ix+1)*nz+iz] - vx[iy*nxz+ix*nz+iz]) + \
(vy[(iy+1)*nxz+ix*nz+iz] - vy[iy*nxz+ix*nz+iz]) + \
(vz[iy*nxz+ix*nz+iz+1] - vz[iy*nxz+ix*nz+iz])) - \
C2*((vx[iy*nxz+(ix+2)*nz+iz] - vx[iy*nxz+(ix-1)*nz+iz]) + \
(vy[(iy+2)*nxz+ix*nz+iz] - vy[(iy-1)*nxz+ix*nz+iz]) + \
(vz[iy*nxz+ix*nz+iz+2] - vz[iy*nxz+ix*nz+iz-1]))
int getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *src, shotPar *shot, bndPar *bnd, int verbose);
int readModel3D(modPar mod, bndPar bnd, float *rox, float *roz, float *l2m);
void vinit();
int updateVelocitiesHalo(float *vx, float *vy, float *vz, float *p, float *ro, int halo, int npx, int npy, int npz);
int updateVelocities(float *vx, float *vy, float *vz, float *p, float *ro, int halo, int npx, int npy, int npz);
int updatePressureHalo(float *vx, float *vy, float *vz, float *p, float *l2m, int halo, int npx, int npy, int npz);
int updatePressure(float *vx, float *vy, float *vz, float *p, float *l2m, int halo, int npx, int npy, int npz);
int exchangeHalo(float *leftRecv, float *leftSend, float *rightRecv, float *rightSend, int size, int leftrank, int rightrank, MPI_Request *reqRecv, MPI_Request *reqSend, int tag);
int newHaloVxVz(float *vx, float *vz, int npx, int npy, int npz, int halo, float *leftRecv, float *rightRecv, float *frontRecv, float *backRecv, float *topRecv, float *bottomRecv);
int newHaloP(float *p, int npx, int npy, int npz, int halo, float *leftRecv, float *rightRecv, float *frontRecv, float *backRecv, float *topRecv, float *bottomRecv);
int copyHaloVxVz(float *vx, float *vz, int npx, int npy, int npz, int halo, float *leftSend, float *rightSend, float *frontSend, float *backSend, float *topSend, float *bottomSend);
int copyHaloP(float *p, int npx, int npy, int npz, int halo, float *leftSend, float *rightSend, float *frontSend, float *backSend, float *topSend, float *bottomSend);
int waitForHalo(MPI_Request *reqRecv, MPI_Request *reqSend);
float gauss2time(float t, float f, float t0);
double wallclock_time(void);
void name_ext(char *filename, char *extension);
/* Self documentation */
char *sdoc[] = {
" ",
" file_rcv=recv.su .. base name for receiver files",
" file_snap=snap.su . base name for snapshot files",
" nx=256 ............ number of samples in x-direction",
" ny=nx ............. number of samples in y-direction",
" nz=nx ............. number of samples in z-direction",
" dx=5 .............. spatial sampling in x-direction",
" dy=5 .............. spatial sampling in y-direction",
" dz=5 .............. spatial sampling in z-direction",
"" ,
" verbose=0 ......... silent mode; =1: display info",
" ",
" Jan Thorbecke 2016",
" Cray / TU Delft",
" E-mail: janth@xs4all.nl ",
"",
NULL};
int main (int argc, char *argv[])
{
modPar mod;
recPar rec;
snaPar sna;
wavPar wav;
srcPar src;
bndPar bnd;
shotPar shot;
float *wavelet;
int nx, ny, nz, dims[3], period[3], reorder, coord[3], ndims=3;
int npx, npy, npz, halo, nt;
int my_rank, size, source, dest, snapwritten;
int left, right, front, back, top, bottom;
int direction, displ, halosizex, halosizey, halosizez;
int ioXx, ioXz, ioYx, ioYz, ioZz, ioZx, ioPx, ioPz;
int it, ix, iy, iz, iyp, ixp, izp, isrc, ixsrc, iysrc, izsrc, c1, c2;
int sizes[3], subsizes[3], starts[3];
int gsizes[3], gsubsizes[3], gstarts[3];
int error, rc, verbose;
float fx, fy, fz, dx, dy, dz, flx, fly, flz;
float *p, *vx, *vy, *vz, *rox, *roz, *roy, *l2m, hcp, hro, fac;
float *leftRecv, *leftSend, *rightRecv, *rightSend;
float *frontRecv, *frontSend, *backRecv, *backSend;
float *topRecv, *topSend, *bottomRecv, *bottomSend;
float dt, src_ampl, fmax, fpeaksrc, t0src, time, snaptime;
double t00, t0, t1, t2, tcomm, tcomp, thcomp, tcopy, ttot, tio;
char err_buffer[MPI_MAX_ERROR_STRING];
int resultlen;
MPI_Comm COMM_CART;
MPI_Request reqSend[6], reqRecv[6];
MPI_Status status[12];
MPI_Datatype local_array, global_array;
MPI_Offset disp;
MPI_Info fileinfo;
MPI_File fh;
char filename[1000], *file_snap, *file_rcv;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
MPI_Comm_size(MPI_COMM_WORLD, &size);
vinit();
t0= wallclock_time();
initargs(argc,argv);
requestdoc(0);
if (!getparint("verbose",&verbose)) verbose=0;
if (!getparstring("file_snap",&file_snap)) file_snap="snap.su";
if (!getparstring("file_rcv",&file_rcv)) file_rcv="recv.su";
getParameters3D(&mod, &rec, &sna, &wav, &src, &shot, &bnd, verbose);
/* divide 3D cube among available processors */
dims[0]=0;
dims[1]=0;
dims[2]=0;
MPI_Dims_create(size, ndims, dims);
/* dims contains the number of MPI-tasks in each direction */
/* set number of grid points based on number of procs in dims */
if (!getparint("nx",&nx)) nx=256;
if (!getparint("ny",&ny)) ny=nx;
if (!getparint("nz",&nz)) nz=nx;
if (!getparfloat("dx",&dx)) dx=5.;
if (!getparfloat("dy",&dy)) dy=5.;
if (!getparfloat("dz",&dz)) dz=5.;
halo = 2;
/* new larger dimensions to fit with the domain-decomposition */
nz=dims[0]*ceil(mod.naz/dims[0]);
nx=dims[1]*ceil(mod.nax/dims[1]);
ny=dims[2]*ceil(mod.nay/dims[2]);
// dt=0.001;
nt=4096;
t0src=0.50;
hcp=2000.0;
hro=1800.0;
tcomm=tcomp=thcomp=tcopy=tio=0.0;
/* for stability 10 points per wavelenght */
fmax=hcp/(mod.dx*8);
dt=0.295*mod.dx/hcp;
fpeaksrc=0.2*fmax; /* Ricker wavelet has peak at 1/3 of fmax */
fac = mod.dt/mod.dx;
fx=-mod.dx*nx/2; fy=-mod.dy*ny/2; fz=0;
npz = 2*halo+nz/dims[0];
npx = 2*halo+nx/dims[1];
npy = 2*halo+ny/dims[2];
wavelet = (float *)calloc(nt,sizeof(float));
/* find which MPI-task has the source position */
snaptime = t0src+1.80*npx*dx*0.5/hcp;
snapwritten=0;
nt = (int) 1.1*(t0src+snaptime)/dt;
nt = (int) (t0src+1.5)/dt;
if (verbose && my_rank==0) {
fprintf(stderr,"fmax=%f fpeak=%f dt=%e\n", fmax, fpeaksrc, dt);
fprintf(stderr,"nx=%d nprocx=%d ny=%d nprocy=%d nz=%d nprocz=%d\n", nx, dims[1], ny, dims[2], nz, dims[0]);
fprintf(stderr,"npx=%d npy=%d npz=%d nt=%d time=%f\n", npx, npy, npz, nt, nt*dt);
fprintf(stderr,"source expected at local boundary at %f seconds\n", npx*dx*0.5/hcp);
fflush(stderr);
}
if (my_rank==0) fprintf(stderr,"memory per MPI task is %ld MB\n", (6*npx*npy*npz*4/(1024*1024)));
/* allocate wavefields and medium properties for local grid domains */
p = (float *)calloc(npx*npy*npz,sizeof(float));
vx = (float *)calloc(npx*npy*npz,sizeof(float));
vy = (float *)calloc(npx*npy*npz,sizeof(float));
vz = (float *)calloc(npx*npy*npz,sizeof(float));
/* read velocity and density files */
readModel3D(mod, bnd, rox, roz, l2m);
/* for 2.5 D model npy=1 */
rox = (float *)calloc(npx*npy*npz,sizeof(float));
roy= (float *)calloc(npx*npy*npz,sizeof(float));
roz= (float *)calloc(npx*npy*npz,sizeof(float));
l2m = (float *)calloc(npx*npy*npz,sizeof(float));
/* define homogeneus model */
for (ix=0; ix<npx*1*npz; ix++) {
rox[ix] = fac/hro;
roy[ix] = fac/hro;
roz[ix] = fac/hro;
l2m[ix] = fac*hcp*hcp*hro;
}
/* create cartesian domain decomposition */
period[0]=0;
period[1]=0;
period[2]=0;
reorder=0;
MPI_Cart_create(MPI_COMM_WORLD, 3, dims, period, reorder, &COMM_CART);
/* find out coordinates of the rank */
MPI_Cart_coords(COMM_CART, my_rank, 3, coord);
flz = fz+(dz*nz/dims[0])*coord[0];
flx = fx+(dx*nx/dims[1])*coord[1];
fly = fy+(dy*ny/dims[2])*coord[2];
if (verbose>=2) fprintf(stderr,"Rank %d coordinates are %d %d %d orig=(%5.2F, %5.2f, %5.2f) \n", my_rank, coord[0], coord[1], coord[2], flx, fly, flz);
fflush(stderr);
/* find out neighbours of the rank, MPI_PROC_NULL is a hard boundary of the model */
displ=1;
MPI_Cart_shift(COMM_CART, 1, 1, &left, &right);
MPI_Cart_shift(COMM_CART, 2, 1, &top, &bottom);
MPI_Cart_shift(COMM_CART, 0, 1, &front, &back);
if (verbose>=2) fprintf(stderr, "Rank %d in direction 0 has LR neighbours %d %d FB %d %d TB %d %d\n", my_rank, left, right, front, back, top, bottom);
fflush(stderr);
/* allocate of halo areas */
halosizex = npy*npz*halo;
leftRecv = (float *)calloc(3*halosizex,sizeof(float));
rightRecv = (float *)calloc(3*halosizex,sizeof(float));
leftSend = (float *)calloc(3*halosizex,sizeof(float));
rightSend = (float *)calloc(3*halosizex,sizeof(float));
halosizey = npx*npz*halo;
frontRecv = (float *)calloc(3*halosizey,sizeof(float));
backRecv = (float *)calloc(3*halosizey,sizeof(float));
frontSend = (float *)calloc(3*halosizey,sizeof(float));
backSend = (float *)calloc(3*halosizey,sizeof(float));
halosizez = npy*npx*halo;
topRecv = (float *)calloc(3*halosizez,sizeof(float));
bottomRecv = (float *)calloc(3*halosizez,sizeof(float));
topSend = (float *)calloc(3*halosizez,sizeof(float));
bottomSend = (float *)calloc(3*halosizez,sizeof(float));
if (my_rank==0) fprintf(stderr,"memory per MPI task for halo exchange is %ld MB\n", ((12*(halosizex+halosizey+halosizez))*4/(1024*1024)));
/* create subarrays(excluding halo areas) to write to file with MPI-IO */
/* data in the local array */
sizes[0]=npz;
sizes[1]=npx;
sizes[2]=npy;
subsizes[0]=sizes[0]-2*halo;
subsizes[1]=sizes[1]-2*halo;
subsizes[2]=sizes[2]-2*halo;
starts[0]=halo;
starts[1]=halo;
starts[2]=halo;
MPI_Type_create_subarray(3, sizes, subsizes, starts, MPI_ORDER_C,
MPI_FLOAT, &local_array);
MPI_Type_commit(&local_array);
/* data in the global array */
gsizes[0]=nz;
gsizes[1]=nx;
gsizes[2]=ny;
gsubsizes[0]=subsizes[0];
gsubsizes[1]=subsizes[1];
gsubsizes[2]=subsizes[2];
gstarts[0]=subsizes[0]*coord[0];
gstarts[1]=subsizes[1]*coord[1];
gstarts[2]=subsizes[2]*coord[2];
MPI_Type_create_subarray(3, gsizes, gsubsizes, gstarts, MPI_ORDER_C,
MPI_FLOAT, &global_array);
MPI_Type_commit(&global_array);
/* compute field of the inner grid excluding halo areas */
ioXx=2;
ioXz=ioXx-1;
ioYx=2;
ioYz=ioYx-1;
ioZz=2;
ioZx=ioZz-1;
ioPx=1;
ioPz=ioPx;
t00 = wallclock_time();
for (it=0; it<nt; it++) {
time = it*dt;
wavelet[it] = gauss2time(time,fpeaksrc,t0src);
}
if (my_rank==0) {
FILE *fp;
fp = fopen("src.bin", "w+");
fwrite( wavelet, sizeof(float), nt, fp);
fflush(fp);
fclose(fp);
}
/*
nt =1;
sprintf(filename,"snap_nz%d_nx%d_ny%d.bin",nz, nx, ny);
for (ix=0; ix<npx*npy*npz; ix++) {
p[ix] = my_rank;
}
MPI_Info_create(&fileinfo);
MPI_Info_set(fileinfo, "striping_factor", STRIPE_COUNT);
MPI_Info_set(fileinfo, "striping_unit", STRIPE_SIZE);
MPI_File_delete(filename, MPI_INFO_NULL);
rc = MPI_File_open(MPI_COMM_WORLD, filename, MPI_MODE_RDWR|MPI_MODE_CREATE, fileinfo, &fh);
if (rc != MPI_SUCCESS) {
fprintf(stderr, "could not open input file\n");
MPI_Abort(MPI_COMM_WORLD, 2);
}
disp = 0;
rc = MPI_File_set_view(fh, disp, MPI_FLOAT, global_array, "native", fileinfo);
if (rc != MPI_SUCCESS) {
fprintf(stderr, "error setting view on results file\n");
MPI_Abort(MPI_COMM_WORLD, 4);
}
rc = MPI_File_write_all(fh, p, 1, local_array, status);
if (rc != MPI_SUCCESS) {
MPI_Error_string(rc,err_buffer,&resultlen);
fprintf(stderr,err_buffer);
MPI_Abort(MPI_COMM_WORLD, 5);
}
MPI_File_close(&fh);
*/
/* Main loop over the number of time steps */
for (it=0; it<nt; it++) {
t0 = wallclock_time();
time = it*dt;
//fprintf(stderr,"modeling time step %d for time %f\n", it, time);
/* define source wavelet */
wavelet[it] = gauss2time(time,fpeaksrc,t0src);
/* update of grid values on halo areas */
updateVelocitiesHalo(vx, vy, vz, p, rox, halo, npx, npy, npz);
t1 = wallclock_time();
thcomp += t1-t0;
/* copy of halo areas */
copyHaloVxVz(vx, vz, npx, npy, npz, halo, leftSend, rightSend, frontSend, backSend, topSend, bottomSend);
t2 = wallclock_time();
tcopy += t2-t1;
/* start a-synchronous communication of halo areas to neighbours */
/* this is done first for Vx,Vz fields only */
exchangeHalo(leftRecv, leftSend, rightRecv, rightSend, 2*halosizex, left, right, &reqRecv[0], &reqSend[0], 0);
exchangeHalo(frontRecv, frontSend, backRecv, backSend, 2*halosizey, front, back, &reqRecv[2], &reqSend[2], 4);
exchangeHalo(topRecv, topSend, bottomRecv, bottomSend, 2*halosizez, top, bottom, &reqRecv[4], &reqSend[4], 8);
t1 = wallclock_time();
tcomm += t1-t2;
/* compute update on grid values excluding halo areas */
updateVelocities(vx, vy, vz, p, rox, halo, npx, npy, npz);
t2 = wallclock_time();
tcomp += t2-t1;
/* wait for Vx.Vz halo exchange */
waitForHalo(&reqRecv[0], &reqSend[0]);
t1 = wallclock_time();
tcomm += t1-t2;
/* copy of halo areas back to arrays */
newHaloVxVz(vx, vz, npx, npy, npz, halo, leftRecv, rightRecv, frontRecv, backRecv, topRecv, bottomRecv);
t2 = wallclock_time();
tcopy += t2-t1;
/* add Force source on the Vz grid */
src_ampl = wavelet[it];
/* check if source position is in local domain */
/* for the moment place a source in the middle of each domain */
ixsrc = npx/2;
iysrc = npy/2;
izsrc = npz/2;
isrc = iysrc*npx*npz+ixsrc*npz+izsrc;
// fprintf(stderr,"npz=%d npx=%d npy=%d isrc=%d\n", npz, npx, npy, isrc);
/* source scaling factor to compensate for discretisation */
src_ampl *= rox[isrc]*l2m[isrc]/(dt);
/* Force source */
//if (my_rank == 0) vz[isrc] += 0.25*src_ampl*ro[isrc]*dz;
vz[isrc] += 0.25*src_ampl*rox[isrc]*dz;
/* compute field on the grid of the halo areas */
updatePressureHalo(vx, vy, vz, p, l2m, halo, npx, npy, npz);
t1 = wallclock_time();
thcomp += t1-t2;
/* copy p-field and sent to neighbours */
copyHaloP(p, npx, npy, npz, halo, leftSend, rightSend, frontSend, backSend, topSend, bottomSend);
exchangeHalo(leftRecv, leftSend, rightRecv, rightSend, halosizex, left, right, &reqRecv[0], &reqSend[0], 0);
exchangeHalo(frontRecv, frontSend, backRecv, backSend, halosizey, front, back, &reqRecv[2], &reqSend[2], 4);
exchangeHalo(topRecv, topSend, bottomRecv, bottomSend, halosizez, top, bottom, &reqRecv[4], &reqSend[4], 8);
t2 = wallclock_time();
tcomm += t2-t1;
/* compute update on grid values excluding halo areas */
updatePressure(vx, vy, vz, p, l2m, halo, npx, npy, npz);
t1 = wallclock_time();
tcomp += t1-t2;
/* wait for P halo exchange */
waitForHalo(&reqRecv[0], &reqSend[0]);
t2 = wallclock_time();
tcomm += t2-t1;
newHaloP(p, npx, npy, npz, halo, leftRecv, rightRecv, frontRecv, backRecv, topRecv, bottomRecv);
t1 = wallclock_time();
tcopy += t1-t2;
// fprintf(stderr,"rank %d did time step %d in %f seconds\n", my_rank, it, t1-t0);
// fflush(stderr);
/* write snapshots to file */
// if (time >= snaptime && !snapwritten) {
if ((it+1)%100==0 ) {
t1 = wallclock_time();
sprintf(filename,"snap_nz%d_nx%d_ny%d_it%4d.bin",nz, nx, ny, it);
MPI_Info_create(&fileinfo);
MPI_Info_set(fileinfo, "striping_factor", STRIPE_COUNT);
MPI_Info_set(fileinfo, "striping_unit", STRIPE_SIZE);
MPI_File_delete(filename, MPI_INFO_NULL);
rc = MPI_File_open(MPI_COMM_WORLD, filename, MPI_MODE_RDWR|MPI_MODE_CREATE, fileinfo, &fh);
if (rc != MPI_SUCCESS) {
fprintf(stderr, "could not open input file\n");
MPI_Abort(MPI_COMM_WORLD, 2);
}
disp = 0;
rc = MPI_File_set_view(fh, disp, MPI_FLOAT, global_array, "native", fileinfo);
if (rc != MPI_SUCCESS) {
fprintf(stderr, "error setting view on results file\n");
MPI_Abort(MPI_COMM_WORLD, 4);
}
rc = MPI_File_write_all(fh, p, 1, local_array, status);
if (rc != MPI_SUCCESS) {
MPI_Error_string(rc,err_buffer,&resultlen);
fprintf(stderr,err_buffer);
MPI_Abort(MPI_COMM_WORLD, 5);
}
MPI_File_close(&fh);
/* MPI_Info_create(&fileinfo);
MPI_File_delete(filename, MPI_INFO_NULL);
MPI_File_open(MPI_COMM_WORLD, filename, MPI_MODE_RDWR|MPI_MODE_CREATE, MPI_INFO_NULL, &fh);
MPI_File_set_view(fh, 0, MPI_FLOAT, global_array, "native", MPI_INFO_NULL);
MPI_File_write_all(fh, p, npz*npx*npy, local_array, status);
MPI_File_close(&fh);
*/
snapwritten+=1;
t2 = wallclock_time();
tio += t2-t1;
}
}
ttot = wallclock_time() - t00;
if (my_rank == 0) {
fprintf(stderr,"rank %d total time in %f seconds\n", my_rank, ttot);
fprintf(stderr,"rank %d comm time in %f seconds\n", my_rank, tcomm);
fprintf(stderr,"rank %d comp time in %f seconds\n", my_rank, tcomp);
fprintf(stderr,"rank %d hcomp time in %f seconds\n", my_rank, thcomp);
fprintf(stderr,"rank %d copy time in %f seconds\n", my_rank, tcopy);
fprintf(stderr,"rank %d io time in %f seconds\n", my_rank, tio);
fprintf(stderr,"rank %d snaphsots written to file\n", snapwritten);
}
MPI_Finalize();
return 0;
}
int updateVelocities(float *vx, float *vy, float *vz, float *p, float *ro, int halo, int npx, int npy, int npz)
{
int ix, iy, iz, iyp, ixp, izp, c1, c2, nxz;
int ixs, ixe, iys, iye, izs, ize;
float DpDx, DpDy, DpDz;
nxz=npx*npz;
c1 = 9.0/8.0;
c2 = -1.0/24.0;
ixs=2*halo; ixe=npx-2*halo;
iys=2*halo; iye=npy-2*halo;
izs=2*halo; ize=npz-2*halo;
/* calculate vx,vy,vz for all grid points except on the virtual boundary*/
#pragma omp for private (iy, ix, iz) nowait
#pragma ivdep
for (iy=iys; iy<iye; iy++) {
iyp=iy*nxz;
for (ix=ixs; ix<ixe; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=izs; iz<ize; iz++) {
DpDx = Dx(p,ix,iyp,iz,npz);
DpDy = Dy(p,ixp,iy,iz,nxz);
DpDz = Dz(p,ixp,iyp,iz);
vz[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDz;
vx[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDx;
vy[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDy;
}
}
}
return 0;
}
int updateVelocitiesHalo(float *vx, float *vy, float *vz, float *p, float *ro, int halo, int npx, int npy, int npz)
{
int ix, iy, iz, iyp, ixp, izp, c1, c2, nxz;
float DpDx, DpDy, DpDz;
nxz=npx*npz;
c1 = 9.0/8.0;
c2 = -1.0/24.0;
/* calculate vx,vy,vz for all halo grid points */
/* compute halo areas at left side */
#pragma omp for private (iy, ix, iz) nowait
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=halo; ix<2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
DpDx = Dx(p,ix,iyp,iz,npz);
DpDy = Dy(p,ixp,iy,iz,nxz);
DpDz = Dz(p,ixp,iyp,iz);
vz[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDz;
vx[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDx;
vy[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDy;
}
}
}
/* compute halo areas at right side */
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=npx-2*halo; ix<npx-halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
DpDx = Dx(p,ix,iyp,iz,npz);
DpDy = Dy(p,ixp,iy,iz,nxz);
DpDz = Dz(p,ixp,iyp,iz);
vz[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDz;
vx[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDx;
vy[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDy;
}
}
}
/* compute halo areas at front side */
for (iy=halo; iy<2*halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
DpDx = Dx(p,ix,iyp,iz,npz);
DpDy = Dy(p,ixp,iy,iz,nxz);
DpDz = Dz(p,ixp,iyp,iz);
vz[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDz;
vx[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDx;
vy[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDy;
}
}
}
/* compute halo areas at back side */
for (iy=npy-2*halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
DpDx = Dx(p,ix,iyp,iz,npz);
DpDy = Dy(p,ixp,iy,iz,nxz);
DpDz = Dz(p,ixp,iyp,iz);
vz[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDz;
vx[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDx;
vy[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDy;
}
}
}
/* compute halo areas at top side */
for (iy=2*halo; iy<npy-2*halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<2*halo; iz++) {
DpDx = Dx(p,ix,iyp,iz,npz);
DpDy = Dy(p,ixp,iy,iz,nxz);
DpDz = Dz(p,ixp,iyp,iz);
vz[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDz;
vx[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDx;
vy[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDy;
}
}
}
/* compute halo areas at bottom side */
for (iy=2*halo; iy<npy-2*halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=npz-2*halo; iz<npz-halo; iz++) {
DpDx = Dx(p,ix,iyp,iz,npz);
DpDy = Dy(p,ixp,iy,iz,nxz);
DpDz = Dz(p,ixp,iyp,iz);
vz[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDz;
vx[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDx;
vy[iyp+ixp+iz] += ro[iyp+ixp+iz]*DpDy;
}
}
}
return 0;
}
int updatePressure(float *vx, float *vy, float *vz, float *p, float *l2m, int halo, int npx, int npy, int npz)
{
int ix, iy, iz, iyp, ixp, izp, c1, c2, nxz;
int ixs, ixe, iys, iye, izs, ize;
nxz=npx*npz;
c1 = 9.0/8.0;
c2 = -1.0/24.0;
ixs=2*halo; ixe=npx-2*halo;
iys=2*halo; iye=npy-2*halo;
izs=2*halo; ize=npz-2*halo;
/* calculate p/tzz for all grid points except on the virtual boundary */
#pragma omp for private (ix, iz)
#pragma ivdep
for (iy=iys; iy<iye; iy++) {
iyp=iy*nxz;
for (ix=ixs; ix<ixe; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=izs; iz<ize; iz++) {
p[iyp+ixp+iz] += l2m[iyp+ixp+iz]*(Dv(vx,vz,ix,iy,iz,npz,nxz));
}
}
}
return 0;
}
int updatePressureHalo(float *vx, float *vy, float *vz, float *p, float *l2m, int halo, int npx, int npy, int npz)
{
int ix, iy, iz, iyp, ixp, izp, c1, c2, nxz;
nxz=npx*npz;
c1 = 9.0/8.0;
c2 = -1.0/24.0;
/* calculate p/tzz for all grid points except on the virtual boundary */
/* compute halo areas at left side */
#pragma omp for private (iy, ix, iz) nowait
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=halo; ix<2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] += l2m[iyp+ixp+iz]*(Dv(vx,vz,ix,iy,iz,npz,nxz));
}
}
}
/* compute halo areas at right side */
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=npx-2*halo; ix<npx-halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] += l2m[iyp+ixp+iz]*(Dv(vx,vz,ix,iy,iz,npz,nxz));
}
}
}
/* compute halo areas at front side */
for (iy=halo; iy<2*halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] += l2m[iyp+ixp+iz]*(Dv(vx,vz,ix,iy,iz,npz,nxz));
}
}
}
/* compute halo areas at back side */
for (iy=npy-2*halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] += l2m[iyp+ixp+iz]*(Dv(vx,vz,ix,iy,iz,npz,nxz));
}
}
}
/* compute halo areas at top side */
for (iy=2*halo; iy<npy-2*halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<2*halo; iz++) {
p[iyp+ixp+iz] += l2m[iyp+ixp+iz]*(Dv(vx,vz,ix,iy,iz,npz,nxz));
}
}
}
/* compute halo areas at bottom side */
for (iy=2*halo; iy<npy-2*halo; iy++) {
iyp=iy*nxz;
for (ix=2*halo; ix<npx-2*halo; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=npz-2*halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] += l2m[iyp+ixp+iz]*(Dv(vx,vz,ix,iy,iz,npz,nxz));
}
}
}
return 0;
}
int exchangeHalo(float *leftRecv, float *leftSend, float *rightRecv, float *rightSend, int size, int leftrank, int rightrank, MPI_Request *reqRecv, MPI_Request *reqSend, int tag)
{
int error, my_rank, ltag;
MPI_Status status;
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
if (leftrank != MPI_PROC_NULL) {
ltag = tag;
error = MPI_Irecv(leftRecv, size, MPI_FLOAT, leftrank, ltag, MPI_COMM_WORLD, &reqRecv[0]);
assert (error == MPI_SUCCESS);
// fprintf(stderr,"rank %d recv data from %d left\n", my_rank, leftrank);
ltag = tag+1;
error = MPI_Isend(leftSend, size, MPI_FLOAT, leftrank, ltag, MPI_COMM_WORLD, &reqSend[0]);
assert (error == MPI_SUCCESS);
// fprintf(stderr,"rank %d send data to %d left\n", my_rank, leftrank);
}
else {
reqRecv[0] = MPI_REQUEST_NULL;
reqSend[0] = MPI_REQUEST_NULL;
}
if (rightrank != MPI_PROC_NULL) {
ltag = tag+1;
error = MPI_Irecv(rightRecv, size, MPI_FLOAT, rightrank, ltag, MPI_COMM_WORLD, &reqRecv[1]);
// fprintf(stderr,"rank %d recv data from %d right\n", my_rank, rightrank);
assert (error == MPI_SUCCESS);
ltag = tag;
error = MPI_Isend(rightSend, size, MPI_FLOAT, rightrank, ltag, MPI_COMM_WORLD, &reqSend[1]);
assert (error == MPI_SUCCESS);
// fprintf(stderr,"rank %d send data to %d right\n", my_rank, rightrank);
}
else {
reqRecv[1] = MPI_REQUEST_NULL;
reqSend[1] = MPI_REQUEST_NULL;
}
return 0;
}
int waitForHalo(MPI_Request *reqRecv, MPI_Request *reqSend)
{
int i;
MPI_Status status;
int error;
for (i=0; i<6; i++) {
error = MPI_Wait(&reqSend[i], &status);
assert (error == MPI_SUCCESS);
}
// MPI_Barrier(MPI_COMM_WORLD);
for (i=0; i<6; i++) {
error = MPI_Wait(&reqRecv[i], &status);
assert (error == MPI_SUCCESS);
}
return 0;
}
int copyHaloVxVz(float *vx, float *vz, int npx, int npy, int npz, int halo, float *leftSend, float *rightSend, float *frontSend, float *backSend, float *topSend, float *bottomSend)
{
int ix, iy, iz, ih, iyp, ixp, halosizex, halosizey, halosizez, nxz;
nxz = npx*npz;
/* copy halo areas at left side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=halo; ix<2*halo; ix++) {
ixp=ix*npz;
ih=(ix-halo)*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
leftSend[iy*npz*halo+ih+iz] = vx[iyp+ixp+iz];
leftSend[halosizex+iy*npz*halo+ih+iz] = vz[iyp+ixp+iz];
}
}
}
/* copy halo areas at right side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=npx-2*halo; ix<npx-halo; ix++) {
ixp=ix*npz;
ih=(ix-(npx-2*halo))*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
rightSend[iy*npz*halo+ih+iz] = vx[iyp+ixp+iz];
rightSend[halosizex+iy*npz*halo+ih+iz] = vz[iyp+ixp+iz];
}
}
}
/* copy halo areas at front side */
halosizey = npx*npz*halo;
for (iy=halo; iy<2*halo; iy++) {
iyp=iy*nxz;
ih=(iy-halo)*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
frontSend[ih+ixp+iz] = vx[iyp+ixp+iz];
frontSend[halosizey+ih+ixp+iz] = vz[iyp+ixp+iz];
}
}
}
/* copy halo areas at back side */
for (iy=npy-2*halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
ih=(iy-(npy-2*halo))*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
backSend[ih+ixp+iz] = vx[iyp+ixp+iz];
backSend[halosizey+ih+ixp+iz] = vz[iyp+ixp+iz];
}
}
}
/* copy halo areas at top side */
halosizez = npy*npx*halo;
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<2*halo; iz++) {
ih=iz-halo;
topSend[iy*npx*halo+ix*halo+ih] = vx[iyp+ixp+iz];
topSend[halosizez+iy*npx*halo+ix*halo+ih] = vz[iyp+ixp+iz];
}
}
}
/* copy halo areas at bottom side */
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=npz-2*halo; iz<npz-halo; iz++) {
ih=(iz-(npz-2*halo));
bottomSend[iy*npx*halo+ix*halo+ih] = vx[iyp+ixp+iz];
bottomSend[halosizez+iy*npx*halo+ix*halo+ih] = vz[iyp+ixp+iz];
}
}
}
return 0;
}
int copyHaloP(float *p, int npx, int npy, int npz, int halo, float *leftSend, float *rightSend, float *frontSend, float *backSend, float *topSend, float *bottomSend)
{
int ix, iy, iz, ih, iyp, ixp, halosizex, halosizey, halosizez, nxz;
nxz = npx*npz;
/* copy halo areas at left side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=halo; ix<2*halo; ix++) {
ixp=ix*npz;
ih=(ix-halo)*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
leftSend[iy*npz*halo+ih+iz] = p[iyp+ixp+iz];
}
}
}
/* copy halo areas at right side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=npx-2*halo; ix<npx-halo; ix++) {
ixp=ix*npz;
ih=(ix-(npx-2*halo))*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
rightSend[iy*npz*halo+ih+iz] = p[iyp+ixp+iz];
}
}
}
/* copy halo areas at front side */
halosizey = npx*npz*halo;
for (iy=halo; iy<2*halo; iy++) {
iyp=iy*nxz;
ih=(iy-halo)*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
frontSend[ih+ixp+iz] = p[iyp+ixp+iz];
}
}
}
/* copy halo areas at back side */
for (iy=npy-2*halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
ih=(iy-(npy-2*halo))*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
backSend[ih+ixp+iz] = p[iyp+ixp+iz];
}
}
}
/* copy halo areas at top side */
halosizez = npy*npx*halo;
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<2*halo; iz++) {
ih=iz-halo;
topSend[iy*npx*halo+ix*halo+ih] = p[iyp+ixp+iz];
}
}
}
/* copy halo areas at bottom side */
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=npz-2*halo; iz<npz-halo; iz++) {
ih=(iz-(npz-2*halo));
bottomSend[iy*npx*halo+ix*halo+ih] = p[iyp+ixp+iz];
}
}
}
return 0;
}
/* copy communicated halo areas back to compute grids */
int newHaloVxVz(float *vx, float *vz, int npx, int npy, int npz, int halo, float *leftRecv, float *rightRecv, float *frontRecv, float *backRecv, float *topRecv, float *bottomRecv)
{
int ix, iy, iz, ih, iyp, ixp, halosizex, halosizey, halosizez, nxz;
nxz = npx*npz;
/* copy halo areas at left side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=0; ix<halo; ix++) {
ixp=ix*npz;
ih=ixp;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
vx[iyp+ixp+iz] = leftRecv[iy*npz*halo+ih+iz];
vz[iyp+ixp+iz] = leftRecv[halosizex+iy*npz*halo+ih+iz];
}
}
}
/* copy halo areas at right side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=npx-halo; ix<npx; ix++) {
ixp=ix*npz;
ih=(ix-(npx-halo))*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
vx[iyp+ixp+iz] = rightRecv[iy*npz*halo+ih+iz];
vz[iyp+ixp+iz] = rightRecv[halosizex+iy*npz*halo+ih+iz];
}
}
}
/* copy halo areas at front side */
halosizey = npx*npz*halo;
for (iy=0; iy<halo; iy++) {
iyp=iy*nxz;
ih=iyp;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
vx[iyp+ixp+iz] = frontRecv[ih+ixp+iz];
vz[iyp+ixp+iz] = frontRecv[halosizey+ih+ixp+iz];
}
}
}
/* copy halo areas at back side */
for (iy=npy-halo; iy<npy; iy++) {
iyp=iy*nxz;
ih=(iy-(npy-halo))*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
vx[iyp+ixp+iz] = backRecv[ih+ixp+iz];
vz[iyp+ixp+iz] = backRecv[halosizey+ih+ixp+iz];
}
}
}
/* copy halo areas at top side */
halosizez = npy*npx*halo;
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=0; iz<halo; iz++) {
ih=iz;
vx[iyp+ixp+iz] = topRecv[iy*npx*halo+ix*halo+ih];
vz[iyp+ixp+iz] = topRecv[halosizez+iy*npx*halo+ix*halo+ih];
}
}
}
/* copy halo areas at bottom side */
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=npz-halo; iz<npz; iz++) {
ih=(iz-(npz-halo));
vx[iyp+ixp+iz] = bottomRecv[iy*npx*halo+ix*halo+ih];
vz[iyp+ixp+iz] = bottomRecv[halosizez+iy*npx*halo+ix*halo+ih];
}
}
}
return 0;
}
/* copy communicated halo areas back to compute grids */
int newHaloP(float *p, int npx, int npy, int npz, int halo, float *leftRecv, float *rightRecv, float *frontRecv, float *backRecv, float *topRecv, float *bottomRecv)
{
int ix, iy, iz, ih, iyp, ixp, halosizex, halosizey, halosizez, nxz;
nxz = npx*npz;
/* copy halo areas at left side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=0; ix<halo; ix++) {
ixp=ix*npz;
ih=ixp;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] = leftRecv[iy*npz*halo+ih+iz];
}
}
}
/* copy halo areas at right side */
halosizex = npy*npz*halo;
for (iy=halo; iy<npy-halo; iy++) {
iyp=iy*nxz;
for (ix=npx-halo; ix<npx; ix++) {
ixp=ix*npz;
ih=(ix-(npx-halo))*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] = rightRecv[iy*npz*halo+ih+iz];
}
}
}
/* copy halo areas at front side */
halosizey = npx*npz*halo;
for (iy=0; iy<halo; iy++) {
iyp=iy*nxz;
ih=iyp;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] = frontRecv[ih+ixp+iz];
}
}
}
/* copy halo areas at back side */
for (iy=npy-halo; iy<npy; iy++) {
iyp=iy*nxz;
ih=(iy-(npy-halo))*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=halo; iz<npz-halo; iz++) {
p[iyp+ixp+iz] = backRecv[ih+ixp+iz];
}
}
}
/* copy halo areas at top side */
halosizez = npy*npx*halo;
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=0; iz<halo; iz++) {
ih=iz;
p[iyp+ixp+iz] = topRecv[iy*npx*halo+ix*halo+ih];
}
}
}
/* copy halo areas at bottom side */
for (iy=0; iy<npy; iy++) {
iyp=iy*nxz;
for (ix=0; ix<npx; ix++) {
ixp=ix*npz;
#pragma ivdep
for (iz=npz-halo; iz<npz; iz++) {
ih=(iz-(npz-halo));
p[iyp+ixp+iz] = bottomRecv[iy*npx*halo+ix*halo+ih];
}
}
}
return 0;
}
float gauss2time(float t, float f, float t0)
{
float value, time;
time = t-t0;
value = ((1.0-2.0*M_PI*M_PI*f*f*time*time)*exp(-M_PI*M_PI*f*f*time*time));
return value;
}
fdelmodc3D/CMWC4096.c 0 → 100644
+ 50
− 0
View file @ 8133dee1
/* http://en.wikipedia.org/wiki/Multiply-with-carry */
#include<stdlib.h>
#include<limits.h>
/* random number generator which can be used as an alternative for drand48() */
/* http://school.anhb.uwa.edu.au/personalpages/kwessen/shared/Marsaglia03.html*/
static unsigned long Q[4096],c=362436; /* choose random initial c<809430660 and */
/* 4096 random 32-bit integers for Q[] */
void seedCMWC4096(void)
{
int i;
for (i=0; i<4096; i++) {
Q[i] = lrand48();
}
return;
}
unsigned long CMWC4096(void)
{
unsigned long long t, a=18782LL;
static unsigned long i=4095;
unsigned long x,r=0xfffffffe;
i=(i+1)&4095;
t=a*Q[i]+c;
c=(t>>32);
x=t+c;
if(x<c){x++;c++;}
return(Q[i]=r-x);
}
/* replace defaults with five random seed values in calling program */
static unsigned long x=123456789,y=362436069,z=521288629,w=88675123,v=886756453;
unsigned long xorshift(void)
{unsigned long t;
t=(x^(x>>7)); x=y; y=z; z=w; w=v;
v=(v^(v<<6))^(t^(t<<13)); return (y+y+1)*v;}
double dcmwc4096(void)
{
double rd;
// rd = ((double)xorshift())/((double)ULONG_MAX);
rd = ((double)CMWC4096())/((double)ULONG_MAX);
return rd;
}
fdelmodc3D/FiguresPaper/Figure10.scr 0 → 100755
+ 162
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1
#PBS -N InterfModeling
#PBS -q fourweeks
#PBS -V
#
#reference and 2 SI results for visco-acoustic media, 2x200 s. + 2x1 hours
export PATH=../../bin:$PATH
makewave w=g1 fmax=30 t0=0.10 dt=0.0008 nt=4096 db=-40 file_out=G1.su verbose=1
makemod sizex=10000 sizez=4100 dx=10 dz=10 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=400,400 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 z=1100,1100,1100,1600,1100,1100,1100 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=2100,2100 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=2600,2600 cp=5500 ro=2200
xsrc1=1000
xsrc2=9000
zsrc1=3300
zsrc2=3900
dxsrc=10
base=fw5000_Q
../fdelmodc \
file_cp=simple_cp.su ischeme=2 \
Qp=15 \
file_den=simple_ro.su \
file_rcv=${base}.su \
file_src=G1.su \
rec_delay=0.1 \
dtrcv=0.008 \
verbose=3 \
tmod=4.108 \
dxrcv=50.0 \
src_random=0 \
wav_random=0 \
fmax=30 \
xsrc=5000 \
zsrc=0.0 \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
makemod sizex=10000 sizez=5000 dx=10 dz=10 cp0=1500 ro0=1000 file_base=hom.su
../fdelmodc \
file_cp=hom_cp.su ischeme=2 \
Qp=15 \
file_den=hom_ro.su \
file_rcv=${base}_D.su \
file_src=G1.su \
rec_delay=0.1 \
dtrcv=0.008 \
verbose=3 \
tmod=4.108 \
dxrcv=50.0 \
src_random=0 \
wav_random=0 \
fmax=30 \
xsrc=5000 \
zsrc=0.0 \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
suop2 ${base}_rvz.su ${base}_D_rvz.su > Reffw5000_Q_rvz.su
suwind s=1 j=1 tmax=4 f1=0.0 < Reffw5000_Q_rvz.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
supsimage perc=99 f1=0 f2=-5000 hbox=4 wbox=3 \
label1='time (s)' label2='lateral position (m)' \
labelsize=10 f2num=-5000 d2num=2500 > shotRefQ_5000_0.eps
tmod=120
tsrc1=0.1
tsrc2=120
tlength=120
nsrc=8000
fmax=30
file_shot=shotRQ_volume_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
echo $file_shot
#volume
zsrc1=500
zsrc2=4090
../fdelmodc \
file_cp=simple_cp.su ischeme=2 \
Qp=15 \
file_den=simple_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=1 \
tmod=$tmod \
dxrcv=10.0 \
plane_wave=0 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su 10
#deep
zsrc1=2700
zsrc2=4090
file_shot=shotRQ_deep_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
../fdelmodc \
file_cp=simple_cp.su ischeme=2 \
file_den=simple_ro.su \
Qp=15 \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=1 \
tmod=$tmod \
dxrcv=10.0 \
plane_wave=0 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su 10
fdelmodc3D/FiguresPaper/Figure10Hom.scr 0 → 100755
+ 141
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -q long
#PBS -V
#
export PATH=../../bin:$PATH
makewave w=g2 fp=10 t0=0.15 dt=0.0010 nt=4096 file_out=G2.su verbose=1
makemod sizex=10000 sizez=4100 dx=10 dz=10 cp0=1500 ro0=1000 file_base=hom.su
xsrc1=100
xsrc2=9900
dxsrc=10
tmod=120
tsrc1=0.1
tsrc2=120
tlength=120
nsrc=8000
fmax=30
for wav_random in 0 1;
do
file_shot=shotHRS${wav_random}_volume_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
echo $file_shot
#volume
zsrc1=500
zsrc2=4090
../fdelmodc \
file_cp=hom_cp.su ischeme=1 \
file_den=hom_ro.su \
file_src=G2.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=1 \
tmod=$tmod \
dxrcv=10.0 \
plane_wave=0 \
src_random=1 \
wav_random=${wav_random} \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su 10
#deep
zsrc1=2700
zsrc2=4090
file_shot=shotHRS${wav_random}_deep_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
fdelmodc \
file_cp=hom_cp.su ischeme=1 \
file_den=hom_ro.su \
file_src=G2.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=1 \
tmod=$tmod \
dxrcv=10.0 \
plane_wave=0 \
src_random=1 \
wav_random=${wav_random} \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su 10
#plane
zsrc1=2700
zsrc2=2700
file_shot=shotHRS${wav_random}_plane_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
fdelmodc \
file_cp=hom_cp.su ischeme=1 \
file_den=hom_ro.su \
file_src=G2.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=1 \
tmod=$tmod \
dxrcv=10.0 \
plane_wave=1 \
xsrc=5000 zsrc=2700 \
src_random=0 \
wav_random=${wav_random} \
fmax=$fmax \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su 10
done
fdelmodc3D/FiguresPaper/Figure11.scr 0 → 100755
+ 55
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -q long
#PBS -V
#
# calls fdelmodc_long.scr, can not be reproduced; software in test phase
echo " This Figure can not be reproduced completely"
echo " The progrom corrsmp used is still in test phase"
echo " Mail j.w.thorbecke@tudelft.nl if you want to reproduce this Figure."
exit;
./fdelmodc_long.scr
export OMP_NUM_THREADS=1
~/src/CorrSMP_FD_Files/corrsmp \
file_base=long/shotAcoustic_T3600_S1500_Dt500_F30_001_rvz.su \
nc=1 nstation=201 fullcorr=1 dt=0.008 \
nt=16384 verbose=1 ntc=4096 file_out=outI.su causal=1
#results long recordings and using corrsmp
suwind s=1 j=1 tmax=4 f1=0.0 key=fldr min=101 max=101 < outI_cc1.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
supsimage perc=99 f1=0 f2=-5000 hbox=4 wbox=3 \
label1='time (s)' label2='lateral position (m)' \
labelsize=10 f2num=-5000 d2num=2500 > long_Corr_I.eps
~/src/CorrSMP_FD_Files/corrsmp \
file_base=long/shotAcoustic_T3600_S1500_Dt500_F30_001_rvz.su \
nc=1 nstation=201 fullcorr=1 dt=0.008 \
nt=16384 verbose=1 ntc=4096 file_out=outII causal=2
suwind s=1 j=1 tmax=4 f1=0.0 key=fldr min=101 max=101 < outII_cc1.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
supsimage perc=99 f1=0 f2=-5000 hbox=4 wbox=3 \
label1='time (s)' label2='lateral position (m)' \
labelsize=10 f2num=-5000 d2num=2500 > long_Corr_II.eps
~/src/CorrSMP_FD_Files/corrsmp \
file_base=long/shotAcoustic_T3600_S1500_Dt500_F30_001_rvz.su \
nc=1 nstation=201 fullcorr=1 dt=0.008 \
nt=16384 verbose=1 ntc=4096 file_out=outIII causal=4
suwind s=1 j=1 tmax=4 f1=0.0 key=fldr min=101 max=101 < outIII_cc1.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
supsimage perc=99 f1=0 f2=-5000 hbox=4 wbox=3 \
label1='time (s)' label2='lateral position (m)' \
labelsize=10 f2num=-5000 d2num=2500 > long_Corr_III.eps
fdelmodc3D/FiguresPaper/Figure12.scr 0 → 100755
+ 46
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -q long
#PBS -V
#
# calls fdelmodc_amplitude.scr, can not be reproduced; software in test phase
echo " This Figure can not be reproduced completely"
echo " The progrom corrsmp used is still in test phase"
echo " Mail j.w.thorbecke@tudelft.nl if you want to reproduce this Figure."
exit;
./fdelmodc_amplitude.scr
#amplitude distribution of sources
supsgraph wbox=3 hbox=4 style=seimic < src_ampl.su \
f1=-250 d1=5.05 linecolor=black labelsize=10 titlesize=10 \
label1=amplitude label2=occurence > amplitudeDistr.eps
export OMP_NUM_THREADS=1
~/src/CorrSMP_FD_Files/corrsmp \
file_base=long/shotAcousticA500_T3600_S1500_Dt500_F30_001_rvz.su \
nc=1 nstation=201 fullcorr=1 dt=0.008 \
nt=16384 verbose=1 ntc=4096 file_out=outA.su causal=1
suwind key=fldr min=101 max=101 < outA.su_cc1.su > outA_Station101_Comp1.su
suwind s=1 j=1 tmax=4 f1=0.0 < outA_Station101_Comp1.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
supsimage perc=99 f1=0 f2=-5000 hbox=4 wbox=3 \
label1='time (s)' label2='lateral position (m)' \
labelsize=10 f2num=-5000 d2num=2500 > longA_Corr_I.eps
suwind s=1 j=1 tmax=4 f1=0.0 < outA_Station101_Comp1.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
sunormalize norm=rms | \
supsimage perc=99 f1=0 f2=-5000 hbox=4 wbox=3 \
label1='time (s)' label2='lateral position (m)' \
labelsize=10 f2num=-5000 d2num=2500 > longA_norm_Corr_I.eps
fdelmodc3D/FiguresPaper/Figure13.scr 0 → 100755
+ 163
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -q long
#PBS -V
#
# amplitude variations on source strength, 3x1500 s.
export PATH=../../bin:$PATH
makemod sizex=10000 sizez=4100 dx=10 dz=10 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=400,400 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 z=1100,1100,1100,1600,1100,1100,1100 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=2100,2100 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=2600,2600 cp=5500 ro=2200
xsrc1=500
xsrc2=9500
zsrc1=500
zsrc2=4090
tmod=120
tsrc1=0.1
tsrc2=120
tlength=120
nsrc=150
fmax=30
#Gaussian amplitude distribution
A=10000
file_shot=shotRS_A${A}_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
fdelmodc \
file_cp=simple_cp.su ischeme=1 \
file_den=simple_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=4 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=$A \
distribution=1 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
#use clip=500
SIrand.scr ${base}_rvz.su 50
#amplitude distribution of sources
f1=`surange < src_ampl.su | grep f1 | awk '{print $2 }'`
d1=`surange < src_ampl.su | grep d1 | awk '{print $2 }'`
supsgraph wbox=3 hbox=1 style=normal < src_ampl.su \
f1=$f1 d1=$d1 d1num=10000 x1end=50000 d2num=1 linecolor=black labelsize=10 titlesize=10 \
label1=amplitude label2=occurence > amplitudeDistrGauss${A}.eps
#Flat amplitude distribution
A=50000
file_shot=shotRS_A${A}_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
fdelmodc \
file_cp=simple_cp.su ischeme=1 \
file_den=simple_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=4 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=$A \
distribution=0 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
#use clip=900
SIrand.scr ${base}_rvz.su 50
#amplitude distribution of sources
#amplitude distribution of sources
f1=$(echo "scale=3; -5*$A" | bc -l)
d1=`surange < src_ampl.su | grep d1 | awk '{print $2 }'`
supsgraph wbox=3 hbox=1 style=normal < src_ampl.su \
f1=$f1 d1=$d1 d1num=10000 x1end=50000 d2num=1 linecolor=black labelsize=10 titlesize=10 \
label1=amplitude label2=occurence > amplitudeDistrFlat${A}.eps
#No amplitude distribution
A=0
file_shot=shotRS_A${A}_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
fdelmodc \
file_cp=simple_cp.su ischeme=1 \
file_den=simple_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=3 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
#use clip=4e-6
SIrand.scr ${base}_rvz.su 50
suspike nt=1000 ntr=1 nspk=1 ix1=1 it1=500 | sugain scale=$nsrc | \
supsgraph wbox=3 hbox=1 style=normal \
f1=-50000 d1=100.1 d1num=10000 d2num=50 x1end=50000 linecolor=black labelsize=10 titlesize=10 \
label1=amplitude label2=occurence > amplitude${A}.eps
fdelmodc3D/FiguresPaper/Figure13Amp.scr 0 → 100755
+ 130
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -q long
#PBS -V
#
# computes only the amplitude distributions pictures, 5 s.
export PATH=../../bin:$PATH
makemod sizex=10000 sizez=4100 dx=10 dz=10 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=400,400 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 z=1100,1100,1100,1600,1100,1100,1100 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=2100,2100 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=2600,2600 cp=5500 ro=2200
xsrc1=500
xsrc2=9500
zsrc1=500
zsrc2=4090
tmod=0.2
tsrc1=0.1
tsrc2=0.2
tlength=0.2
nsrc=150
fmax=30
#Gaussian amplitude distribution
A=10000
file_shot=shotRS_A${A}_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
../fdelmodc \
file_cp=simple_cp.su ischeme=1 \
file_den=simple_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=4 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=$A \
distribution=1 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
#amplitude distribution of sources
f1=`surange < src_ampl.su | grep f1 | awk '{print $2 }'`
d1=`surange < src_ampl.su | grep d1 | awk '{print $2 }'`
supsgraph wbox=3 hbox=1 style=normal < src_ampl.su \
f1=$f1 d1=$d1 d1num=10000 x1end=50000 d2num=1 linecolor=black labelsize=10 titlesize=10 \
label1=amplitude label2=occurence > amplitudeDistrGauss${A}.eps
#Flat amplitude distribution
A=50000
file_shot=shotRS_A${A}_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
fdelmodc \
file_cp=simple_cp.su ischeme=1 \
file_den=simple_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0010 \
verbose=4 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=$A \
distribution=0 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
#amplitude distribution of sources
#amplitude distribution of sources
f1=$(echo "scale=3; -5*$A" | bc -l)
d1=`surange < src_ampl.su | grep d1 | awk '{print $2 }'`
supsgraph wbox=3 hbox=1 style=normal < src_ampl.su \
f1=$f1 d1=$d1 d1num=10000 x1end=50000 d2num=1 linecolor=black labelsize=10 titlesize=10 \
label1=amplitude label2=occurence > amplitudeDistrFlat${A}.eps
#No amplitude distribution
A=0
file_shot=shotRS_A${A}_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
suspike nt=1000 ntr=1 nspk=1 ix1=1 it1=500 | sugain scale=$nsrc | \
supsgraph wbox=3 hbox=1 style=normal \
f1=-50000 d1=100.1 d1num=10000 d2num=50 x1end=50000 linecolor=black labelsize=10 titlesize=10 \
label1=amplitude label2=occurence > amplitude${A}.eps
fdelmodc3D/FiguresPaper/Figure14-15.scr 0 → 100755
+ 172
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1
#PBS -N InterfModeling
#PBS -q fourweeks
#PBS -V
#
# receivers and source placed on model with topography, 161 hours
echo " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
echo " This Figure can not be reproduced completely"
echo " The progrom corrsmp used is still in test phase"
echo " Mail j.w.thorbecke@tudelft.nl if you want to reproduce this Figure."
echo " Only the reference response will be modeled"
export PATH=../../bin:$PATH
dt=0.0004
ntap=120
fmax=45
makemod sizex=10000 sizez=4100 dx=5 dz=5 cp0=0 ro0=1000 file_base=real2.su \
orig=0,-800 gradunit=0 \
intt=def poly=2 cp=2450 ro=1000 gradcp=14 grad=0 \
x=0,1000,1700,1800,2000,3000,4000,4500,6000,6800,7000,7500,8100,8800,10000 \
z=-100,-200,-250,-200,-200,-120,-300,-600,-650,-500,-350,-200,-200,-150,-200 \
intt=rough var=200,3.2,1 poly=2 x=0,3000,8000,10000 \
z=400,250,300,500 cp=4500,4200,4800,4500 ro=1400 gradcp=5 grad=0 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=1100,1100,1100,1600,1100,1100,1100 cp=4000 ro=2000 gradcp=8 grad=0 \
intt=def poly=0 x=0,10000 z=1750,2050 cp=4500,5100 ro=1500 gradcp=13 grad=0 \
intt=def poly=0 x=0,10000 z=1850,2150 cp=6000,4200 ro=1500 gradcp=14 grad=0 \
intt=def poly=0 x=0,10000 z=1950,2250 cp=4800,4800 ro=1500 gradcp=5 grad=0 \
intt=def poly=0 x=0,10000 z=2000,2300 cp=6100,5000 ro=1500 gradcp=13 grad=0 \
intt=def poly=0 x=0,10000 z=2100,2400 cp=3800,5000 ro=1500 gradcp=20 grad=0 \
intt=def poly=0 x=0,10000 z=2150,2450 cp=5000 ro=1500 gradcp=14 grad=0 \
intt=def poly=0 x=0,10000 z=2350,2650 cp=5800 ro=1500 gradcp=5 grad=0 \
intt=def poly=0 x=0,10000 z=2600,2600 cp=5500 ro=2200 gradcp=5 grad=0
sushw key=f1 a=0 < real2_cp.su | \
sushw key=f1 a=0 | \
supsimage hbox=6 wbox=8 labelsize=12 f2num=-5000 d2num=1000 \
wrgb=0,0,1.0 grgb=0,1.0,0 brgb=1.0,0,0 legend=1 lnice=1 lstyle=vertright \
bclip=7053.02 wclip=0 label1="depth [m]" label2="lateral position [m]" \
> model2_cp.eps
makewave w=g2 fmax=45 t0=0.10 dt=$dt nt=4096 db=-40 file_out=G2.su verbose=1
extendModel file_in=real2_ro.su nafter=$ntap nbefore=$ntap nabove=0 nbelow=$ntap > vel2_edge_ro.su
extendModel file_in=real2_cp.su nafter=$ntap nbefore=$ntap nabove=0 nbelow=$ntap > vel2_edge_cp.su
#reference
fdelmodc \
file_cp=vel2_edge_cp.su ischeme=1 \
file_den=vel2_edge_ro.su \
file_rcv=shot_real2_x5000_topo.su \
file_src=G2.su \
dtrcv=0.004 \
verbose=4 \
tmod=3.004 \
dxrcv=20.0 \
zrcv1=-800 \
zrcv2=-800 \
xrcv1=0 \
xrcv2=10000 \
sinkdepth=1 \
src_random=0 \
wav_random=0 \
xsrc1=5000 \
xsrc2=5000 \
zsrc1=-800 \
tsrc1=0.0 \
dipsrc=1 \
ntaper=$ntap \
left=4 right=4 top=1 bottom=4
sushw key=f1,delrt a=0.0,0.0 < shot_real2_x5000_topo_rvz.su | \
basop choice=1 shift=-0.1 | \
supsimage clip=2e-12 f1=0 f2=-5000 x1end=3.004 hbox=8 wbox=6 \
label1="time (s)" label2="lateral position (m)" \
labelsize=10 f2num=-5000 d2num=1000 d1num=0.5 > shot_real2_x5000_topo.eps
#clip minclip=0 < SrcRecPositions.su > nep.su
#addmul file_in1=nep.su file_in2=vel2_edge_cp.su a=8000 | \
exit;
tlength=50
nsrc=1500
xsrc1=100
xsrc2=9900
zsrc1=1750
zsrc2=2600
fmax=30
tmod=3600
tsrc1=0.1
tsrc2=3600
file_shot=shot_real2_T${tmod}_S${nsrc}_Dt${tlength}_F${fmax}.su
fdelmodc \
file_cp=vel2_edge_cp.su ischeme=1 \
file_den=vel2_edge_ro.su \
file_rcv=$file_shot \
dtrcv=0.004 \
dt=$dt \
verbose=4 \
rec_ntsam=15000 \
tmod=$tmod \
dxrcv=20.0 \
zrcv1=-800 \
zrcv2=-800 \
xrcv1=0 \
xrcv2=10000 \
sinkdepth=1 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
dipsrc=0 \
ntaper=$ntap \
left=4 right=4 top=1 bottom=4
suop2 SrcRecPositions.su vel2_edge_cp.su w1=8000 w2=1.0 op=sum | \
sugain nclip=0 | \
supsimage hbox=6 wbox=8 labelsize=10 \
f2=-5600 f2num=-5000 d2num=1000 f1num=-500 d1num=500 \
x2beg=-5000 x2end=5000 legend=1 \
bclip=8200 wclip=0 label1="depth (m)" label2="lateral position (m)" \
> real_sources_cp.eps
supsimage < vel2_edge_cp.su hbox=6 wbox=8 labelsize=10 \
f2=-5600 f2num=-5000 d2num=1000 f1num=-500 d1num=500 \
x2beg=-5000 x2end=5000 legend=1 bps=24 lstyle=vertright \
bclip=6500 label1="depth (m)" label2="lateral position (m)" \
brgb=0.1,0.1,0.1 grgb=0.9,0.9,0.9 wrgb=1.0,1.0,1.0 \
curve=SrcPositions1500.txt npair=1500 curvecolor=black curvedash=-1 curvewidth=3 \
> real_sources_gray_cp.eps
supsimage hbox=6 wbox=8 labelsize=12 < SrcRecPositions.su \
f2=-5600 f2num=50000 d2num=1000 f1num=5000 d1num=500 \
legend=1 lnice=1 \
label1=" " label2=" " bclip=1 wclip=0 \
> rec_pos_real2.eps
export OMP_NUM_THREADS=4
export MKL_SERIAL=yes
export DFTI_NUMBER_OF_USER_THREADS=1
$HOME/src/CorrSMP_FD_Files/corrsmp file_base=shot_real2_T3600_S1500_Dt50_F30_001_rvz.su ntc=2048 fmax=60 \
verbose=3 nstation=501 \
file_out=corrsmp_out.su nc=1 fullcorr=1 dt=0.004 nt=15000
suwind key=fldr min=251 max=251 < corrsmp_out_cc1.su | \
supsimage clip=5e-13 f1=0 f2=-5000 d2=20 x1end=3.004 hbox=8 wbox=6 \
label1="time (s)" label2="lateral position (m)" \
labelsize=10 f2num=-5000 d2num=1000 d1num=0.5 > corr_real2_x5000_topo.eps
fdelmodc3D/FiguresPaper/Figure17_19AppendixA.scr 0 → 100755
+ 92
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1
#PBS -N InterfModeling
#PBS -q fourweeks
#PBS -V
#
export PATH=../../bin:$PATH
which fdelmodc
makemod sizex=10000 sizez=4100 dx=10 dz=10 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=400,400 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=1100,1100,1100,1600,1100,1100,1100 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=2100,2100 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=2600,2600 cp=5500 ro=2200
makewave w=g2 fmax=45 t0=0.10 dt=0.001 nt=4096 db=-40 file_out=G2.su verbose=1
xsrc1=100
xsrc2=9900
zsrc1=2100
zsrc2=4000
file_shot=shotRandomPos${xsrc}_${zsrc1}.su
fdelmodc \
file_cp=simple_cp.su ischeme=1 \
file_den=simple_ro.su \
file_rcv=$file_shot \
dtrcv=0.008 \
dt=0.0010 \
verbose=4 \
tmod=10.000 \
dxrcv=20.0 \
zrcv1=10 \
zrcv2=10 \
xrcv1=0 \
xrcv2=10000 \
src_random=1 \
wav_random=1 \
fmax=30 \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=5.0 \
nsrc=20 \
dipsrc=0 \
ntaper=45 \
tsnap1=0.1 tsnap2=6.0 dtsnap=0.2 \
left=4 right=4 top=1 bottom=4 \
nxmax=2500 nzmax=1400 ntmax=10000
psgraph < srcTimeLengthN\=10001.bin n1=10001 \
labelsize=12 d1=0.001 style=normal linecolor=blue \
label1="start time (s)" label2="source duration (s)" \
d1num=1 d2num=1 wbox=8 hbox=4 x1end=5 > srcTimeLength.eps
supswigp < src_nwav.su \
labelsize=10 label1='time (s)' label2='source number' x1end=6 \
d2=1 d2num=1 hbox=4 wbox=6 fill=0 \
titlesize=-1 > src_nwav.eps
suwind < src_nwav.su key=tracl min=11 max=11 | \
supsgraph hbox=2 wbox=4 style=normal \
labelsize=10 label2='amplitude' label1='time (s)' \
titlesize=-1 d1num=1.0 > src11_wiggle.eps
suwind < src_nwav.su key=tracl min=11 max=11 | \
supsgraph hbox=2 wbox=4 style=normal \
labelsize=10 label2='amplitude' label1='time (s)' \
titlesize=-1 x1end=0.05 > src11_wiggle_zbeg.eps
suwind < src_nwav.su key=tracl min=11 max=11 | \
supsgraph hbox=2 wbox=4 style=normal \
labelsize=10 label2='amplitude' label1='time (s)' \
titlesize=-1 x1beg=3.60 x1end=3.65 > src11_wiggle_zend.eps
suwind < src_nwav.su key=tracl min=11 max=11 | \
sufft | suamp| supsgraph hbox=2 wbox=4 style=normal \
labelsize=10 label2='amplitude' label1='frequency (Hz)' \
titlesize=-1 x1end=100 d1num=10 > src11_ampl.eps
fconv file_in1=src_nwav.su auto=1 shift=1 mode=cor1 | \
sugain qbal=1 | \
supswigp x1beg=-1 x1end=1 d2num=1 hbox=4 wbox=6 \
labelsize=10 label2='source number' label1='time (s)' \
titlesize=-1 fill=0 > src_nwav_autoCorr_Norm.eps
fdelmodc3D/FiguresPaper/Figure2.scr 0 → 100755
+ 104
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
# starts fdelmodc only to compute the source positions, 1 s.
export PATH=../../bin:$PATH
makemod sizex=10000 sizez=5000 dx=5 dz=5 cp0=0 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=900,900 cp=1500 ro=1000 \
intt=def poly=0 x=0,10000 z=1300,1300 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=2000,2000,2000,2500,2000,2000,2000 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=3000,3000 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=3500,3500 cp=5500 ro=2200
suwind itmin=181 < simple_cp.su | sushw key=f1 a=0 > vel_cp.su
suwind itmin=181 < simple_ro.su | sushw key=f1 a=0 > vel_ro.su
xsrc1=100
xsrc2=9900
dxsrc=10
#volume
zsrc1=500
zsrc2=4090
tmod=0.01
tsrc2=0.1
tlength=0.01
nsrc=1000
fmax=30
#Figure 2b,c,d,e,f
file_shot=shotR_T${tmod}_S${nsrc}_Dt${tlength}_F${fmax}.su
echo $file_shot
#dummy modeling just to generate the source positions to be used in the Figure
../fdelmodc \
file_cp=vel_cp.su ischeme=1 \
file_den=vel_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0005 \
verbose=4 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=0 \
xsrc=5000 zsrc=2700 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
ntaper=45 \
left=4 right=4 top=1 bottom=4
#Figure 2
suop2 SrcRecPositions.su vel_cp.su w1=7000 w2=1.0 op=sum | \
sugain nclip=1500 | \
sushw key=f1,f2 a=0,-5000 | \
supsimage hbox=6 wbox=8 labelsize=10 \
f2num=-5000 d2num=1000 f1num=-500 d1num=500 \
x2beg=-5000 x2end=5000 x1beg=-900 x1end=4000 \
bclip=7000 wclip=0 label1="depth (m)" label2="lateral position (m)" \
x2beg=-5000 x2end=5000 \
> simple_sources_cp.eps
# laternative color scheme for above picture
# wrgb=0,0,1.0 grgb=0,1.0,0 brgb=1.0,0,0 x2beg=-5000 x2end=5000 \
supsimage hbox=6 wbox=8 labelsize=18 < vel_cp.su \
f2=-5000 f2num=-5000 d2num=1000 f1num=-500 d1num=500 \
x2beg=-5000 x2end=5000 x1beg=-900 x1end=4000 \
label1="depth (m)" label2="lateral position (m)" \
wrgb=0,0,1.0 grgb=0,1.0,0 brgb=1.0,0,0 x2beg=-5000 x2end=5000 \
> simple_cp.eps
# use adapted psimage from su to plot points using curve function
supsimage hbox=6 wbox=8 labelsize=10 < simple_cp.su \
f1=-900 f2=-5000 f2num=-5000 d2num=1000 f1num=-500 d1num=500 \
x2beg=-5000 x2end=5000 x1beg=-900 x1end=4000 \
label1="depth (m)" label2="lateral position (m)" \
wrgb=0,0,1.0 grgb=0,1.0,0 brgb=1.0,0,0 x2beg=-5000 x2end=5000 \
curve=SrcPositions1000.txt npair=1000 curvecolor=black curvedash=-1 curvewidth=1 \
> simple_srcpos_cp.eps
supsimage hbox=6 wbox=8 labelsize=18 < SrcRecPositions.su\
f2num=50000 d2num=1000 f1num=50000 d1num=500 \
bclip=1 wclip=0 label1=" " label2=" "\
> sources.eps
fdelmodc3D/FiguresPaper/Figure20AppendixA.scr 0 → 100755
+ 148
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
cp=2000
rho=1000
dx=2.5
dt=0.0001
halfdt=`perl -e "print 0.5*$dt;"`
export PATH=../../bin:$PATH
makemod sizex=2000 sizez=2000 dx=$dx dz=$dx cp0=$cp ro0=$rho orig=-1000,0 file_base=simple.su
makewave fp=15 dt=$dt file_out=wave.su nt=4096 t0=0.1
######### MONOPOLE #######
../fdelmodc \
file_cp=simple_cp.su ischeme=1 iorder=4 \
file_den=simple_ro.su \
file_src=wave.su \
file_rcv=shot_fd.su \
src_type=1 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.0005 \
verbose=2 \
tmod=0.5115 \
dxrcv=5.0 \
xrcv1=-500 xrcv2=500 \
zrcv1=500 zrcv2=500 \
xsrc=0 zsrc=1000 \
ntaper=80 \
left=4 right=4 top=4 bottom=4
makewave fp=15 dt=0.0005 file_out=wave.su nt=4096 t0=0.1
green c=$cp rho=$rho file_src=wave.su zsrc1=500 xrcv=-500,500 dxrcv=5 nt=4096 dip=0 | suwind nt=1024 > shot_green_rp.su
green c=$cp rho=$rho file_src=wave.su zsrc1=500 xrcv=-500,500 dxrcv=5 nt=4096 dip=0 p_vz=1 | suwind nt=1024 > shot_green_rvz.su
# rp
(suwind key=tracl min=101 max=101 < shot_fd_rp.su | basop choice=shift shift=$halfdt; suwind key=tracl min=101 max=101 < shot_green_rp.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=4 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 d2num=0.04 x2end=0.080 f2num=-0.4 > mon_rp.eps
(suwind key=tracl min=101 max=101 < shot_fd_rp.su | basop choice=shift shift=$halfdt; suwind key=tracl min=101 max=101 < shot_green_rp.su;) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=2 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green x1beg=0.255 x1end=0.258 x2beg=0.0785 x2end=0.0797 d2num=0.0002 f2num=0.0786 > mon_zoom_rp.eps
suwind key=tracl min=101 max=101 < shot_fd_rp.su | basop choice=shift shift=$halfdt > trace_fd.su
suwind key=tracl min=101 max=101 < shot_green_rp.su > trace_green.su
sumax < trace_green.su outpar=nep
Mmax=`cat nep | awk '{print $1}'`
a=`perl -e "print 100.0/$Mmax;"`
echo $a
sudiff trace_green.su trace_fd.su | basop choice=shift shift=-0.1 | sugain scale=$a | suop op=abs | supsgraph style=normal labelsize=12 wbox=4 hbox=2 label1='time in seconds' label2="Relative error in percentage of peak" linewidth=0.1 linecolor=red x2beg=0.0 x2end=1.0 f2num=0.0 d2num=0.5 > mon_diff_dx${dx}_rp.eps
(suwind key=tracl min=101 max=101 < shot_fd_rp.su | basop choice=shift shift=$halfdt; suwind key=tracl min=101 max=101 < shot_green_rp.su ) | basop choice=shift shift=-0.1 | suxgraph
# rvz
(suwind key=tracl min=101 max=101 < shot_fd_rvz.su ; suwind key=tracl min=101 max=101 < shot_green_rvz.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=4 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 d2num=2e-8 f2num=-4e-8 x2end=4e-8 x2beg=-2.7e-8 > mon_rvz.eps
(suwind key=tracl min=101 max=101 < shot_fd_rvz.su ; suwind key=tracl min=101 max=101 < shot_green_rvz.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=2 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green d2num=2e-10 f2num=3.86e-8 x2end=3.96e-8 x2beg=3.86e-8 x1beg=0.255 x1end=0.258 > mon_zoom_rvz.eps
suwind key=tracl min=101 max=101 < shot_fd_rvz.su > trace_fd.su
suwind key=tracl min=101 max=101 < shot_green_rvz.su > trace_green.su
sumax < trace_green.su outpar=nep
Mmax=`cat nep | awk '{print $1}'`
a=`perl -e "print 100.0/$Mmax;"`
echo $a
sudiff trace_green.su trace_fd.su | basop choice=shift shift=-0.1 | sugain scale=$a | suop op=abs | supsgraph style=normal labelsize=12 wbox=4 hbox=2 label1='time in seconds' label2="Relative error in percentage of peak" linewidth=0.1 linecolor=red x2beg=0.0 x2end=1.0 f2num=0.0 d2num=0.5 > mon_diff_dx${dx}_rvz.eps
(suwind key=tracl min=101 max=101 < shot_fd_rvz.su ; suwind key=tracl min=101 max=101 < shot_green_rvz.su) | basop choice=shift shift=-0.1 | suxgraph
exit;
######### DIPOLE #######
makewave fp=15 dt=$dt file_out=wave.su nt=4096 t0=0.1
fdelmodc \
file_cp=simple_cp.su ischeme=1 iorder=4 \
file_den=simple_ro.su \
file_src=wave.su \
file_rcv=shot_fd_dip.su \
src_type=1 \
src_orient=2 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.0005 \
verbose=2 \
tmod=0.5115 \
dxrcv=5.0 \
xrcv1=-500 xrcv2=500 \
zrcv1=500 zrcv2=500 \
xsrc=0 zsrc=1000 \
ntaper=80 \
left=4 right=4 top=4 bottom=4
makewave fp=15 dt=0.0005 file_out=wave.su nt=4096 t0=0.1
green c=$cp rho=$rho file_src=wave.su zsrc1=500 xrcv=-500,500 dxrcv=5 nt=4096 dip=1 | suwind nt=1024 > shot_green_dip_rp.su
green c=$cp rho=$rho file_src=wave.su zsrc1=500 xrcv=-500,500 dxrcv=5 nt=4096 p_vz=1 dip=1 | suwind nt=1024 > shot_green_dip_rvz.su
shift=`perl -e "print 0.5*$dx/$cp;"`
# rp
(suwind key=tracl min=101 max=101 < shot_fd_dip_rp.su | basop choice=shift shift=-$shift; suwind key=tracl min=101 max=101 < shot_green_dip_rp.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal titlesize=-1 labelsize=10 wbox=4 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green d2num=0.001 x2end=0.004001 f2num=-0.03 > dip_rp.eps
(suwind key=tracl min=101 max=101 < shot_fd_dip_rp.su | basop choice=shift shift=-$shift; suwind key=tracl min=101 max=101 < shot_green_dip_rp.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=2 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green x2beg=0.00379 x2end=0.00388 d2num=0.00002 f2num=0.00380 x1beg=0.2430 x1end=0.24601 f1num=0.243 d1num=0.002 > dip_zoom_rp.eps
suwind key=tracl min=101 max=101 < shot_fd_dip_rp.su | basop choice=shift shift=-$shift > trace_fd.su
suwind key=tracl min=101 max=101 < shot_green_dip_rp.su > trace_green.su
sumax < trace_green.su outpar=nep
Mmax=`cat nep | awk '{print $1}'`
a=`perl -e "print 100.0/$Mmax;"`
echo $a
sudiff trace_green.su trace_fd.su | basop choice=shift shift=-0.1 | sugain scale=$a | suop op=abs | supsgraph style=normal labelsize=12 wbox=4 hbox=2 label1='time in seconds' label2="Relative error in percentage of peak" linewidth=0.1 linecolor=red x2beg=0.0 x2end=1.0 f2num=0.0 d2num=0.5 > dip_diff_dx${dx}_rp.eps
(suwind key=tracl min=101 max=101 < shot_fd_dip_rp.su | basop choice=shift shift=-$shift ; suwind key=tracl min=101 max=101 < shot_green_dip_rp.su ) | basop choice=shift shift=-0.1 | suxgraph
# rvz
(suwind key=tracl min=101 max=101 < shot_fd_dip_rvz.su | basop choice=shift shift=-$shift; suwind key=tracl min=101 max=101 < shot_green_dip_rvz.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=4 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green d2num=6e-10 f2num=-1.2e-9 x2end=2e-9 x2beg=-1.3e-9 > dip_rvz.eps
(suwind key=tracl min=101 max=101 < shot_fd_dip_rvz.su | basop choice=shift shift=-$shift; suwind key=tracl min=101 max=101 < shot_green_dip_rvz.su ) | basop choice=shift shift=-0.1 | supsgraph style=normal labelsize=10 wbox=2 hbox=2 label1='time in seconds' label2="Amplitude" linewidth=0.1 linecolor=red,green d2num=2e-11 f2num=1.85e-9 x2end=1.944e-9 x2beg=1.85e-9 x1beg=0.2430 x1end=0.24601 f1num=0.243 d1num=0.002 > dip_zoom_rvz.eps
suwind key=tracl min=101 max=101 < shot_fd_dip_rvz.su | basop choice=shift shift=-$shift > trace_fd.su
suwind key=tracl min=101 max=101 < shot_green_dip_rvz.su > trace_green.su
sumax < trace_green.su outpar=nep
Mmax=`cat nep | awk '{print $1}'`
a=`perl -e "print 100.0/$Mmax;"`
echo $a
sudiff trace_green.su trace_fd.su | basop choice=shift shift=-0.1 | sugain scale=$a | suop op=abs | supsgraph style=normal labelsize=12 wbox=4 hbox=2 label1='time in seconds' label2="Relative error in percentage of peak" linewidth=0.1 linecolor=red x2beg=0.0 x2end=1.0 f2num=0.0 d2num=0.5 > dip_diff_dx${dx}_rvz.eps
(suwind key=tracl min=101 max=101 < shot_fd_dip_rvz.su | basop choice=shift shift=-$shift; suwind key=tracl min=101 max=101 < shot_green_dip_rvz.su ) | basop choice=shift shift=-0.1 | suxgraph
fdelmodc3D/FiguresPaper/Figure3.scr 0 → 100755
+ 106
− 0
View file @ 8133dee1
#!/bin/bash
#
#Forward model all source positions one by one (this takes a very long time)
#
# calls Simple_model_base, and Simple_model_sides.scr, 122 hours!
export PATH=../../bin:$PATH
./Simple_model_base.scr
./Simple_model_sides.scr
# Correlation horizontal layer
export dxsrc=10
export xsrc1=1000
export xsrc2=9000
export xrcv1=0
export xrcv2=10000
export file_out=/tmp/corr.su
export file_in=T_simple_dxrcv50_dx10.su
export dxrcv=50
(( nshots = ($xsrc2-$xsrc1)/$dxsrc + 1 ))
(( nrecv = ($xrcv2-$xrcv1)/$dxrcv + 1 ))
(( middle = $xsrc1 + (($nshots-1)/2)*$dxsrc ))
(( xa = ($nrecv+1)/2 ))
(( xsrc_sim = $xrcv1 + ($xa-1)*$dxrcv ))
echo xsrc_sim=$xsrc_sim
echo nshot=$nshots
echo nrec=$nrecv
echo xa=$xa
echo middle=$middle
# first do the correlation with the middle trace of each shot
rm -rf $file_out
for (( xsrc = $xsrc1; xsrc<=$xsrc2; xsrc+=$dxsrc )) do
echo $xsrc $middle;
sxkey=${xsrc}000
suwind key=sx min=$sxkey max=$sxkey < $file_in > shot.su
suwind key=tracl min=$xa max=$xa < shot.su | \
fconv mode=cor2 file_in1=shot.su shift=1 >> $file_out
done;
#then sort on constant gx values and stack all traces belonging to the same gx (=sum over sx)
(( f2 = -($middle-$xrcv1) ))
susort < $file_out gx > /tmp/corr_sort.su
sustack key=gx < /tmp/corr_sort.su | sushw key=f2 a=$f2 > corr_stack${dxrcv}_${xsrc_sim}_${dxsrc}.su
#Correlation sides
export file_in=T_simple_sides_dxrcv50_dx10.su
nshots=720
export fldr1=902
export fldr2=1621
echo $xsrc_sim
echo $nshots
echo $nrecv
echo $xa
echo $middle
rm -rf $file_out
for (( fldr = $fldr1; fldr<=$fldr2; fldr+=1 )) do
echo $fldr;
suwind key=fldr min=$fldr max=$fldr < $file_in > shot.su
suwind key=tracl min=$xa max=$xa < shot.su | \
fconv mode=cor2 file_in1=shot.su shift=1 >> $file_out
done;
(( f2 = -($middle-$xrcv1) ))
susort < $file_out gx > corr_sort.su
sustack key=gx < corr_sort.su | sushw key=f2 a=$f2 > corr_stack${dxrcv}_sides_$xsrc_sim.su
rm $file_out /tmp/corr_sort.su shot.su
# sources at left and right sides of the model receivers on top gray scales dxrc
cat corr_stack50_sides_5000.su | \
suwind tmax=4 tmin=0 f1=-8.008 | \
sushw key=f1,delrt a=0.0,0.0 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
supsimage clip=2e-15 f1=0 f2=-5000 hbox=4 wbox=3 \
labelsize=10 f2num=-5000 d2num=2500 \
label1='time (s)' label2='lateral position (m)' > SimCorr50_sides_sx5000_simple.eps
cat corr_stack50_5000_10.su | \
suwind tmax=4 tmin=0 f1=-8.008 | \
sushw key=f1,delrt a=0.0,0.0 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
supsimage clip=2e-15 f1=0 f2=-5000 hbox=4 wbox=3 \
labelsize=10 f2num=-5000 d2num=2500 \
label1='time (s)' label2='lateral position (m)'> SimCorr50_sx5000_simple.eps
susum corr_stack50_5000_10.su corr_stack50_sides_5000.su | \
suwind tmax=4 tmin=0 f1=-8.008 | \
sushw key=f1,delrt a=0.0,0.0 | \
sufilter amps=0,0.5,1,1,0 f=0,2,3,50,60 | \
supsimage clip=2e-15 f1=0 f2=-5000 hbox=4 wbox=3 \
labelsize=10 f2num=-5000 d2num=2500 \
label1='time (s)' label2='lateral position (m)'> SimCorr50_add_sides_sx5000_simple.eps
fdelmodc3D/FiguresPaper/Figure3_nofree.scr 0 → 100755
+ 58
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
# direct modeled reference result Figure 3d, 500 s.
export PATH=../../bin:$PATH
makewave w=g1 fmax=45 t0=0.15 dt=0.0005 nt=4096 db=-40 file_out=G1.su verbose=1
cp G1.su G1_c.su
fconv file_in1=G1.su file_in2=G1_c.su mode=cor1 verbose=1 > corrw.su
basop choice=shift shift=0.2 dx=1 file_in=corrw.su > G1corr.su
makemod sizex=10000 sizez=5000 dx=5 dz=5 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=1300,1300 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=2000,2000,2000,2500,2000,2000,2000 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=3000,3000 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=3500,3500 cp=5500 ro=2200
sushw key=f1 a=-900 < simple_cp.su > vel_cp.su
sushw key=f1 a=-900 < simple_ro.su > vel_ro.su
xsrc=5000
file_shot=shotRefNoFree_x${xsrc}.su
../fdelmodc \
file_cp=vel_cp.su ischeme=1 \
file_den=vel_ro.su \
file_src=G1corr.su \
fmax=45 \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
verbose=4 \
rec_delay=0.2 \
tmod=4.2 \
dxrcv=50.0 \
zrcv1=0 \
zrcv2=0 \
xrcv1=0 xrcv2=10000 \
plane_wave=0 \
amplitude=0 \
xsrc=${xsrc} zsrc=0 \
ntaper=45 \
left=4 right=4 top=4 bottom=4
suwind s=1 j=1 tmax=4.008 f1=0.0 < shotRefNofree_x${xsrc}_rvz.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
sugain scale=-1 | \
supsimage f1=0 f2=-5000 hbox=4 wbox=3 x1end=4.0 \
labelsize=10 f2num=-5000 d2num=2500 verbose=1 clip=2e-7 \
label1='time (s)' label2='lateral position (m)' > shotRefNofree_5000_0.eps
fdelmodc3D/FiguresPaper/Figure3_ref.scr 0 → 100755
+ 55
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
# direct modeled reference result Figure 3d, 500 s.
export PATH=../../bin:$PATH
makewave w=g1 fmax=45 t0=0.15 dt=0.0005 nt=4096 db=-40 file_out=G1.su verbose=1
cp G1.su G1_c.su
fconv file_in1=G1.su file_in2=G1_c.su mode=cor1 verbose=1 > corrw.su
basop choice=shift shift=0.2 dx=1 file_in=corrw.su > G1corr.su
makemod sizex=10000 sizez=5000 dx=5 dz=5 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=1300,1300 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=2000,2000,2000,2500,2000,2000,2000 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=3000,3000 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=3500,3500 cp=5500 ro=2200
suwind itmin=181 < simple_cp.su | sushw key=f1 a=0 > vel_cp.su
suwind itmin=181 < simple_ro.su | sushw key=f1 a=0 > vel_ro.su
xsrc=5000
file_shot=shotRef_x${xsrc}.su
../fdelmodc \
file_cp=vel_cp.su ischeme=1 \
file_den=vel_ro.su \
file_src=G1corr.su \
fmax=45 \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
verbose=4 \
rec_delay=0.2 \
tmod=0.5 \
dxrcv=50.0 \
plane_wave=0 \
amplitude=0 \
xsrc=${xsrc} zsrc=5 \
ntaper=45 \
left=4 right=4 top=1 bottom=4 verbose=1
suwind s=1 j=1 tmax=4.008 f1=0.0 < shotRef_x${xsrc}_rvz.su | \
sushw key=f1,delrt,d2 a=0.0,0.0,50 | \
sugain scale=-1 | \
supsimage f1=0 f2=-5000 hbox=4 wbox=3 x1end=4.0 \
labelsize=10 f2num=-5000 d2num=2500 verbose=1 clip=2e-7 \
label1='Time (s)' label2='Lateral position (m)' > shotRef_5000_0.eps
fdelmodc3D/FiguresPaper/Figure4.scr 0 → 100755
+ 75
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
# inside the length loop the script SIrand.scr is called to compute the retrieved response from the modeled data
# 5 different source signature lengths, 5x3.5 hours
export PATH=../../bin:$PATH
makemod sizex=10000 sizez=5000 dx=5 dz=5 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=1300,1300 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=2000,2000,2000,2500,2000,2000,2000 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=3000,3000 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=3500,3500 cp=5500 ro=2200
suwind itmin=181 < simple_cp.su | sushw key=f1 a=0 > vel_cp.su
suwind itmin=181 < simple_ro.su | sushw key=f1 a=0 > vel_ro.su
xsrc1=100
xsrc2=9900
dxsrc=10
#volume
zsrc1=500
zsrc2=4090
tmod=120
tsrc2=120
tlength=120
nsrc=1000
fmax=30
#Figure 7b,c,d,e,f
for tlength in 120 60 30 10 5;
do
file_shot=shotR_T${tmod}_S${nsrc}_Dt${tlength}_F${fmax}.su
echo $file_shot
fdelmodc \
file_cp=vel_cp.su ischeme=1 \
file_den=vel_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0005 \
verbose=2 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=0 \
xsrc=5000 zsrc=2700 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su 50
done
fdelmodc3D/FiguresPaper/Figure5.scr 0 → 100755
+ 66
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
# simulates 8000 short (2.5 s) sources, 3.5 hours
export PATH=../../bin:$PATH
makemod sizex=10000 sizez=5000 dx=5 dz=5 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=1300,1300 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=2000,2000,2000,2500,2000,2000,2000 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=3000,3000 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=3500,3500 cp=5500 ro=2200
suwind itmin=181 < simple_cp.su | sushw key=f1 a=0 > vel_cp.su
suwind itmin=181 < simple_ro.su | sushw key=f1 a=0 > vel_ro.su
xsrc1=100
xsrc2=9900
#volume
tmod=120
tsrc2=120
fmax=30
zsrc1=500
zsrc2=4090
tlength=5
nsrc=8000
file_shot=shotR_T${tmod}_S${nsrc}_Dt${tlength}_F${fmax}.su
fdelmodc \
file_cp=vel_cp.su ischeme=1 \
file_den=vel_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0005 \
verbose=2 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=0 \
xsrc=5000 zsrc=2700 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su
fdelmodc3D/FiguresPaper/Figure6.scr 0 → 100755
+ 102
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
# 5 different number of random sources, 5x3.5 hours
export PATH=../../bin:$PATH
#makewave w=g1 fp=10 t0=0.15 dt=0.0010 nt=4096 file_out=G2.su verbose=1
makemod sizex=10000 sizez=5000 dx=5 dz=5 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=1300,1300 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=2000,2000,2000,2500,2000,2000,2000 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=3000,3000 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=3500,3500 cp=5500 ro=2200
suwind itmin=181 < simple_cp.su | sushw key=f1 a=0 > vel_cp.su
suwind itmin=181 < simple_ro.su | sushw key=f1 a=0 > vel_ro.su
xsrc1=100
xsrc2=9900
#volume
zsrc1=500
zsrc2=4090
tmod=120
tsrc2=120
tlength=120
fmax=30
#Figure 6a,b,c,d,e
rm Trace.su
for nsrc in 8000 1000 500 100 50;
do
file_shot=shotR_T${tmod}_S${nsrc}_Dt${tsrc2}_F${fmax}.su
echo $file_shot
fdelmodc \
file_cp=vel_cp.su ischeme=1 \
file_den=vel_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0005 \
verbose=2 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=0 \
xsrc=5000 zsrc=2700 \
src_random=1 \
wav_random=1 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su
ntraces=`surange < ${base}_rvz.su | grep traces| awk '{print $1 }'`
echo $ntraces
middle=$(echo "scale=0; ($ntraces+1)/2"| bc -l)
echo $middle
susum causal.su noncausal.su | \
suwind s=1 j=1 tmax=4 f1=0.0 | \
sushw key=f1,delrt,d2,fldr a=0.0,0.0,50,8000 | \
suwind key=tracl min=$middle max=$middle >> Trace.su
done
#add reference trace to trace comparison
suwind key=tracf min=101 max=101 < shotRef_x5000_rvz.su | sugain scale=-1 >> Trace.su
suwind tmin=1.5 tmax=3.0 f1=0.0 < Trace.su | \
sunormalize norm=max | \
supsgraph \
hbox=4 wbox=3 labelsize=10 x1beg=1.5 x1end=3.0 \
linecolor=red,green,blue,emerald,black f1=1.5 d1num=0.5 \
label1='time (s)' > shotTraces_T120_S_Dt120_F30.eps
suwind tmin=1.5 tmax=3.1 f1=0.0 < Trace.su | \
sunormalize norm=max | \
supswigp \
hbox=4 wbox=3 labelsize=10 x1beg=1.5 x1end=3.0 \
f1=1.5 d1num=0.5 linewidth=1 fill=0 f2=1 d2=1 d2num=1 \
label1='time (s)' label2='number of sources' > shotWiggles_T120_S_Dt120_F30.eps
fdelmodc3D/FiguresPaper/Figure6f.scr 0 → 100755
+ 59
− 0
View file @ 8133dee1
#!/bin/bash
#
# make postscript file of middle trace after Figure6.scr, 1 s.
shot=shotR_T120_S8000_Dt120_F30_rvz.su
base=`echo $shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
ntraces=`surange < $shot | grep traces| awk '{print $1 }'`
echo $ntraces
middle=$(echo "scale=0; ($ntraces+1)/2"| bc -l)
echo $middle
./SIrand.scr $shot
susum causal.su noncausal.su | \
suwind s=1 j=1 tmax=4 f1=0.0 | \
sushw key=f1,delrt,d2,fldr a=0.0,0.0,50,8000 | \
suwind key=tracl min=$middle max=$middle > Trace.su
shot=shotR_T120_S1000_Dt120_F30_rvz.su
./SIrand.scr $shot
susum causal.su noncausal.su | \
suwind s=1 j=1 tmax=4 f1=0.0 | \
sushw key=f1,delrt,d2,fldr a=0.0,0.0,50,1000 | \
suwind key=tracl min=$middle max=$middle >> Trace.su
shot=shotR_T120_S500_Dt120_F30_rvz.su
./SIrand.scr $shot
susum causal.su noncausal.su | \
suwind s=1 j=1 tmax=4 f1=0.0 | \
sushw key=f1,delrt,d2,fldr a=0.0,0.0,50,500 | \
suwind key=tracl min=$middle max=$middle >> Trace.su
shot=shotR_T120_S100_Dt120_F30_rvz.su
./SIrand.scr $shot
susum causal.su noncausal.su | \
suwind s=1 j=1 tmax=4 f1=0.0 | \
sushw key=f1,delrt,d2,fldr a=0.0,0.0,50,100 | \
suwind key=tracl min=$middle max=$middle >> Trace.su
suwind s=1 j=1 tmax=4 f1=0.0 < shotRef_5000_0_rvz.su | \
sushw key=f1,delrt,d2,fldr a=0.0,0.0,50,100 | \
suwind key=tracl min=$middle max=$middle >> Trace.su
suwind tmin=1.5 tmax=3.0 f1=0.0 < Trace.su | \
sunormalize norm=max | \
supsgraph \
hbox=4 wbox=3 labelsize=10 x1beg=1.5 x1end=3.0 \
linecolor=red,green,blue,emerald,black f1=1.5 d1num=0.5 \
> shotTraces_T120_S_Dt120_F30.eps
suwind tmin=1.5 tmax=3.1 f1=0.0 < Trace.su | \
sunormalize norm=max | \
supswigp \
hbox=4 wbox=3 labelsize=10 x1beg=1.5 x1end=3.0 \
f1=1.5 d1num=0.5 linewidth=1 fill=0 f2=1 d2=1 \
label1='time (s)' label2='number of sources' > shotWiggles_T120_S_Dt120_F30.eps
fdelmodc3D/FiguresPaper/Figure6length.scr 0 → 100755
+ 100
− 0
View file @ 8133dee1
#!/bin/bash
#PBS -l nodes=1:ppn=2
#PBS -N InterfModeling
#PBS -V
#
# alternative not used in paper, fixed source signature length, 5x3.5 hours
export PATH=../../bin:$PATH
#makewave w=g1 fp=10 t0=0.15 dt=0.0010 nt=4096 file_out=G2.su verbose=1
makemod sizex=10000 sizez=5000 dx=5 dz=5 cp0=1500 ro0=1000 file_base=simple.su \
intt=def poly=0 x=0,10000 z=1300,1300 cp=2000 ro=1400 \
intt=def poly=2 x=0,2000,3000,5000,7000,8000,10000 \
z=2000,2000,2000,2500,2000,2000,2000 cp=4000 ro=2000 \
intt=def poly=0 x=0,10000 z=3000,3000 cp=3000 ro=1500 \
intt=def poly=0 x=0,10000 z=3500,3500 cp=5500 ro=2200
suwind itmin=181 < simple_cp.su | sushw key=f1 a=0 > vel_cp.su
suwind itmin=181 < simple_ro.su | sushw key=f1 a=0 > vel_ro.su
xsrc1=100
xsrc2=9900
#volume
zsrc1=500
zsrc2=4090
tmod=120
tsrc2=120
tlength=120
fmax=30
#Figure 6a,b,c,d,e
rm Trace.su
for nsrc in 8000 1000 500 100 50;
do
file_shot=shotRVt_T${tmod}_S${nsrc}_Dt${tlength}_F${fmax}.su
echo $file_shot
fdelmodc \
file_cp=vel_cp.su ischeme=1 \
file_den=vel_ro.su \
file_rcv=$file_shot \
rec_type_p=0 \
dtrcv=0.008 \
rec_ntsam=16384 \
dt=0.0005 \
verbose=4 \
tmod=$tmod \
dxrcv=50.0 \
plane_wave=0 \
amplitude=0 \
xsrc=5000 zsrc=2700 \
src_random=1 \
wav_random=1 \
length_random=0 \
fmax=$fmax \
xsrc1=$xsrc1 \
xsrc2=$xsrc2 \
zsrc1=$zsrc1 \
zsrc2=$zsrc2 \
tsrc1=0.0 \
tsrc2=$tsrc2 \
tlength=$tlength \
nsrc=$nsrc \
ntaper=45 \
left=4 right=4 top=1 bottom=4
base=`echo $file_shot | awk 'BEGIN { FS = "." } ; { print $1 }'`
echo $base
SIrand.scr ${base}_rvz.su
ntraces=`surange < ${base}_rvz.su | grep traces| awk '{print $1 }'`
echo $ntraces
middle=$(echo "scale=0; ($ntraces+1)/2"| bc -l)
echo $middle
susum causal.su noncausal.su | \
suwind s=1 j=1 tmax=4 f1=0.0 | \
sushw key=f1,delrt,d2,fldr a=0.0,0.0,50,8000 | \
suwind key=tracl min=$middle max=$middle >> Trace.su
done
suwind tmin=1.5 tmax=3.0 f1=0.0 < Trace.su | \
sunormalize norm=max | \
supsgraph \
hbox=4 wbox=3 labelsize=10 x1beg=1.5 x1end=3.0 \
linecolor=red,green,blue,emerald,black f1=1.5 d1num=0.5 \
label1='time (s)' > shotTracesVt_T120_S_Dt120_F30.eps
suwind tmin=1.5 tmax=3.1 f1=0.0 < Trace.su | \
sunormalize norm=max | \
supswigp \
hbox=4 wbox=3 labelsize=10 x1beg=1.5 x1end=3.0 \
f1=1.5 d1num=0.5 linewidth=1 fill=0 f2=1 d2=1 \
label1='time (s)' label2='number of sources' > shotWigglesVt_T120_S_Dt120_F30.eps
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment