Skip to content
Snippets Groups Projects
Commit 44c3dd2d authored by Joeri Brackenhoff's avatar Joeri Brackenhoff
Browse files

Merge branch 'master' of https://github.com/JanThorbecke/OpenSource

parents ce16c33f f653696e
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 5573 additions and 113 deletions
.gitignore
+ 1
− 1
View file @ 44c3dd2d
......@@ -4,7 +4,7 @@ bin/*
*.su
*.bin
*/.DS*
*/*/.DS*
**/.DS*
.DS*
Make_include
fdelmodc_cuda.tgz
......
FFTlib/cc1fft.c
+ 16
− 9
View file @ 44c3dd2d
......@@ -61,10 +61,17 @@ void cc1fft(complex *data, int n, int sign)
REAL scl;
complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
#endif
int id;
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
pe = mp_my_threadnum();
......@@ -99,26 +106,26 @@ void cc1fft(complex *data, int n, int sign)
}
acmlcc1fft(sign, scl, inpl, n, data, 1, y, 1, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
if (n != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n);
Status = DftiCreateDescriptor(&handle[id], DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
nprev[id] = n;
}
if (sign < 0) {
Status = DftiComputeBackward(handle, data);
Status = DftiComputeBackward(handle[id], data);
}
else {
Status = DftiComputeForward(handle, data);
Status = DftiComputeForward(handle[id], data);
}
#else
cc1_fft(data, n, sign);
......
FFTlib/ccmfft.c
+ 16
− 9
View file @ 44c3dd2d
......@@ -64,11 +64,18 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
REAL scl;
complex *y;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
int j;
#endif
int id;
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
if (n1 != nprev) {
......@@ -99,29 +106,29 @@ void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
}
acmlccmfft(sign, scl, inpl, n2, n1, data, 1, ld1, y, 1, ld1, work, &isys);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
if (n1 != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n1);
Status = DftiCreateDescriptor(&handle[id], DFTI_SINGLE, DFTI_COMPLEX, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
nprev[id] = n1;
}
if (sign < 0) {
for (j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, &data[j*ld1]);
Status = DftiComputeBackward(handle[id], &data[j*ld1]);
}
}
else {
for (j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &data[j*ld1]);
Status = DftiComputeForward(handle[id], &data[j*ld1]);
}
}
#else
......
FFTlib/cr1fft.c
+ 18
− 10
View file @ 44c3dd2d
......@@ -57,12 +57,20 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
static REAL *work, *table, scale=1.0;
REAL scl;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
REAL *tmp;
MKL_LONG Status;
int i;
#endif
int id;
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
if (n != nprev) {
......@@ -98,34 +106,34 @@ void cr1fft(complex *cdata, REAL *rdata, int n, int sign)
}
acmlcrfft(one, n, rdata, work, &isys);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
if (n != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
Status = DftiCreateDescriptor(&handle[id], DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
Status = DftiSetValue(handle[id], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
nprev[id] = n;
}
tmp = (float *)malloc(n*sizeof(float));
Status = DftiComputeBackward(handle, (MKL_Complex8 *)cdata, tmp);
Status = DftiComputeBackward(handle[id], (MKL_Complex8 *)cdata, tmp);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeBackward FAIL\n");
......
FFTlib/crmfft.c
+ 19
− 11
View file @ 44c3dd2d
......@@ -71,12 +71,20 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
static REAL *work;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
REAL *tmp;
MKL_LONG Status;
int i, j;
#endif
int id;
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
nmp = mp_my_threadnum();
......@@ -138,37 +146,37 @@ void crmfft(complex *cdata, REAL *rdata, int n1, int n2, int ldc, int ldr, int s
memcpy(&rdata[j*ldr],&data[j*n1],n1*sizeof(REAL));
}
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
if (n1 != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
Status = DftiCreateDescriptor(&handle[id], DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
Status = DftiSetValue(handle[id], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
//This options is what we would like, but is depreciated in the future
//Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_REAL);
//Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_REAL);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
nprev[id] = n1;
}
tmp = (float *)malloc(n1*sizeof(float));
for (j=0; j<n2; j++) {
Status = DftiComputeBackward(handle, (MKL_Complex8 *)&cdata[j*ldc], tmp);
Status = DftiComputeBackward(handle[id], (MKL_Complex8 *)&cdata[j*ldc], tmp);
rdata[j*ldr] = tmp[0];
if (sign < 0) {
for (i=1; i<n1; i++) rdata[j*ldr+i] = -sign*tmp[n1-i];
......
FFTlib/rc1fft.c
+ 18
− 18
View file @ 44c3dd2d
......@@ -57,11 +57,18 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
static REAL *work, *table, scale=1.0;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
int i;
#endif
int id;
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
if (n != nprev) {
......@@ -102,40 +109,33 @@ void rc1fft(REAL *rdata, complex *cdata, int n, int sign)
cdata[n/2].i=0.0;
free(data);
#elif defined(MKL)
if (n != nprev) {
DftiFreeDescriptor(&handle);
if (n != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
Status = DftiCreateDescriptor(&handle[id], DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
Status = DftiSetValue(handle[id], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
/*
Status = DftiSetValue(handle, DFTI_FORWARD_DOMAIN, DFTI_REAL);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
*/
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n;
nprev[id] = n;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
Status = DftiComputeForward(handle[id], rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
......
FFTlib/rcmfft.c
+ 20
− 12
View file @ 44c3dd2d
......@@ -69,11 +69,19 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
static REAL *work;
REAL scl, *data;
#elif defined(MKL)
static DFTI_DESCRIPTOR_HANDLE handle=0;
static int nprev=0;
static DFTI_DESCRIPTOR_HANDLE handle[MAX_NUMTHREADS];
static int nprev[MAX_NUMTHREADS];
MKL_LONG Status;
int i,j;
int i, j;
#endif
int id;
#ifdef _OPENMP
id = omp_get_thread_num();
#else
id = 0;
#endif
#if defined(HAVE_LIBSCS)
nmp = mp_my_threadnum();
......@@ -130,39 +138,39 @@ void rcmfft(REAL *rdata, complex *cdata, int n1, int n2, int ldr, int ldc, int s
}
free(data);
#elif defined(MKL)
if (n1 != nprev) {
DftiFreeDescriptor(&handle);
if (n1 != nprev[id]) {
DftiFreeDescriptor(&handle[id]);
Status = DftiCreateDescriptor(&handle, DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
Status = DftiCreateDescriptor(&handle[id], DFTI_SINGLE, DFTI_REAL, 1, (MKL_LONG)n1);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCreateDescriptor FAIL\n");
}
Status = DftiSetValue(handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);
Status = DftiSetValue(handle[id], DFTI_PLACEMENT, DFTI_NOT_INPLACE);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiSetValue(handle, DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
Status = DftiSetValue(handle[id], DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);
if (! DftiErrorClass(Status, DFTI_NO_ERROR)) {
dfti_status_print(Status);
printf(" DftiSetValue FAIL\n");
}
Status = DftiCommitDescriptor(handle);
Status = DftiCommitDescriptor(handle[id]);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiCommitDescriptor FAIL\n");
}
nprev = n1;
nprev[id] = n1;
}
Status = DftiComputeForward(handle, rdata, (MKL_Complex8 *)cdata);
Status = DftiComputeForward(handle[id], rdata, (MKL_Complex8 *)cdata);
if(! DftiErrorClass(Status, DFTI_NO_ERROR)){
dfti_status_print(Status);
printf(" DftiComputeForward FAIL\n");
}
for (j=0; j<n2; j++) {
Status = DftiComputeForward(handle, &rdata[j*ldr], (MKL_Complex8 *)&cdata[j*ldc]);
Status = DftiComputeForward(handle[id], &rdata[j*ldr], (MKL_Complex8 *)&cdata[j*ldc]);
for (i=1; i<((n1-1)/2)+1; i++) {
cdata[j*ldc+i].i *= -sign;
}
......
fdelmodc3D/.vscode/settings.json 0 → 100644
+ 5
− 0
View file @ 44c3dd2d
{
"files.associations": {
"*.tcc": "c"
}
}
\ No newline at end of file
fdelmodc3D/acoustic4_3D.c 0 → 100644
+ 174
− 0
View file @ 44c3dd2d
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<assert.h>
#include"fdelmodc3D.h"
#define MIN(x,y) ((x) < (y) ? (x) : (y))
long applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, long verbose);
long storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, long verbose);
long reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, long ixsrc, long izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, long verbose);
long boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
long boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, long itime, long verbose);
long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float **src_nwav, float *vx, float *vy, float *vz, float *p, float *rox, float *roy, float *roz, float *l2m, long verbose)
{
/*********************************************************************
COMPUTATIONAL OVERVIEW OF THE 4th ORDER STAGGERED GRID:
The captial symbols T (=P) Txz,Vx,Vz represent the actual grid
The indices ix,iz are related to the T grid, so the capital T
symbols represent the actual modelled grid.
one cel (iz,ix)
|
V extra column of vx,txz
|
------- V
| txz vz| txz vz txz vz txz vz txz vz txz vz txz
| |
| vx t | vx t vx t vx t vx t vx t vx
-------
txz vz txz vz txz vz txz vz txz vz txz vz txz
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
| | | | | | |
txz vz txz Vz--Txz-Vz--Txz-Vz Txz-Vz txz vz txz
| | | | | | |
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
| | | | | | |
txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
| | | | | | |
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
| | | | | | |
txz vz txz Vz Txz-Vz Txz-Vz Txz-Vz txz vz txz
| | | | | | |
vx t vx T---Vx--T---Vx--T---Vx--T vx t vx
txz vz txz vz txz vz txz vz txz vz txz vz txz
vx t vx t vx t vx t vx t vx t vx
txz vz txz vz txz vz txz vz txz vz txz vz txz <--|
|
extra row of txz/vz |
AUTHOR:
Jan Thorbecke (janth@xs4all.nl)
The Netherlands
***********************************************************************/
float c1, c2;
long ix, iy, iz;
long n1, n2;
long ioXx, ioXy, ioXz, ioYx, ioYy, ioYz, ioZx, ioZy, ioZz, ioPx, ioPy, ioPz;
c1 = 9.0/8.0;
c2 = -1.0/24.0;
n1 = mod.naz;
n2 = mod.nay;
/* calculate vx for all grid points except on the virtual boundary*/
#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
#pragma ivdep
for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
for (iz=mod.ioXz; iz<mod.ieXz; iz++) {
vx[iy*n2*n1+ix*n1+iz] -= rox[iy*n2*n1+ix*n1+iz]*(
c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+(ix-1)*n1+iz]) +
c2*(p[iy*n2*n1+(ix+1)*n1+iz] - p[iy*n2*n1+(ix-2)*n1+iz]));
}
}
}
/* calculate vy for all grid points except on the virtual boundary*/
#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
for (iy=mod.ioYy; iy<mod.ieYy; iy++) {
#pragma ivdep
for (ix=mod.ioYx; ix<mod.ieYx; ix++) {
for (iz=mod.ioYz; iz<mod.ieYz; iz++) {
vy[iy*n2*n1+ix*n1+iz] -= roy[iy*n2*n1+ix*n1+iz]*(
c1*(p[iy*n2*n1+ix*n1+iz] - p[(iy-1)*n2*n1+ix*n1+iz]) +
c2*(p[(iy+1)*n2*n1+ix*n1+iz] - p[(iy-2)*n2*n1+ix*n1+iz]));
}
}
}
/* calculate vz for all grid points except on the virtual boundary */
#pragma omp for private (ix, iy, iz) schedule(guided,1)
for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
#pragma ivdep
for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
for (iz=mod.ioZz; iz<mod.ieZz; iz++) {
vz[iy*n2*n1+ix*n1+iz] -= roz[iy*n2*n1+ix*n1+iz]*(
c1*(p[iy*n2*n1+ix*n1+iz] - p[iy*n2*n1+ix*n1+iz-1]) +
c2*(p[iy*n2*n1+ix*n1+iz+1] - p[iy*n2*n1+ix*n1+iz-2]));
}
}
}
/* boundary condition clears velocities on boundaries */
boundariesP3D(mod, bnd, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, NULL, NULL, itime, verbose);
/* Add force source */
if (src.type > 5) {
applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, src_nwav, verbose);
}
/* this is needed because the P fields are not using tapered boundaries (bnd....=4) */
if (bnd.top==2) mod.ioPz += bnd.npml;
if (bnd.bot==2) mod.iePz -= bnd.npml;
if (bnd.fro==2) mod.ioPy += bnd.npml;
if (bnd.bac==2) mod.iePy -= bnd.npml;
if (bnd.lef==2) mod.ioPx += bnd.npml;
if (bnd.rig==2) mod.iePx -= bnd.npml;
/* calculate p/tzz for all grid points except on the virtual boundary */
#pragma omp for private (ix, iy, iz) schedule(guided,1)
//#pragma omp for private (ix, iz) schedule(dynamic)
#pragma ivdep
for (ix=mod.ioPx; ix<mod.iePx; ix++) {
#pragma ivdep
for (iy=mod.ioPy; iy<mod.iePy; iy++) {
for (iz=mod.ioPz; iz<mod.iePz; iz++) {
p[iy*n1*n2+ix*n1+iz] -= l2m[iy*n1*n2+ix*n1+iz]*(
c1*(vx[iy*n1*n2+(ix+1)*n1+iz] - vx[iy*n1*n2+ix*n1+iz]) +
c2*(vx[iy*n1*n2+(ix+2)*n1+iz] - vx[iy*n1*n2+(ix-1)*n1+iz]) +
c1*(vy[(iy+1)*n1*n2+ix*n1+iz] - vy[iy*n1*n2+ix*n1+iz]) +
c2*(vy[(iy+2)*n1*n2+ix*n1+iz] - vy[(iy-1)*n1*n2+ix*n1+iz]) +
c1*(vz[iy*n1*n2+ix*n1+iz+1] - vz[iy*n1*n2+ix*n1+iz]) +
c2*(vz[iy*n1*n2+ix*n1+iz+2] - vz[iy*n1*n2+ix*n1+iz-1]));
}
}
}
if (bnd.top==2) mod.ioPz -= bnd.npml;
if (bnd.bot==2) mod.iePz += bnd.npml;
if (bnd.fro==2) mod.ioPy -= bnd.npml;
if (bnd.bac==2) mod.iePy += bnd.npml;
if (bnd.lef==2) mod.ioPx -= bnd.npml;
if (bnd.rig==2) mod.iePx += bnd.npml;
/* Add stress source */
if (src.type < 6) {
applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, src_nwav, verbose);
}
/* Free surface: calculate free surface conditions for stresses */
/* check if there are sources placed on the free surface */
storeSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, verbose);
/* Free surface: calculate free surface conditions for stresses */
boundariesV3D(mod, bnd, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, NULL, NULL, itime, verbose);
/* restore source positions on the edge */
reStoreSourceOnSurface3D(mod, src, bnd, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, verbose);
return 0;
}
fdelmodc3D/applySource3D.c 0 → 100644
+ 365
− 0
View file @ 44c3dd2d
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<assert.h>
#include"fdelmodc3D.h"
void vmess(char *fmt, ...);
#define c1 (9.0/8.0)
#define c2 (-1.0/24.0)
/*********************************************************************
*
* Add's the source amplitude(s) to the grid.
*
* For the acoustic schemes, the source-type must not be txx tzz or txz.
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
*
**********************************************************************/
long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, long ixsrc, long iysrc, long izsrc, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *txy, float *tyz, float *rox, float *roy, float *roz, float *l2m, float **src_nwav, long verbose)
{
long is0, ibndz, ibndy, ibndx;
long isrc, ix, iy, iz, n1, n2;
long id1, id2, id3;
float src_ampl, time, scl, dt, sdx;
float Mxx, Myy, Mzz, Mxz, Myz, Mxy, rake;
static long first=1;
if (src.type==6) {
ibndz = mod.ioXz;
ibndy = mod.ioXy;
ibndx = mod.ioXx;
}
else if (src.type==7) {
ibndz = mod.ioZz;
ibndy = mod.ioZy;
ibndx = mod.ioZx;
}
else if (src.type==2) {
ibndz = mod.ioTz;
ibndy = mod.ioTy;
ibndx = mod.ioTx;
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
}
else {
ibndz = mod.ioPz;
ibndy = mod.ioPy;
ibndx = mod.ioPx;
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
if (bnd.bac==4 || bnd.fro==2) ibndy += bnd.ntap;
if (bnd.top==4 || bnd.top==2) ibndz += bnd.ntap;
}
n1 = mod.naz;
n2 = mod.nax;
dt = mod.dt;
sdx = 1.0/mod.dx;
/* special txz source activated? */
if ((bnd.top==1) && (src.type==2)) {
iz = izsrc + ibndz;
if (iz==ibndz) {
if (src.orient != 1) {
if (first) {
vmess("Only monopole Txz source allowed at surface. Reset to monopole");
first = 0;
}
src.orient=1;
}
}
}
/*
* for plane wave sources the sources are placed
* around the central shot position
* the first source position has an offset in x of is0
*
* itime = 0 corresponds with time=0
* itime = 1 corresponds with time=dt
* src[0] (the first sample) corresponds with time = 0
*/
is0 = -1*floor((src.n-1)/2);
#pragma omp for private (isrc, src_ampl, ix, iy, iz, time, id1, id2, id3, scl)
for (isrc=0; isrc<src.n; isrc++) {
src_ampl=0.0;
/* calculate the source position */
if (src.random || src.multiwav) {
ix = src.x[isrc] + ibndx;
iy = src.y[isrc] + ibndy;
iz = src.z[isrc] + ibndz;
}
else { /* plane wave and point sources */
ix = ixsrc + ibndx + is0 + isrc;
iy = iysrc + ibndy + is0 + isrc;
iz = izsrc + ibndz;
}
time = itime*dt - src.tbeg[isrc];
id1 = floor(time/dt);
id2 = id1+1;
/* delay not reached or no samples left in source wavelet? */
if ( (time < 0.0) || ( (itime*dt) >= src.tend[isrc]) ) continue;
// fprintf(stderr,"isrc=%li ix=%li iz=%li src.x=%li src.z=%li\n", isrc, ix, iz, src.x[isrc], src.z[isrc]);
if (!src.multiwav) { /* only one wavelet for all sources */
src_ampl = src_nwav[0][id1]*(id2-time/dt) + src_nwav[0][id2]*(time/dt-id1);
}
else { /* multi-wavelet sources */
src_ampl = src_nwav[isrc][id1]*(id2-time/dt) + src_nwav[isrc][id2]*(time/dt-id1);
}
if (src_ampl==0.0) continue;
if ( ((ix-ibndx)<0) || ((ix-ibndx)>mod.nx) ) continue; /* source outside grid */
if ( ((iy-ibndy)<0) || ((iy-ibndy)>mod.ny) ) continue; /* source outside grid */
if (verbose>=4 && itime==0) {
vmess("Source %li positioned at grid ix=%li iy=%li iz=%li",isrc, ix, iy, iz);
}
/* cosine squared windowing to reduce edge effects on shot arrays */
if ( (src.n>1) && src.window) {
scl = 1.0;
if (isrc < src.window) {
scl = cos(0.5*M_PI*(src.window - isrc)/src.window);
}
else if (isrc > src.n-src.window+1) {
scl = cos(0.5*M_PI*(src.window - (src.n-isrc+1))/src.window);
}
src_ampl *= scl*scl;
}
/* source scaling factor to compensate for discretisation */
/* old amplitude setting does not obey reciprocity */
// src_ampl *= rox[ix*n1+iz]*l2m[ix*n1+iz]/(dt);
/* in older version added factor 2.0 to be compliant with defined Green's functions in Marchenko algorithm */
/* this is now set to 1.0 */
src_ampl *= (1.0/mod.dx)*l2m[iy*n1*n2+ix*n1+iz];
if (verbose>5) {
vmess("Source %li at grid [ix=%li,iy=%li,iz=%li] at itime %li has value %e",isrc, ix, iy, iz, itime, src_ampl);
}
/* Force source */
if (src.type == 6) {
vx[iy*n1*n2+ix*n1+iz] += src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
//vx[ix*n1+iz] = 0.5*(vx[(ix+1)*n1+iz]+vx[(ix-1)*n1+iz])+src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
}
else if (src.type == 7) {
vz[iy*n1*n2+ix*n1+iz] += src_ampl*roz[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
/* stable implementation changes amplitude and more work is needed */
//vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
//vz[ix*n1+iz] = 0.25*(vz[ix*n1+iz-2]+vz[ix*n1+iz-1]+vz[ix*n1+iz]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
} /* src.type */
/* Stress source */
if (mod.ischeme <= 2) { /* Acoustic scheme */
/* Compressional source */
if (src.type == 1) {
if (src.orient != 1) src_ampl=src_ampl/mod.dx;
if (src.orient==1) { /* monopole */
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
}
else if (src.orient==2) { /* dipole +/- */
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
}
else if (src.orient==3) { /* dipole - + */
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
}
else if (src.orient==4) { /* dipole +/0/- */
if (iz > ibndz)
tzz[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
if (iz < mod.nz+ibndz-1)
tzz[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
}
else if (src.orient==5) { /* dipole + - */
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
}
}
}
else { /* Elastic scheme */
/* Compressional source */
if (src.type == 1) {
if (src.orient==1) { /* monopole */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
}
else if (src.orient==2) { /* dipole +/- */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
txx[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
tzz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
}
else if (src.orient==3) { /* dipole - + */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
txx[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
tzz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
}
else if (src.orient==4) { /* dipole +/0/- */
if (iz > ibndz) {
txx[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
tzz[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
}
if (iz < mod.nz+ibndz-1) {
txx[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
tzz[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
}
}
else if (src.orient==5) { /* dipole + - */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
txx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
tzz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
}
}
else if (src.type == 2) {
/* Txz source */
if ((iz == ibndz) && bnd.top==1) {
txz[iy*n1*n2+(ix-1)*n1+iz-1] += src_ampl;
txz[iy*n1*n2+ix*n1+iz-1] += src_ampl;
}
else {
txz[iy*n1*n2+ix*n1+iz] += src_ampl;
}
/* possible dipole orientations for a txz source */
if (src.orient == 2) { /* dipole +/- */
txz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
}
else if (src.orient == 3) { /* dipole - + */
txz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
}
else if (src.orient == 4) { /* dipole +/O/- */
/* correction: subtrace previous value to prevent z-1 values. */
txz[iy*n1*n2+ix*n1+iz] -= 2.0*src_ampl;
txz[iy*n1*n2+ix*n1+iz+1] += src_ampl;
}
else if (src.orient == 5) { /* dipole + - */
txz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
}
}
/* Tzz source */
else if(src.type == 3) {
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
}
/* Txx source */
else if(src.type == 4) {
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
}
/***********************************************************************
* pure potential shear S source (experimental)
* Curl S-pot = CURL(F) = dF_x/dz - dF_z/dx
***********************************************************************/
else if(src.type == 5) {
src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
if (src.orient == 3) src_ampl = -src_ampl;
/* first order derivatives */
vx[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vx[(iy-1)*n1*n2+ix*n1+iz-1] -= src_ampl*sdx;
vy[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vy[iy*n1*n2+(ix-1)*n1+iz-1] -= src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vz[(iy-1)*n1*n2+(ix-1)*n1+iz] += src_ampl*sdx;
/* second order derivatives */
/*
vx[ix*n1+iz] += c1*src_ampl*sdx;
vx[ix*n1+iz-1] -= c1*src_ampl*sdx;
vx[ix*n1+iz+1] += c2*src_ampl*sdx;
vx[ix*n1+iz-2] -= c2*src_ampl*sdx;
vz[ix*n1+iz] -= c1*src_ampl*sdx;
vz[(ix-1)*n1+iz] += c1*src_ampl*sdx;
vz[(ix+1)*n1+iz] -= c2*src_ampl*sdx;
vz[(ix-2)*n1+iz] += c2*src_ampl*sdx;
*/
/* determine second position of dipole */
if (src.orient == 2) { /* dipole +/- vertical */
iz += 1;
vx[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vx[(iy-1)*n1*n2+ix*n1+iz-1] += src_ampl*sdx;
vy[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vy[iy*n1*n2+(ix-1)*n1+iz-1] += src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vz[(iy-1)*n1*n2+(ix-1)*n1+iz] -= src_ampl*sdx;
}
else if (src.orient == 3) { /* dipole - + horizontal */
ix += 1;
vx[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vx[(iy-1)*n1*n2+ix*n1+iz-1] += src_ampl*sdx;
vy[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vy[iy*n1*n2+(ix-1)*n1+iz-1] += src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vz[(iy-1)*n1*n2+(ix-1)*n1+iz] -= src_ampl*sdx;
}
}
/***********************************************************************
* pure potential pressure P source (experimental)
* Divergence P-pot = DIV(F) = dF_x/dx + dF_z/dz
***********************************************************************/
else if(src.type == 8) {
src_ampl = src_ampl*rox[iy*n1*n2+ix*n1+iz]/(l2m[iy*n1*n2+ix*n1+iz]);
if (src.orient == 3) src_ampl = -src_ampl;
vx[iy*n1*n2+(ix+1)*n1+iz] += src_ampl*sdx;
vx[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vy[(iy+1)*n1*n2+ix*n1+iz] += src_ampl*sdx;
vy[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz+1] += src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz] -= src_ampl*sdx;
/* determine second position of dipole */
if (src.orient == 2) { /* dipole +/- */
iz += 1;
vx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl*sdx;
vx[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vy[(iy+1)*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vy[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz+1] -= src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
}
else if (src.orient == 3) { /* dipole - + */
ix += 1;
vx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl*sdx;
vx[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vy[(iy+1)*n1*n2+ix*n1+iz] -= src_ampl*sdx;
vy[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz+1] -= src_ampl*sdx;
vz[iy*n1*n2+ix*n1+iz] += src_ampl*sdx;
}
}
else if(src.type == 9) {
rake = 0.5*M_PI;
Mxx = -1.0*(sin(src.dip)*cos(rake)*sin(2.0*src.strike)+sin(src.dip*2.0)*sin(rake)*sin(src.strike)*sin(src.strike));
Mxz = -1.0*(cos(src.dip)*cos(rake)*cos(src.strike)+cos(src.dip*2.0)*sin(rake)*sin(src.strike));
Mzz = sin(src.dip*2.0)*sin(rake);
txx[iy*n1*n2+ix*n1+iz] -= Mxx*src_ampl;
tzz[iy*n1*n2+ix*n1+iz] -= Mzz*src_ampl;
txz[iy*n1*n2+ix*n1+iz] -= Mxz*src_ampl;
} /* src.type */
} /* ischeme */
} /* loop over isrc */
return 0;
}
fdelmodc3D/boundaries3D.c 0 → 100644
+ 2253
− 0
View file @ 44c3dd2d
This diff is collapsed.
fdelmodc3D/fdelmodc3D.c
+ 48
− 39
View file @ 44c3dd2d
......@@ -589,51 +589,41 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
// vx, vz, tzz, rox, roz, l2m, verbose);
// break;
case -1 : /* Acoustic dissipative media FD kernel */
acoustic4_qr(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, rox, roz, l2m, verbose);
vmess("Acoustic dissipative not yet available");
break;
case 1 : /* Acoustic FD kernel */
if (mod.iorder==2) {
acoustic2(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, rox, roz, l2m, verbose);
vmess("Acoustic order 2 not yet available");
}
else if (mod.iorder==4) {
if (mod.sh) {
acousticSH4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, rox, roz, l2m, verbose);
vmess("SH order 4 not yet available");
}
else {
acoustic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, rox, roz, l2m, verbose);
acoustic4_3D(mod, src, wav, bnd, it, ixsrc, iysrc, izsrc, src_nwav,
vx, vy, vz, tzz, rox, roy, roz, l2m, verbose);
}
}
else if (mod.iorder==6) {
acoustic6(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, rox, roz, l2m, verbose);
vmess("Acoustic order 6 not yet available");
}
break;
case 2 : /* Visco-Acoustic FD kernel */
viscoacoustic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, rox, roz, l2m, tss, tep, q, verbose);
vmess("Visco-Acoustic order 4 not yet available");
break;
case 3 : /* Elastic FD kernel */
if (mod.iorder==4) {
elastic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
vmess("Elastic order 4 not yet available");
}
else if (mod.iorder==6) {
elastic6(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
vmess("Elastic order 6 not yet available");
}
break;
case 4 : /* Visco-Elastic FD kernel */
viscoelastic4(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul,
tss, tep, tes, r, q, p, verbose);
vmess("Visco-Elastic order 4 not yet available");
break;
case 5 : /* Elastic FD kernel with S-velocity set to zero*/
elastic4dc(mod, src, wav, bnd, it, ixsrc, izsrc, src_nwav,
vx, vz, tzz, txx, txz, rox, roz, l2m, lam, mul, verbose);
vmess("DC-Elastic order 4 not yet available");
break;
}
......@@ -649,29 +639,33 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
/* Note that time step it=0 (t=0 for t**-fields t=-1/2 dt for v*-field) is not recorded */
isam = (it-rec.delay-itwritten)/rec.skipdt+1;
/* store time at receiver positions */
getRecTimes(mod, rec, bnd, it, isam, vx, vz, tzz, txx, txz,
l2m, rox, roz,
rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
getRecTimes3D(mod, rec, bnd, it, isam, vx, vy, vz,
tzz, tyy, txx, txz, txy, tyz,
l2m, rox, roy, roz,
rec_vx, rec_vy, rec_vz,
rec_txx, rec_tyy, rec_tzz, rec_txz, rec_txy, rec_tyz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
/* at the end of modeling a shot, write receiver array to output file(s) */
if (writeToFile && (it+rec.skipdt <= it1-1) ) {
fileno = ( ((it-rec.delay)/rec.skipdt)+1)/rec.nt;
writeRec(rec, mod, bnd, wav, ixsrc, izsrc, isam+1, ishot, fileno,
rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
writeRec3D(rec, mod, bnd, wav, ixsrc, iysrc, izsrc, isam+1, ishot, fileno,
rec_vx, rec_vy, rec_vz, rec_txx, rec_tyy, rec_tzz, rec_txz, rec_tyz, rec_txy,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
}
}
/* write snapshots to output file(s) */
if (sna.nsnap) {
writeSnapTimes(mod, sna, bnd, wav, ixsrc, izsrc, it, vx, vz, tzz, txx, txz, verbose);
writeSnapTimes3D(mod, sna, bnd, wav, ixsrc, iysrc, izsrc, it,
vx, vy, vz, tzz, tyy, txx, txz, tyz, txy, verbose);
}
/* calculate beams */
if(sna.beam) {
getBeamTimes(mod, sna, vx, vz, tzz, txx, txz,
beam_vx, beam_vz, beam_txx, beam_tzz, beam_txz,
getBeamTimes3D(mod, sna, vx, vy, vz, tzz, tyy, txx, txz, tyz, txy,
beam_vx, beam_vy, beam_vz, beam_txx, beam_tyy, beam_tzz, beam_txz, beam_tyz, beam_txy,
beam_p, beam_pp, beam_ss, verbose);
}
}
......@@ -696,27 +690,29 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
if (fileno) fileno++;
if (rec.scale==1) { /* scale receiver with distance src-rcv */
float xsrc, zsrc, Rrec, rdx, rdz;
float xsrc, ysrc, zsrc, Rrec, rdx, rdy, rdz;
long irec;
xsrc=mod.x0+mod.dx*ixsrc;
ysrc=mod.y0+mod.dy*iysrc;
zsrc=mod.z0+mod.dz*izsrc;
for (irec=0; irec<rec.n; irec++) {
rdx=mod.x0+rec.xr[irec]-xsrc;
rdy=mod.y0+rec.yr[irec]-ysrc;
rdz=mod.z0+rec.zr[irec]-zsrc;
Rrec = sqrt(rdx*rdx+rdz*rdz);
fprintf(stderr,"Rec %li is scaled with distance %f R=%.2f,%.2f S=%.2f,%.2f\n", irec, Rrec,rdx,rdz,xsrc,zsrc);
Rrec = sqrt(rdx*rdx+rdy*rdy+rdz*rdz);
fprintf(stderr,"Rec %li is scaled with distance %f R=%.2f,%.2f,%.2f S=%.2f,%.2f,%.2f\n", irec, Rrec,rdx,rdy,rdz,xsrc,ysrc,zsrc);
for (it=0; it<rec.nt; it++) {
rec_p[irec*rec.nt+it] *= sqrt(Rrec);
}
}
}
writeRec(rec, mod, bnd, wav, ixsrc, izsrc, isam+1, ishot, fileno,
rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
writeRec3D(rec, mod, bnd, wav, ixsrc, iysrc, izsrc, isam+1, ishot, fileno,
rec_vx, rec_vy, rec_vz, rec_txx, rec_tyy, rec_tzz, rec_txz, rec_tyz, rec_txy,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
writeBeams(mod, sna, ixsrc, izsrc, ishot, fileno,
beam_vx, beam_vz, beam_txx, beam_tzz, beam_txz,
beam_p, beam_pp, beam_ss, verbose);
writeBeams3D(mod, sna, ixsrc, iysrc, izsrc, ishot, fileno,
beam_vx, beam_vy, beam_vz, beam_txx, beam_tyy, beam_tzz, beam_txz, beam_tyz, beam_txy,
beam_p, beam_pp, beam_ss, verbose);
} /* end of loop over number of shots */
......@@ -731,20 +727,26 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
initargs(argc,argv); /* this will free the arg arrays declared */
free(rox);
free(roy);
free(roz);
free(l2m);
free(src_nwav[0]);
free(src_nwav);
free(vx);
free(vy);
free(vz);
free(tzz);
freeStoreSourceOnSurface();
if (rec.type.vz) free(rec_vz);
if (rec.type.vy) free(rec_vy);
if (rec.type.vx) free(rec_vx);
if (rec.type.p) free(rec_p);
if (rec.type.txx) free(rec_txx);
if (rec.type.tyy) free(rec_tyy);
if (rec.type.tzz) free(rec_tzz);
if (rec.type.txz) free(rec_txz);
if (rec.type.tyz) free(rec_tyz);
if (rec.type.txy) free(rec_txy);
if (rec.type.pp) free(rec_pp);
if (rec.type.ss) free(rec_ss);
if (rec.type.ud) {
......@@ -753,11 +755,15 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
}
if(sna.beam) {
if (sna.type.vz) free(beam_vz);
if (sna.type.vy) free(beam_vy);
if (sna.type.vx) free(beam_vx);
if (sna.type.p) free(beam_p);
if (sna.type.txx) free(beam_txx);
if (sna.type.tyy) free(beam_tyy);
if (sna.type.tzz) free(beam_tzz);
if (sna.type.txz) free(beam_txz);
if (sna.type.tyz) free(beam_tyz);
if (sna.type.txy) free(beam_txy);
if (sna.type.pp) free(beam_pp);
if (sna.type.ss) free(beam_ss);
}
......@@ -771,7 +777,10 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, iysrc, izsrc, it, src_nwav, verbos
free(lam);
free(mul);
free(txz);
free(tyz);
free(txy);
free(txx);
free(tyy);
}
if (mod.ischeme==4) {
free(tss);
......
fdelmodc3D/getBeamTimes3D.c 0 → 100644
+ 246
− 0
View file @ 44c3dd2d
#define _FILE_OFFSET_BITS 64
#define _LARGEFILE_SOURCE
#define _LARGEFILE64_SOURCE
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include <string.h>
#include "segy.h"
#include "fdelmodc3D.h"
/**
* getBeamTimes: stores energy fields (beams) in arrays at certain time steps
* writeBeams: writes the stored fields to output file(s)
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
FILE *fileOpen(char *file, char *ext, int append);
int traceWrite(segy *hdr, float *data, int n, FILE *fp);
void name_ext(char *filename, char *extension);
void vmess(char *fmt, ...);
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
long getBeamTimes3D(modPar mod, snaPar sna, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *tyz, float *txy,
float *beam_vx, float *beam_vy, float *beam_vz, float *beam_txx, float *beam_tyy, float *beam_tzz, float *beam_txz, float *beam_tyz, float *beam_txy,
float *beam_p, float *beam_pp, float *beam_ss, long verbose)
{
long n1, n2, ibndx, ibndy, ibndz, ixs, iys, izs, ize, i, j, l;
long ix, iy, iz, ix2, iy2, iz2;
float sdx, s, p;
ibndx = mod.ioPx;
ibndy = mod.ioPy;
ibndz = mod.ioPz;
n1 = mod.naz;
n2 = mod.nax;
sdx = 1.0/mod.dx;
izs = sna.z1+ibndx;
ize = sna.z2+ibndz;
for (iys=sna.y1, l=0; iys<=sna.y2; iys+=sna.skipdy, l++) {
iy = iys+ibndy;
iy2 = iy+1;
for (ixs=sna.x1, i=0; ixs<=sna.x2; ixs+=sna.skipdx, i++) {
ix = ixs+ibndx;
ix2 = ix+1;
if (sna.type.vx) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_vx[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(vx[iy*n1*n2+ix2*n1+iz]*vx[iy*n1*n2+ix2*n1+iz]);
}
}
if (sna.type.vy) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_vy[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(vy[iy2*n1*n2+ix*n1+iz]*vx[iy2*n1*n2+ix*n1+iz]);
}
}
if (sna.type.vz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_vz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(vz[iy*n1*n2+ix*n1+iz+1]*vz[iy*n1*n2+ix*n1+iz+1]);
}
}
if (sna.type.p) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_p[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tzz[iy*n1*n2+ix*n1+iz]*tzz[iy*n1*n2+ix*n1+iz]);
}
}
if (sna.type.tzz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_tzz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tzz[iy*n1*n2+ix*n1+iz]*tzz[iy*n1*n2+ix*n1+iz]);
}
}
if (sna.type.tyy) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_tyy[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tyy[iy*n1*n2+ix*n1+iz]*tyy[iy*n1*n2+ix*n1+iz]);
}
}
if (sna.type.txx) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_txx[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(txx[iy*n1*n2+ix*n1+iz]*txx[iy*n1*n2+ix*n1+iz]);
}
}
if (sna.type.txz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_txz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(txz[iy*n1*n2+ix2*n1+iz+1]*txz[iy*n1*n2+ix2*n1+iz+1]);
}
}
if (sna.type.tyz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_tyz[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(tyz[iy2*n1*n2+ix*n1+iz+1]*tyz[iy2*n1*n2+ix*n1+iz+1]);
}
}
if (sna.type.txz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
beam_txy[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(txy[iy2*n1*n2+ix2*n1+iz]*txy[iy2*n1*n2+ix2*n1+iz]);
}
}
/* calculate divergence of velocity field */
if (sna.type.pp) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
iz2 = iz+1;
p = sdx*((vx[iy*n1*n2+ix2*n1+iz]-vx[iy*n1*n2+ix*n1+iz])+
(vy[iy2*n1*n2+ix*n1+iz]-vy[iy*n1*n2+ix*n1+iz])+
(vz[iy*n1*n2+ix*n1+iz2]-vz[iy*n1*n2+ix*n1+iz]));
beam_pp[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(p*p);
}
}
/* calculate rotation of velocity field */
if (sna.type.ss) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
iz2 = iz+1;
s = sdx*((vx[iy2*n1*n2+ix2*n1+iz2]-vx[iy*n1*n2+ix2*n1+iz])-
(vy[iy2*n1*n2+ix2*n1+iz2]-vy[iy2*n1*n2+ix*n1+iz])-
(vz[iy2*n1*n2+ix2*n1+iz2]-vz[iy*n1*n2+ix*n1+iz2]));
beam_ss[l*sna.nz*sna.nz+i*sna.nz+j] += sqrt(s*s);
}
}
}
}
return 0;
}
long writeBeams3D(modPar mod, snaPar sna, long ixsrc, long iysrc, long izsrc, long ishot, long fileno,
float *beam_vx, float *beam_vy, float *beam_vz, float *beam_txx, float *beam_tyy, float *beam_tzz, float *beam_txz, float *beam_tyz, float *beam_txy,
float *beam_p, float *beam_pp, float *beam_ss, long verbose)
{
FILE *fpvx, *fpvy, *fpvz, *fptxx, *fptyy, *fptzz, *fptxz, *fptyz, *fptxy, *fpp, *fppp, *fpss;
long append;
long ix, iy;
char number[16], filename[1024];
segy hdr;
if (sna.beam==0) return 0;
/* all beam snapshots are written to the same output file(s) */
if (ishot) append=1;
else append=0;
strcpy(filename, sna.file_beam);
if (fileno) {
sprintf(number,"_%03d",fileno);
name_ext(filename, number);
}
if (verbose>2) vmess("Writing beam data to file %s", filename);
if (sna.type.vx) fpvx = fileOpen(filename, "_bvx", (int)append);
if (sna.type.vy) fpvy = fileOpen(filename, "_bvy", (int)append);
if (sna.type.vz) fpvz = fileOpen(filename, "_bvz", (int)append);
if (sna.type.p) fpp = fileOpen(filename, "_bp", (int)append);
if (sna.type.txx) fptxx = fileOpen(filename, "_btxx", (int)append);
if (sna.type.tyy) fptyy = fileOpen(filename, "_btyy", (int)append);
if (sna.type.tzz) fptzz = fileOpen(filename, "_btzz", (int)append);
if (sna.type.txz) fptxz = fileOpen(filename, "_btxz", (int)append);
if (sna.type.tyz) fptyz = fileOpen(filename, "_btyz", (int)append);
if (sna.type.txy) fptxy = fileOpen(filename, "_btxy", (int)append);
if (sna.type.pp) fppp = fileOpen(filename, "_bpp", (int)append);
if (sna.type.ss) fpss = fileOpen(filename, "_bss", (int)append);
memset(&hdr,0,TRCBYTES);
hdr.dt = 1000000*(mod.dt);
hdr.scalco = -1000;
hdr.scalel = -1000;
hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
hdr.sy = 1000*(mod.y0+iysrc*mod.dy);
hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
hdr.fldr = ishot+1;
hdr.trid = 1;
hdr.ns = sna.nz;
hdr.trwf = sna.nx*sna.ny;
hdr.ntr = sna.nx*sna.ny;
hdr.f1 = sna.z1*mod.dz+mod.z0;
hdr.f2 = sna.x1*mod.dx+mod.x0;
hdr.d1 = mod.dz*sna.skipdz;
hdr.d2 = mod.dx*sna.skipdx;
for (iy=0; iy<sna.ny; iy++) {
for (ix=0; ix<sna.nx; ix++) {
hdr.tracf = iy*sna.nx+ix+1;
hdr.tracl = iy*sna.nx+ix+1;
hdr.gx = 1000*(mod.x0+(sna.x1+ix)*mod.dx);
hdr.gy = 1000*(mod.y0+(sna.y1+ix)*mod.dy);
if (sna.type.vx) {
traceWrite( &hdr, &beam_vx[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpvx) ;
}
if (sna.type.vy) {
traceWrite( &hdr, &beam_vy[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpvy) ;
}
if (sna.type.vz) {
traceWrite( &hdr, &beam_vz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpvz) ;
}
if (sna.type.p) {
traceWrite( &hdr, &beam_p[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpp) ;
}
if (sna.type.tzz) {
traceWrite( &hdr, &beam_tzz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptzz) ;
}
if (sna.type.tyy) {
traceWrite( &hdr, &beam_tyy[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptyy) ;
}
if (sna.type.txx) {
traceWrite( &hdr, &beam_txx[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptxx) ;
}
if (sna.type.txz) {
traceWrite( &hdr, &beam_txz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptxz) ;
}
if (sna.type.txy) {
traceWrite( &hdr, &beam_txy[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptxy) ;
}
if (sna.type.tyz) {
traceWrite( &hdr, &beam_tyz[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fptyz) ;
}
if (sna.type.pp) {
traceWrite( &hdr, &beam_pp[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fppp) ;
}
if (sna.type.ss) {
traceWrite( &hdr, &beam_ss[iy*sna.nx*sna.nz+ix*sna.nz], (int)sna.nz, fpss) ;
}
}
}
if (sna.type.vx) fclose(fpvx);
if (sna.type.vy) fclose(fpvy);
if (sna.type.vz) fclose(fpvz);
if (sna.type.p) fclose(fpp);
if (sna.type.txx) fclose(fptxx);
if (sna.type.tyy) fclose(fptyy);
if (sna.type.tzz) fclose(fptzz);
if (sna.type.txz) fclose(fptxz);
if (sna.type.tyz) fclose(fptyz);
if (sna.type.txy) fclose(fptxy);
if (sna.type.pp) fclose(fppp);
if (sna.type.ss) fclose(fpss);
return 0;
}
\ No newline at end of file
fdelmodc3D/getRecTimes3D.c 0 → 100644
+ 514
− 0
View file @ 44c3dd2d
This diff is collapsed.
fdelmodc3D/writeRec3D.c 0 → 100644
+ 255
− 0
View file @ 44c3dd2d
#define _FILE_OFFSET_BITS 64
#define _LARGEFILE_SOURCE
#define _LARGEFILE64_SOURCE
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include <string.h>
#include "segy.h"
#include "fdelmodc3D.h"
#include <genfft.h>
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
float r,i;
} complex;
typedef struct _dcomplexStruct { /* complex number */
double r,i;
} dcomplex;
#endif/* complex */
/**
* Writes the receiver array(s) to output file(s)
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
FILE *fileOpen(char *file, char *ext, int append);
int traceWrite(segy *hdr, float *data, int n, FILE *fp) ;
void name_ext(char *filename, char *extension);
void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
int nkx, float dx, int nt, float dt, float fmin, float fmax,
float cp, float rho, int vznorm, int verbose);
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
long writeRec3D(recPar rec, modPar mod, bndPar bnd, wavPar wav, long ixsrc, long iysrc, long izsrc, long nsam, long ishot, long fileno,
float *rec_vx, float *rec_vy, float *rec_vz, float *rec_txx, float *rec_tyy, float *rec_tzz, float *rec_txz, float *rec_tyz, float *rec_txy,
float *rec_p, float *rec_pp, float *rec_ss, float *rec_udp, float *rec_udvz, long verbose)
{
FILE *fpvx, *fpvy, *fpvz, *fptxx, *fptyy, *fptzz, *fptxz, *fptyz, *fptxy, *fpp, *fppp, *fpss, *fpup, *fpdown;
float *rec_up, *rec_down, *trace, *rec_vze, *rec_pe;
float dx, dy, dt, cp, rho, fmin, fmax;
complex *crec_vz, *crec_p, *crec_up, *crec_dw;
long irec, ntfft, nfreq, nkx, nky, xorig, yorig, ix, iy, iz, it, ibndx, ibndy;
long append, vznorm, sx, sy;
double ddt;
char number[16], filename[1024];
segy hdr;
if (!rec.n) return 0;
if (ishot) append=1;
else append=0;
/* if the total number of samples exceeds rec_ntsam then a new (numbered) file is opened */
/* fileno has a non-zero value (from fdelmodc.c) if the number of samples exceeds rec_ntsam. */
strcpy(filename, rec.file_rcv);
if (fileno) {
sprintf(number,"_%03d",fileno);
name_ext(filename, number);
}
#ifdef MPI
sx = (long)mod.x0+ixsrc*mod.dx;
sprintf(number,"_%06d",sx);
name_ext(filename, number);
#endif
if (verbose>2) vmess("Writing receiver data to file %s", filename);
if (nsam != rec.nt && verbose) vmess("Number of samples written to last file = %li",nsam);
memset(&hdr,0,TRCBYTES);
ddt = (double)mod.dt;/* to avoid rounding in 32 bit precision */
dt = (float)ddt*rec.skipdt;
dx = (rec.x[1]-rec.x[0])*mod.dx;
dy = (rec.y[1]-rec.y[0])*mod.dy;
hdr.dt = (unsigned short)lround((((double)1.0e6*ddt*rec.skipdt)));
hdr.scalco = -1000;
hdr.scalel = -1000;
hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
hdr.sy = 1000*(mod.y0+ixsrc*mod.dy);
hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
hdr.selev = (int)(-1000.0*(mod.z0+izsrc*mod.dz));
hdr.fldr = ishot+1;
hdr.trid = 1;
hdr.ns = nsam;
hdr.trwf = rec.n;
hdr.ntr = rec.n;
if (mod.grid_dir) { /* reverse time modeling */
hdr.f1 = (-mod.nt+1)*mod.dt;
}
else {
hdr.f1 = 0.0;
}
hdr.d1 = mod.dt*rec.skipdt;
hdr.d2 = (rec.x[1]-rec.x[0])*mod.dx;
hdr.f2 = mod.x0+rec.x[0]*mod.dx;
if (rec.type.vx) fpvx = fileOpen(filename, "_rvx", (int)append);
if (rec.type.vy) fpvy = fileOpen(filename, "_rvy", (int)append);
if (rec.type.vz) fpvz = fileOpen(filename, "_rvz", (int)append);
if (rec.type.p) fpp = fileOpen(filename, "_rp", (int)append);
if (rec.type.txx) fptxx = fileOpen(filename, "_rtxx", (int)append);
if (rec.type.tyy) fptyy = fileOpen(filename, "_rtyy", (int)append);
if (rec.type.tzz) fptzz = fileOpen(filename, "_rtzz", (int)append);
if (rec.type.txz) fptxz = fileOpen(filename, "_rtxz", (int)append);
if (rec.type.tyz) fptyz = fileOpen(filename, "_rtyz", (int)append);
if (rec.type.txy) fptxy = fileOpen(filename, "_rtxy", (int)append);
if (rec.type.pp) fppp = fileOpen(filename, "_rpp", (int)append);
if (rec.type.ss) fpss = fileOpen(filename, "_rss", (int)append);
/* decomposed wavefield */
if (rec.type.ud && (mod.ischeme==1 || mod.ischeme==2) ) {
fpup = fileOpen(filename, "_ru", (int)append);
fpdown = fileOpen(filename, "_rd", (int)append);
ntfft = optncr(nsam);
nfreq = ntfft/2+1;
fmin = 0.0;
fmax = wav.fmax;
nkx = optncc(2*mod.nax);
nky = optncc(2*mod.nay);
ibndx = mod.ioPx;
ibndy = mod.ioPy;
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
if (bnd.fro==4 || bnd.fro==2) ibndy += bnd.ntap;
cp = rec.cp;
rho = rec.rho;
if (rec.type.ud==2) vznorm=1;
else vznorm=0;
if (verbose) vmess("Decomposition array at z=%.2f with cp=%.2f rho=%.2f", rec.zr[0]+mod.z0, cp, rho);
rec_up = (float *)calloc(ntfft*nkx*nky,sizeof(float));
rec_down= (float *)calloc(ntfft*nkx*nky,sizeof(float));
crec_vz = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
crec_p = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
crec_up = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
crec_dw = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
rec_vze = rec_up;
rec_pe = rec_down;
/* copy input data into extended arrays with padded zeroes */
for (iy=0; iy<mod.nay; iy++) {
for (ix=0; ix<mod.nax; ix++) {
memcpy(&rec_vze[iy*mod.nax*ntfft+ix*ntfft],&rec_udvz[iy*mod.nax*rec.nt+ix*rec.nt],nsam*sizeof(float));
memcpy(&rec_pe[iy*mod.nax*ntfft+ix*ntfft], &rec_udp[iy*mod.nax*rec.nt+ix*rec.nt], nsam*sizeof(float));
}
}
/* transform from t-x to kx-w */
xorig = ixsrc+ibndx;
xt2wkx(rec_vze, crec_vz, ntfft, nkx, ntfft, nkx, xorig);
xt2wkx(rec_pe, crec_p, ntfft, nkx, ntfft, nkx, xorig);
/* apply decomposition operators */
kxwdecomp(crec_p, crec_vz, crec_up, crec_dw,
nkx, mod.dx, nsam, dt, fmin, fmax, cp, rho, vznorm, verbose);
/* transform back to t-x */
wkx2xt(crec_up, rec_up, ntfft, nkx, nkx, ntfft, xorig);
wkx2xt(crec_dw, rec_down, ntfft, nkx, nkx, ntfft, xorig);
/* reduce array to rec.nt samples rec.n traces */
for (irec=0; irec<rec.n; irec++) {
ix = rec.x[irec]+ibndx;
iy = rec.y[irec]+ibndy;
for (it=0; it<rec.nt; it++) {
rec_up[irec*rec.nt+it] = rec_up[iy*nkx*ntfft+ix*ntfft+it];
rec_down[irec*rec.nt+it] = rec_down[iy*nkx*ntfft+ix*ntfft+it];
}
}
free(crec_vz);
free(crec_p);
free(crec_up);
free(crec_dw);
}
if (rec.type.ud && (mod.ischeme==3 || mod.ischeme==4) ) {
}
for (irec=0; irec<rec.n; irec++) {
hdr.tracf = irec+1;
hdr.tracl = ishot*rec.n+irec+1;
hdr.gx = 1000*(mod.x0+rec.x[irec]*mod.dx);
hdr.gy = 1000*(mod.y0+rec.y[irec]*mod.dy);
hdr.offset = (rec.x[irec]-ixsrc)*mod.dx;
hdr.gelev = (int)(-1000*(mod.z0+rec.z[irec]*mod.dz));
if (rec.type.vx) {
traceWrite( &hdr, &rec_vx[irec*rec.nt], (int)nsam, fpvx) ;
}
if (rec.type.vy) {
traceWrite( &hdr, &rec_vy[irec*rec.nt], (int)nsam, fpvy) ;
}
if (rec.type.vz) {
traceWrite( &hdr, &rec_vz[irec*rec.nt], (int)nsam, fpvz) ;
}
if (rec.type.p) {
traceWrite( &hdr, &rec_p[irec*rec.nt], (int)nsam, fpp) ;
}
if (rec.type.txx) {
traceWrite( &hdr, &rec_txx[irec*rec.nt], (int)nsam, fptxx) ;
}
if (rec.type.tyy) {
traceWrite( &hdr, &rec_tyy[irec*rec.nt], (int)nsam, fptyy) ;
}
if (rec.type.tzz) {
traceWrite( &hdr, &rec_tzz[irec*rec.nt], (int)nsam, fptzz) ;
}
if (rec.type.txz) {
traceWrite( &hdr, &rec_txz[irec*rec.nt], (int)nsam, fptxz) ;
}
if (rec.type.txy) {
traceWrite( &hdr, &rec_txy[irec*rec.nt], (int)nsam, fptxy) ;
}
if (rec.type.tyz) {
traceWrite( &hdr, &rec_tyz[irec*rec.nt], (int)nsam, fptyz) ;
}
if (rec.type.pp) {
traceWrite( &hdr, &rec_pp[irec*rec.nt], (int)nsam, fppp) ;
}
if (rec.type.ss) {
traceWrite( &hdr, &rec_ss[irec*rec.nt], (int)nsam, fpss) ;
}
if (rec.type.ud && mod.ischeme==1) {
traceWrite( &hdr, &rec_up[irec*rec.nt], (int)nsam, fpup) ;
traceWrite( &hdr, &rec_down[irec*rec.nt], (int)nsam, fpdown) ;
}
}
if (rec.type.vx) fclose(fpvx);
if (rec.type.vy) fclose(fpvy);
if (rec.type.vz) fclose(fpvz);
if (rec.type.p) fclose(fpp);
if (rec.type.txx) fclose(fptxx);
if (rec.type.tyy) fclose(fptyy);
if (rec.type.tzz) fclose(fptzz);
if (rec.type.txz) fclose(fptxz);
if (rec.type.txy) fclose(fptxy);
if (rec.type.tyz) fclose(fptyz);
if (rec.type.pp) fclose(fppp);
if (rec.type.ss) fclose(fpss);
if (rec.type.ud) {
fclose(fpup);
fclose(fpdown);
free(rec_up);
free(rec_down);
}
return 0;
}
fdelmodc3D/writeSnapTimes3D.c 0 → 100644
+ 261
− 0
View file @ 44c3dd2d
#define _FILE_OFFSET_BITS 64
#define _LARGEFILE_SOURCE
#define _LARGEFILE64_SOURCE
#define ISODD(n) ((n) & 01)
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <math.h>
#include <string.h>
#include "par.h"
#include "segy.h"
#include "fdelmodc3D.h"
/**
* Writes gridded wavefield(s) at a desired time to output file(s)
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
FILE *fileOpen(char *file, char *ext, int append);
int traceWrite(segy *hdr, float *data, int n, FILE *fp);
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
long writeSnapTimes3D(modPar mod, snaPar sna, bndPar bnd, wavPar wav, long ixsrc, long iysrc, long izsrc, long itime, float *vx, float *vy, float *vz, float *tzz, float *tyy, float *txx, float *txz, float *tyz, float *txy, long verbose)
{
FILE *fpvx, *fpvy, *fpvz, *fptxx, *fptyy, *fptzz, *fptxz, *fptyz, *fptxy, *fpp, *fppp, *fpss;
long append, isnap;
static long first=1;
long n1, n2, ibndx, ibndy, ibndz, ixs, iys, izs, ize, i, j, l;
long ix, iy, iz, ix2, iy2;
float *snap, sdx, stime;
segy hdr;
if (sna.nsnap==0) return 0;
ibndx = mod.ioXx;
ibndy = mod.ioXy;
ibndz = mod.ioXz;
n1 = mod.naz;
n2 = mod.nax;
sdx = 1.0/mod.dx;
if (sna.withbnd) {
sna.nz=mod.naz;
sna.z1=0;
sna.z2=mod.naz-1;
sna.skipdz=1;
sna.ny=mod.nax;
sna.y1=0;
sna.y2=mod.nay-1;
sna.skipdy=1;
sna.nx=mod.nax;
sna.x1=0;
sna.x2=mod.nax-1;
sna.skipdx=1;
}
/* check if this itime is a desired snapshot time */
if ( (((itime-sna.delay) % sna.skipdt)==0) &&
(itime >= sna.delay) &&
(itime <= sna.delay+(sna.nsnap-1)*sna.skipdt) ) {
isnap = NINT((itime-sna.delay)/sna.skipdt);
if (mod.grid_dir) stime = (-wav.nt+1+itime+1)*mod.dt; /* reverse time modeling */
else stime = itime*mod.dt;
if (verbose) vmess("Writing snapshot(%li) at time=%.4f", isnap+1, stime);
if (first) {
append=0;
first=0;
}
else {
append=1;
}
if (sna.type.vx) fpvx = fileOpen(sna.file_snap, "_svx", (int)append);
if (sna.type.vy) fpvy = fileOpen(sna.file_snap, "_svy", (int)append);
if (sna.type.vz) fpvz = fileOpen(sna.file_snap, "_svz", (int)append);
if (sna.type.p) fpp = fileOpen(sna.file_snap, "_sp", (int)append);
if (sna.type.txx) fptxx = fileOpen(sna.file_snap, "_stxx", (int)append);
if (sna.type.tyy) fptyy = fileOpen(sna.file_snap, "_styy", (int)append);
if (sna.type.tzz) fptzz = fileOpen(sna.file_snap, "_stzz", (int)append);
if (sna.type.txz) fptxz = fileOpen(sna.file_snap, "_stxz", (int)append);
if (sna.type.tyz) fptyz = fileOpen(sna.file_snap, "_styz", (int)append);
if (sna.type.txy) fptxy = fileOpen(sna.file_snap, "_stxy", (int)append);
if (sna.type.pp) fppp = fileOpen(sna.file_snap, "_spp", (int)append);
if (sna.type.ss) fpss = fileOpen(sna.file_snap, "_sss", (int)append);
memset(&hdr,0,TRCBYTES);
hdr.dt = 1000000*(sna.skipdt*mod.dt);
hdr.ungpow = (sna.delay*mod.dt);
hdr.scalco = -1000;
hdr.scalel = -1000;
hdr.sx = 1000*(mod.x0+ixsrc*mod.dx);
hdr.sy = 1000*(mod.y0+iysrc*mod.dy);
hdr.sdepth = 1000*(mod.z0+izsrc*mod.dz);
hdr.fldr = isnap+1;
hdr.trid = 1;
hdr.ns = sna.nz;
hdr.trwf = sna.nx*sna.nx;
hdr.ntr = (isnap+1)*sna.nx;
hdr.f1 = sna.z1*mod.dz+mod.z0;
hdr.f2 = sna.x1*mod.dx+mod.x0;
hdr.d1 = mod.dz*sna.skipdz;
hdr.d2 = mod.dx*sna.skipdx;
if (sna.withbnd) {
if ( !ISODD(bnd.top)) hdr.f1 = mod.z0 - bnd.ntap*mod.dz;
if ( !ISODD(bnd.lef)) hdr.f2 = mod.x0 - bnd.ntap*mod.dx;
//if ( !ISODD(bnd.rig)) ;
//if ( !ISODD(bnd.bot)) store=1;
}
/***********************************************************************
* vx velocities have one sample less in x-direction
* vz velocities have one sample less in z-direction
* txz stresses have one sample less in z-direction and x-direction
***********************************************************************/
snap = (float *)malloc(sna.nz*sizeof(float));
/* Decimate, with skipdx and skipdz, the number of gridpoints written to file
and write to file. */
for (iys=sna.y1, l=0; iys<=sna.x2; iys+=sna.skipdx, l++) {
for (ixs=sna.x1, i=0; ixs<=sna.x2; ixs+=sna.skipdx, i++) {
hdr.tracf = l*sna.nx+i+1;
hdr.tracl = isnap*sna.nx*sna.ny+l*sna.nx+i+1;
hdr.gx = 1000*(mod.x0+ixs*mod.dx);
hdr.gy = 1000*(mod.y0+ixs*mod.dy);
ix = ixs+ibndx;
ix2 = ix+1;
iy = iys+ibndy;
iy2 = iy+1;
izs = sna.z1+ibndz;
ize = sna.z2+ibndz;
if (sna.withbnd) {
izs = 0;
ize = sna.z2;
ix = ixs;
ix2 = ix;
iy = iys;
iy2 = iy;
if (sna.type.vz || sna.type.txz || sna.type.tyz) izs = -1;
if ( !ISODD(bnd.lef)) hdr.gx = 1000*(mod.x0 - bnd.ntap*mod.dx);
if ( !ISODD(bnd.fro)) hdr.gy = 1000*(mod.y0 - bnd.ntap*mod.dy);
}
if (sna.type.vx) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = vx[iy*n1*n2+ix2*n1+iz];
}
traceWrite(&hdr, snap, (int)sna.nz, fpvx);
}
if (sna.type.vy) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = vy[iy2*n1*n2+ix*n1+iz];
}
traceWrite(&hdr, snap, (int)sna.nz, fpvy);
}
if (sna.type.vz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = vz[iy*n1*n2+ix*n1+iz+1];
}
traceWrite(&hdr, snap, (int)sna.nz, fpvz);
}
if (sna.type.p) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = tzz[iy*n1*n2+ix*n1+iz];
}
traceWrite(&hdr, snap, (int)sna.nz, fpp);
}
if (sna.type.tzz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = tzz[iy*n1*n2+ix*n1+iz];
}
traceWrite(&hdr, snap, (int)sna.nz, fptzz);
}
if (sna.type.tyy) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = tyy[iy*n1*n2+ix*n1+iz];
}
traceWrite(&hdr, snap, (int)sna.nz, fptyy);
}
if (sna.type.txx) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = txx[iy*n1*n2+ix*n1+iz];
}
traceWrite(&hdr, snap, (int)sna.nz, fptxx);
}
if (sna.type.txz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = txz[iy*n1*n2+ix2*n1+iz+1];
}
traceWrite(&hdr, snap, (int)sna.nz, fptxz);
}
if (sna.type.txy) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = txy[iy2*n1*n2+ix2*n1+iz];
}
traceWrite(&hdr, snap, (int)sna.nz, fptxy);
}
if (sna.type.tyz) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = tyz[iy2*n1*n2+ix*n1+iz+1];
}
traceWrite(&hdr, snap, (int)sna.nz, fptyz);
}
/* calculate divergence of velocity field */
if (sna.type.pp) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = sdx*((vx[iy*n1*n2+(ix+1)*n1+iz]-vx[iy*n1*n2+ix*n1+iz])+
(vy[(iy+1)*n1*n2+ix*n1+iz]-vy[iy*n1*n2+ix*n1+iz])+
(vz[iy*n1*n2+ix*n1+iz+1]-vz[iy*n1*n2+ix*n1+iz]));
}
traceWrite(&hdr, snap, (int)sna.nz, fppp);
}
/* calculate rotation of velocity field */
if (sna.type.ss) {
for (iz=izs, j=0; iz<=ize; iz+=sna.skipdz, j++) {
snap[j] = sdx*((vx[iy*n1*n2+ix*n1+iz]-vx[(iy-1)*n1*n2+ix*n1+iz-1])-
(vy[iy*n1*n2+ix*n1+iz]-vy[iy*n1*n2+(ix-1)*n1+iz-1])-
(vz[iy*n1*n2+ix*n1+iz]-vz[(iy-1)*n1*n2+(ix-1)*n1+iz]));
}
traceWrite(&hdr, snap, (int)sna.nz, fpss);
}
}
}
if (sna.type.vx) fclose(fpvx);
if (sna.type.vy) fclose(fpvy);
if (sna.type.vz) fclose(fpvz);
if (sna.type.p) fclose(fpp);
if (sna.type.txx) fclose(fptxx);
if (sna.type.tyy) fclose(fptyy);
if (sna.type.tzz) fclose(fptzz);
if (sna.type.txz) fclose(fptxz);
if (sna.type.tyz) fclose(fptyz);
if (sna.type.txy) fclose(fptxy);
if (sna.type.pp) fclose(fppp);
if (sna.type.ss) fclose(fpss);
free(snap);
}
return 0;
}
fdelmodc3D/writeSrcRecPos3D.c 0 → 100644
+ 160
− 0
View file @ 44c3dd2d
#include<stdlib.h>
#include<stdio.h>
#include<math.h>
#include<assert.h>
#include"par.h"
#include"fdelmodc3D.h"
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
/**
* Writes the source and receiver positions into a gridded file,
* which has the same size as the input gridded model files.
* Source positions have a value +1 and receivers -1.
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
long writesufile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2);
long writeSrcRecPos3D(modPar *mod, recPar *rec, srcPar *src, shotPar *shot)
{
FILE *fp;
float *dum, sub_x0, sub_y0, sub_z0, dx, dy, dz;
long is, nx, ny, nz, is0, ish, ix, iy, iz, ndot, idx, idy, idz;
char tmpname[1024];
ndot = 2;
nx = mod->nx;
ny = mod->ny;
nz = mod->nz;
dx = mod->dx;
dy = mod->dy;
dz = mod->dz;
sub_x0 = mod->x0;
sub_y0 = mod->y0;
sub_z0 = mod->z0;
/* write velocity field with positions of the sources */
dum = (float *)calloc(nx*ny*nz, sizeof(float));
vmess("Positions: shot=%li src=%li rec=%li", shot->n, src->n, rec->n);
/* source positions for random shots */
if (src->random) {
sprintf(tmpname,"SrcPositions%li.txt",src->n);
fp = fopen(tmpname, "w+");
for (is=0; is<src->n; is++) {
for (idy=0; idy<=ndot; idy++) {
for (idx=0; idx<=ndot; idx++) {
for (idz=0; idz<=ndot; idz++) {
dum[(MAX(0,src->y[is]-idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
dum[(MAX(0,src->y[is]-idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
dum[(MAX(0,src->y[is]-idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
dum[(MAX(0,src->y[is]-idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MAX(0,src->x[is]-idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MIN(nz-1,src->z[is]+idz)] = 1.0;
dum[(MIN(ny-1,src->y[is]+idy))*nz*nx+(MIN(nx-1,src->x[is]+idx))*nz+MAX(0,src->z[is]-idz)] = 1.0;
}
}
}
fprintf(fp, "%f %f %f\n", src->z[is]*dz+sub_z0, src->y[is]*dy+sub_y0, src->x[is]*dx+sub_x0);
}
fclose(fp);
}
/* source positions for single shot sources with plane waves */
else if (src->plane) {
is0 = -1*floor((src->n-1)/2);
sprintf(tmpname,"SrcPositions%li.txt",shot->n);
fp = fopen(tmpname, "w+");
for (ish=0; ish<shot->n; ish++) {
for (is=0; is<src->n; is++) {
ix = shot->x[ish] + 1 + is0 + is;
iy = shot->y[ish] + 1 + is0 + is;
iz = shot->z[ish] + 1;
dum[iy*nx*nz+ix*nz+iz] = 1.0;
dum[(MAX(0,iy-1))*nx*nz+ix*nz+iz] = 1.0;
dum[(MIN(ny-1,iy+1))*nx*nz+ix*nz+iz] = 1.0;
dum[iy*nx*nz+(MAX(0,ix-1))*nz+iz] = 1.0;
dum[iy*nx*nz+(MIN(nx-1,ix+1))*nz+iz] = 1.0;
dum[iy*nx*nz+ix*nz+MAX(0,iz-1)] = 1.0;
dum[iy*nx*nz+ix*nz+MIN(nz-1,iz+1)] = 1.0;
fprintf(fp, "(%f, %f, %f)\n", ix*dx+sub_x0, iy*dy+sub_y0, iz*dz+sub_z0);
}
}
fclose(fp);
}
else if (src->multiwav) {
/* source positions for single shot sources with multiple wavelets */
sprintf(tmpname,"SrcPositions%li.txt",shot->n);
fp = fopen(tmpname, "w+");
for (ish=0; ish<shot->n; ish++) {
for (is=0; is<src->n; is++) {
ix = src->x[is];
iy = src->y[is];
iz = src->z[is];
dum[iy*nx*nz+ix*nz+iz] = 1.0;
dum[(MAX(0,iy-1))*nx*nz+ix*nz+iz] = 1.0;
dum[(MIN(ny-1,iy+1))*nx*nz+ix*nz+iz] = 1.0;
dum[iy*nx*nz+(MAX(0,ix-1))*nz+iz] = 1.0;
dum[iy*nx*nz+(MIN(nx-1,ix+1))*nz+iz] = 1.0;
dum[iy*nx*nz+ix*nz+MAX(0,iz-1)] = 1.0;
dum[iy*nx*nz+ix*nz+MIN(nz-1,iz+1)] = 1.0;
fprintf(fp, "(%f, %f, %f)\n", ix*dx+sub_x0, iy*dy+sub_y0, iz*dz+sub_z0);
}
}
fclose(fp);
}
else {
sprintf(tmpname,"SrcPositions%li.txt",shot->n);
fp = fopen(tmpname, "w+");
for (is=0; is<shot->n; is++) {
for (idy=0; idy<=ndot; idy++) {
for (idx=0; idx<=ndot; idx++) {
for (idz=0; idz<=ndot; idz++) {
dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
dum[(MAX(0,shot->y[is]-idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MAX(0,shot->x[is]-idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MIN(nz-1,shot->z[is]+idz)] = 1.0;
dum[(MIN(ny-1,shot->y[is]+idy))*nz*nx+(MIN(nx-1,shot->x[is]+idx))*nz+MAX(0,shot->z[is]-idz)] = 1.0;
}
}
}
fprintf(fp, "%f %f %f\n", shot->z[is]*dz+sub_z0, shot->y[is]*dy+sub_y0, shot->x[is]*dx+sub_x0);
}
fclose(fp);
}
/* receiver positions */
sprintf(tmpname,"RcvPositions%li.txt",rec->n);
fp = fopen(tmpname, "w+");
for (is=0; is<rec->n; is++) {
dum[rec->y[is]*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
dum[(MAX(0,rec->y[is]-1))*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
dum[(MIN(ny-1,rec->y[is]+1))*nx*nz+rec->x[is]*nz+rec->z[is]] = -1.0;
dum[rec->y[is]*nx*nz+(MAX(0,rec->x[is]-1))*nz+rec->z[is]] = -1.0;
dum[rec->y[is]*nx*nz+(MIN(nx-1,rec->x[is]+1))*nz+rec->z[is]] = -1.0;
dum[rec->y[is]*nx*nz+rec->x[is]*nz+MAX(0,rec->z[is]-1)] = -1.0;
dum[rec->y[is]*nx*nz+rec->x[is]*nz+MIN(nz-1,rec->z[is]+1)] = -1.0;
if (rec->int_vx==3) {
fprintf(fp, "(%f, %f, %f)\n", rec->xr[is]*dx+sub_x0, rec->yr[is]*dy+sub_y0, rec->zr[is]*dz+sub_z0);
}
else {
fprintf(fp, "(%f, %f, %f)\n", rec->x[is]*dx+sub_x0, rec->y[is]*dy+sub_y0, rec->z[is]*dz+sub_z0);
}
}
fclose(fp);
writesufile3D("SrcRecPositions.su", dum, nz, nx*ny, sub_z0, sub_x0, dz, dx);
free(dum);
return 0;
}
fdelmodc3D/writesufile3D.c 0 → 100644
+ 162
− 0
View file @ 44c3dd2d
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
#include "par.h"
#include "fdelmodc3D.h"
#include "SUsegy.h"
#include "segy.h"
/**
* Writes an 2D array to a SU file
*
* AUTHOR:
* Jan Thorbecke (janth@xs4all.nl)
* The Netherlands
**/
#define TRCBYTES 240
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#define ISODD(n) ((n) & 01)
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
long writesufile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2)
{
FILE *file_out;
size_t nwrite, itrace;
long ns;
segy *hdr;
/* Read in parameters */
if (n1 > USHRT_MAX) {
vwarn("Output file %s: number of samples is truncated from %li to USHRT_MAX.", filename, n1);
}
ns = MIN(n1,USHRT_MAX);
file_out = fopen( filename, "w+" );
assert( file_out );
hdr = (segy *)calloc(1,TRCBYTES);
hdr->ns = ns;
hdr->dt = NINT(1000000*(d1));
hdr->d1 = d1;
hdr->d2 = d2;
hdr->f1 = f1;
hdr->f2 = f2;
hdr->fldr = 1;
hdr->trwf = n2;
for (itrace=0; itrace<n2; itrace++) {
hdr->tracl = itrace+1;
nwrite = fwrite( hdr, 1, TRCBYTES, file_out );
assert (nwrite == TRCBYTES);
nwrite = fwrite( &data[itrace*n1], sizeof(float), ns, file_out );
assert (nwrite == ns);
}
fclose(file_out);
free(hdr);
return 0;
}
/**
* Writes an 2D array to a SU file
* special routine for src_nwav array which has a different number of samples for each shot
*
**/
long writesufilesrcnwav3D(char *filename, float **src_nwav, wavPar wav, long n1, long n2, float f1, float f2, float d1, float d2)
{
FILE *file_out;
size_t nwrite, itrace;
float *trace;
long ns;
segy *hdr;
/* Read in parameters */
if (n1 > USHRT_MAX) {
vwarn("Output file %s: number of samples is truncated from %li to USHRT_MAX.", filename, n1);
}
ns = MIN(n1,USHRT_MAX);
file_out = fopen( filename, "w+" );
assert( file_out );
trace = (float *)malloc(n1*sizeof(float));
hdr = (segy *)calloc(1,TRCBYTES);
hdr->ns = ns;
hdr->dt = NINT(1000000*(d1));
hdr->d1 = d1;
hdr->d2 = d2;
hdr->f1 = f1;
hdr->f2 = f2;
hdr->fldr = 1;
hdr->trwf = n2;
for (itrace=0; itrace<n2; itrace++) {
hdr->tracl = itrace+1;
nwrite = fwrite( hdr, 1, TRCBYTES, file_out );
assert (nwrite == TRCBYTES);
memset(trace, 0, n1*sizeof(float));
memcpy(trace, &src_nwav[itrace][0], wav.nsamp[itrace]*sizeof(float));
nwrite = fwrite( &trace[0], sizeof(float), ns, file_out );
assert (nwrite == ns);
}
fclose(file_out);
free(hdr);
free(trace);
return 0;
}
/**
* Writes an 2D array to a SU file
* special routine which used segyhdrs which have ns defined as integer (32 bit)
* to handle more than 2^16 samples per trace.
*
**/
long writeSUfile3D(char *filename, float *data, long n1, long n2, float f1, float f2, float d1, float d2)
{
FILE *file_out;
size_t nwrite, itrace;
SUsegy *SUhdr;
char *ptr;
/* Read in parameters */
ptr = strstr(filename, " ");
*ptr = '\0';
file_out = fopen( filename, "w+" );
assert( file_out );
SUhdr = (SUsegy *)calloc(1,TRCBYTES);
SUhdr->ns = n1;
SUhdr->dt = NINT(1000000*(d1));
SUhdr->d1 = d1;
SUhdr->d2 = d2;
SUhdr->f1 = f1;
SUhdr->f2 = f2;
SUhdr->fldr = 1;
SUhdr->trwf = n2;
for (itrace=0; itrace<n2; itrace++) {
SUhdr->tracl = itrace+1;
nwrite = fwrite( SUhdr, 1, TRCBYTES, file_out );
assert (nwrite == TRCBYTES);
nwrite = fwrite( &data[itrace*n1], sizeof(float), n1, file_out );
assert (nwrite == n1);
}
fclose(file_out);
free(SUhdr);
return 0;
}
marchenko/Makefile
+ 24
− 4
View file @ 44c3dd2d
This diff is collapsed.
marchenko/marchenko_primaries.c 0 → 100644
+ 998
− 0
View file @ 44c3dd2d
This diff is collapsed.
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