Commit 2a316380 authored by niklas.thiel's avatar niklas.thiel

wavefield seperation added

parent 94121d66
......@@ -65,6 +65,13 @@
"XREC2, YREC2" : "93.0 , 0.2",
"NGEOPH" : "80",
"Wavefield Separation" : "comment",
"WAVESEP" : "0",
"VEL" : "1500.0",
"DENS" : "1000.0",
"ANGTAPI" : "0.0",
"ANGTAPO" : "70.0",
"Seismograms" : "comment",
"NDT" : "1",
"SEIS_FORMAT" : "1",
......
......@@ -65,6 +65,14 @@
"XREC2, YREC2" : "93.0 , 0.2",
"NGEOPH" : "80",
"Wavefield Separation" : "comment",
"WAVESEP" : "0",
"VEL" : "1500",
"DENS" : "1000",
"ANGTAPI" : "0",
"ANGTAPO" : "70",
"Seismograms" : "comment",
"NDT" : "1",
"SEIS_FORMAT" : "1",
......
......@@ -699,21 +699,17 @@ int main(int argc, char **argv){
fulldata_curl = matrix(1,ntr_glob,1,NT);
break;
case 5 : /* everything except curl and div*/
switch (WAVETYPE) {
case 1:
if (WAVETYPE==1) {
fulldata_vx = matrix(1,ntr_glob,1,NT);
fulldata_vy = matrix(1,ntr_glob,1,NT);
break;
case 2:
}
if (WAVETYPE==2) {
fulldata_vz = matrix(1,ntr_glob,1,NT);
break;
case 3:
}
if (WAVETYPE==3) {
fulldata_vx = matrix(1,ntr_glob,1,NT);
fulldata_vy = matrix(1,ntr_glob,1,NT);
fulldata_vz = matrix(1,ntr_glob,1,NT);
break;
}
fulldata_p = matrix(1,ntr_glob,1,NT);
break;
......@@ -1442,7 +1438,18 @@ int main(int argc, char **argv){
}
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
/* wavefield separation*/
if (WAVESEP==1) {
if (MYID == 0) {
fprintf(FP,"\n Testposition 0 \n");
pup(fulldata_p, fulldata_vz, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
}
/* overwrite sectionp with pup data */
inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter);
}
break;
} /* end of switch (SEISMO) */
......@@ -1985,7 +1992,15 @@ int main(int argc, char **argv){
catseis(sectionvz, fulldata_vz, recswitch, ntr_glob, MPI_COMM_WORLD);
}
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
if (MYID==0) saveseis_glob(FP,fulldata_vx,fulldata_vy,fulldata_vz,fulldata_p,fulldata_curl,fulldata_div,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
/* wavefield separation */
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n Testposition out pup");
if (MYID == 0) pup(fulldata_p, fulldata_vy, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,1);
/* overwrite sectionp with pup data */
inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter);
}
break;
} /* end of switch (SEISMO) */
......@@ -3485,6 +3500,15 @@ int main(int argc, char **argv){
}
/* wavefield separation */
if (WAVESEP==1) {
catseis(sectionvz, fulldata_vz, recswitch, ntr_glob, MPI_COMM_WORLD);
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
if (MYID == 0) pup(fulldata_p, fulldata_vz, FP, ntr_glob,recpos,recpos_loc,ntr_glob,srcpos,ishot,ns,iter,2);
/* overwrite sectionp with pup data */
inseis(fprec,ishot,sectionp,ntr_glob,ns,3,iter);
}
/*-------------------------------------------------------------------------------*/
/*------------------ end loop over timesteps (forward model) step length -------*/
/*-------------------------------------------------------------------------------*/
......
......@@ -163,7 +163,10 @@ IFOS2D= \
readmod_viscac.c \
time_window_glob.c \
create_trkill_table.c \
filter_frequencies.c
filter_frequencies.c \
pup.c \
fft2_filt.c \
fft2d.c
SNAPMERGE_OBJ = $(SNAPMERGE_SCR:%.c=%.o)
......
......@@ -117,6 +117,9 @@ void exchange_par(void){
extern int WOLFE_TRY_OLD_STEPLENGTH;
extern float WOLFE_C1_SL;
extern float WOLFE_C2_SL;
extern int WAVESEP;
extern float VEL, DENS, ANGTAPI, ANGTAPO;
/* definition of local variables */
......@@ -218,6 +221,12 @@ void exchange_par(void){
fdum[68]=TRKILL_OFFSET_LOWER;
fdum[69]=TRKILL_OFFSET_UPPER;
fdum[70]=VEL;
fdum[71]=DENS;
fdum[72]=ANGTAPI;
fdum[73]=ANGTAPO;
/***********/
/* Integer */
/***********/
......@@ -371,6 +380,10 @@ void exchange_par(void){
idum[114]=TRKILL_OFFSET;
idum[115]=TRKILL_STF_OFFSET;
idum[116]=WAVESEP;
} /** if (MYID == 0) **/
......@@ -498,6 +511,12 @@ void exchange_par(void){
TRKILL_STF_OFFSET_UPPER=fdum[67];
TRKILL_OFFSET_LOWER=fdum[68];
TRKILL_OFFSET_UPPER=fdum[69];
VEL=fdum[70];
DENS=fdum[71];
ANGTAPI=fdum[72];
ANGTAPO=fdum[73];
/***********/
/* Integer */
......@@ -655,6 +674,8 @@ void exchange_par(void){
TRKILL_OFFSET=idum[114];
TRKILL_STF_OFFSET=idum[115];
WAVESEP=idum[116];
MPI_Bcast(&FL[1],L,MPI_FLOAT,0,MPI_COMM_WORLD);
......
......@@ -10,6 +10,7 @@
#include <stddef.h>
#include <string.h>
#include <time.h>
#include <unistd.h>
#include <mpi.h>
#define iround(x) ((int)(floor)(x+0.5))
......@@ -24,8 +25,32 @@
#define REQUEST_COUNT 4
#define WORKFLOW_MAX_VAR 12
/* declaration of functions */
/***************/
/* FFT2D stuff */
/***************/
typedef struct {float re; float im;} COMPLEX;
typedef struct {double re; double im;} DCOMPLEX;
#ifndef ERROR
#define ERROR -1
#define NO_ERROR 0
#endif
#define FFT_FORWARD 0
#define FFT_INVERSE 1
#define TWOPI 6.2831853071795865 /* 2.0 * PI */
#define HALFPI 1.5707963267948966 /* PI / 2.0 */
#define PI8 0.392699081698724 /* PI / 8.0 */
#define RT2 1.4142135623731 /* sqrt(2.0) */
#define IRT2 0.707106781186548 /* 1.0/sqrt(2.0) */
extern int forward_fft2f(COMPLEX *array, int rows, int cols); /* Perform forward 2D transform on a COMPLEX array. */
extern int inverse_fft2f(COMPLEX *array, int rows, int cols); /* Perform inverse 2D transform on a COMPLEX array. */
/****************************/
/* declaration of functions */
/****************************/
void window_cos(float **win, int npad, int nsrc, float it1, float it2, float it3, float it4);
......@@ -125,6 +150,8 @@ void FFT_filt(float ** data, float freqshift, int ntr, int ns,int method);
void FFT(int isign, unsigned long nlog, float *re, float *im); /* NR version*/
void fft2(float **array, float **arrayim, int NYG, int NXG, int dir); /* 2d fft filtering*/
float *filter_frequencies(int *nfrq);
float *holbergcoeff(void);
......@@ -201,6 +228,7 @@ void psource(int nt, float ** sxx, float ** syy, float ** sp,
void psource_rsg(int nt, float ** sxx, float ** syy,
float ** srcpos_loc, float ** signals, int nsrc);
void pup(float **data_p, float **data_vz, FILE *fp, int ntr_glob, int **recpos, int **recpos_loc,int ntr, float ** srcpos, int ishot, int ns, int iter, int sw );
float *rd_sour(int *nts,FILE* fp_source);
......@@ -464,6 +492,7 @@ void err(char err_text[]);
void warning(char warn_text[]);
double maximum(float **a, int nx, int ny);
float *vector(int nl, int nh);
COMPLEX *cplxvector(int nstart, int nend);
int *ivector(int nl, int nh);
double *dvector(int nl, int nh);
float **fmatrix(int nrl, int nrh, int ncl, int nch);
......@@ -473,6 +502,7 @@ float **matrix(int nrl, int nrh, int ncl, int nch);
int **imatrix(int nrl, int nrh, int ncl, int nch);
float ***f3tensor(int nrl, int nrh, int ncl, int nch,int ndl, int ndh);
void free_vector(float *v, int nl, int nh);
void free_cplxvector(COMPLEX *v, int nstart, int nend);
void free_dvector(double *v, int nl, int nh);
void free_ivector(int *v, int nl, int nh);
void free_matrix(float **m, int nrl, int nrh, int ncl, int nch);
......
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2016 For the list of authors, see file AUTHORS.
*
* This file is part of IFOS.
*
* IFOS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* IFOS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* 2D-FFT preparation
* ----------------------------------------------------------------------*/
#include "fd.h"
void fft2(float **array, float **arrayim, int NYG, int NXG, int dir) {
extern FILE *FP;
COMPLEX *C;
int j, i, k, nx, ny, nxy;
float **RE, **IM;
fprintf(FP,"\n Testposition fft2_filt 1\n ");
usleep(900000);
/****************************************************/
/* recompute the dimensions to become values of 2^m */
nx = (int)(pow(2,(ceil(log(NXG)/log(2)))));
ny = (int)(pow(2,(ceil(log(NYG)/log(2)))));
nxy = max(nx,ny);
nxy *= 2;
fprintf(FP,"\n Testposition fft2_filt 1, nxy: %i\n ",nxy);
usleep(900000);
RE = matrix(1, nxy, 1, nxy);
IM = matrix(1, nxy, 1, nxy);
C = cplxvector(1, nxy*nxy);
/* copy array data to the temp-array */
for (i=1;i<=NXG;i++)
for (j=1;j<=NYG;j++){
RE[j][i] = array[j][i];
IM[j][i] = arrayim[j][i];
}
fprintf(FP,"\n Testposition fft2_filt 3\n ");
usleep(900000);
if (dir==1) {
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
C[k].re = RE[j][i];
C[k].im = 0.0;
k++;
}
forward_fft2f(C, nxy, nxy);
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
RE[j][i] = C[k].re;
IM[j][i] = C[k].im;
k++;
}
} else {
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
C[k].re = RE[j][i];
C[k].im = IM[j][i];
k++;
}
inverse_fft2f(C, nxy, nxy);
k = 1;
for (i=1;i<=nxy;i++)
for (j=1;j<=nxy;j++) {
RE[j][i] = C[k].re;
IM[j][i] = 0.0;
k++;
}
}
/* write result back into array*/
for (i=1;i<=NXG;i++)
for (j=1;j<=NYG;j++){
array[j][i]=RE[j][i];
arrayim[j][i]=IM[j][i];
}
free_cplxvector(C, 1, nxy*nxy);
free_matrix(RE, 1,nxy,1,nxy);
free_matrix(IM, 1,nxy,1,nxy);
}
/********************************************************/
/* 2D-FFTSHIFT */
/********************************************************/
// void fft2shift(float **RE, float **IM, int ny, int nx) {
//
// int j, i;
// float **A;
//
//
// A = matrix(1, ny, 1, nx);
//
// /* shift the real part */
// for (i=1;i<=nx;i++)
// for (j=1;j<=ny;j++) A[j][i] = RE[j][i];
// for (i=1;i<=nx/2;i++)
// for (j=1;j<=ny/2;j++) {
// RE[j][i] = A[j+ny/2][i+nx/2];
// RE[j+ny/2][i+nx/2] = A[j][i];
// RE[j+ny/2][i] = A[j][i+nx/2];
// RE[j][i+nx/2] = A[j+ny/2][i];
// }
//
// /* shift the imaginary part */
// for (i=1;i<=nx;i++)
// for (j=1;j<=ny;j++) A[j][i] = IM[j][i];
// for (i=1;i<=nx/2;i++)
// for (j=1;j<=ny/2;j++) {
// IM[j][i] = A[j+ny/2][i+nx/2];
// IM[j+ny/2][i+nx/2] = A[j][i];
// IM[j+ny/2][i] = A[j][i+nx/2];
// IM[j][i+nx/2] = A[j+ny/2][i];
// }
//
// free_matrix(A, 1, ny, 1, nx);
// }
\ No newline at end of file
/*
------------------------------------------------------------------------------
*
* Copyright (C) 2013, see file AUTHORS for list of authors/contributors
*
* This file is part of PROTEUS.
*
* PROTEUS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
*
* PROTEUS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with PROTEUS. See file COPYING and/or
* <http://www.gnu.org/licenses/gpl-2.0.html>.
*
------------------------------------------------------------------------------
*/
/**************************************************************
* kube-gustavson-fft.c
*
* Forward and inverse discrete 2D Fourier transforms.
*
* Originally a FORTRAN implementation published in:
* Programming for Digital Signal Processing, IEEE Press 1979,
* Section 1, by G. D. Bergland and M. T. Dolan.
* Ported long ago from Fortran to old-fashioned Standard C by
* someone who failed to put his or her name in the source.
* Ported from Standard C to ANSI C by Stefan Gustavson
* (stegu@itn.liu.se) 2003-10-20.
*
* This is a pretty fast FFT algorithm. Not super-optimized
* lightning-fast, but very good considering its age and
* relative simplicity. There are several places in the code
* where clock cycles could be saved, for example by getting rid
* of some copying of data back and forth, but the savings
* would not be that great. If you want optimally fast FFTs,
* consider the excellent FFTW package. (http://www.fftw.org)
*
* This software is in the public domain.
*
**************************************************************
*
* int forward_fft2f(COMPLEX *array, int rows, int cols)
* int inverse_fft2f(COMPLEX *array, int rows, int cols)
* int forward_fft2d(DCOMPLEX *array, int rows, int cols)
* int inverse_fft2d(DCOMPLEX *array, int rows, int cols)
*
* These functions compute the forward and inverse DFT's, respectively,
* of a single-precision COMPLEX or DCOMPLEX array by means of an
* FFT algorithm.
*
* The result is a COMPLEX/DCOMPLEX array of the same size, returned
* in the same space as the input array. That is, the original array
* is overwritten and destroyed.
*
* Rows and columns must each be an integral power of 2.
*
* These routines return integer value ERROR if an error was detected,
* NO_ERROR otherwise.
*
* This particular implementation of the DFT uses the transform pair
* defined as follows:
* Let there be two complex arrays each with n rows and m columns.
* Index them as
* f(x,y): 0 <= x <= m - 1, 0 <= y <= n - 1
* F(u,v): -m/2 <= u <= m/2 - 1, -n/2 <= v <= n/2 - 1
*
* Then the forward and inverse transforms are related as follows.
* Forward:
* F(u,v) = \sum_{x=0}^{m-1} \sum_{y=0}^{n-1}
* f(x,y) \exp{-2\pi i (ux/m + vy/n)}
*
* Inverse:
* f(x,y) = 1/(mn) \sum_{u=-m/2}^{m/2-1} \sum_{v=-n/2}^{n/2-1}
* F(u,v) \exp{2\pi i (ux/m + vy/n)}
*
* Therefore, the transforms have these properties:
* 1. \sum_x \sum_y f(x,y) = F(0,0)
* 2. m n \sum_x \sum_y |f(x,y)|^2 = \sum_u \sum_v |F(u,v)|^2
*
*/
#include "fd.h"
/*
* These somewhat awkward macros move data between the input/output array
* and stageBuff, while also performing any conversions float<->double.
*/
#define LoadRow(buffer,row,cols) if (1){\
register int j,k;\
for (j = row*cols, k = 0 ; k < cols ; j++, k++){\
stageBuff[k].re = buffer[j].re;\
stageBuff[k].im = buffer[j].im;\
}\
} else {}
#define StoreRow(buffer,row,cols) if (1){\
register int j,k;\
for (j = row*cols, k = 0 ; k < cols ; j++, k++){\
buffer[j].re = stageBuff[k].re;\
buffer[j].im = stageBuff[k].im;\
}\
} else {}
#define LoadCol(buffer,col,rows,cols) if (1){\
register int j,k;\
for (j = i,k = 0 ; k < rows ; j+=cols, k++){\
stageBuff[k].re = buffer[j].re;\
stageBuff[k].im = buffer[j].im;\
}\
} else {}
#define StoreCol(buffer,col,rows,cols) if (1){\
register int j,k;\
for (j = i,k = 0 ; k < rows ; j+=cols, k++){\
buffer[j].re = stageBuff[k].re;\
buffer[j].im = stageBuff[k].im;\
}\
} else {}
/* Do something useful with an error message */
#define handle_error(msg) fprintf(stderr,msg)
DCOMPLEX *stageBuff; /* buffer to hold a row or column at a time */
COMPLEX *bigBuff; /* a pointer to a float input array */
DCOMPLEX *bigBuffd; /* a pointer to a double input array */
/* Allocate space for stageBuff */
int allocateBuffer(int size)
{
stageBuff = (DCOMPLEX*) malloc(size*sizeof(DCOMPLEX));
if(stageBuff==(DCOMPLEX*)NULL)
{
handle_error("Insufficient storage for fft buffers");
return(ERROR);
}
else return(NO_ERROR);
}
/* Free space for stageBuff */
void freeBuffer(void)
{
free((char*)stageBuff);
}
/* See if exactly one bit is set in integer argument */
int power_of_2(int n)
{
int bits=0;
while(n) {
bits += n & 1;
n >>= 1;
}
return(bits==1);
}
/* Get binary log of integer argument - exact if n is a power of 2 */
int fastlog2(int n)
{
int log = -1;
while(n) {
log++;
n >>= 1;
}
return(log);
}
/* Radix-2 iteration subroutine */
void R2TX(int nthpo, DCOMPLEX *c0, DCOMPLEX *c1)
{
int k,kk;
double *cr0, *ci0, *cr1, *ci1, r1, fi1;
cr0 = &(c0[0].re);
ci0 = &(c0[0].im);
cr1 = &(c1[0].re);
ci1 = &(c1[0].im);
for(k=0; k<nthpo; k+=2)
{
kk = k*2;
r1 = cr0[kk] + cr1[kk];
cr1[kk] = cr0[kk] - cr1[kk];
cr0[kk] = r1;
fi1 = ci0[kk] + ci1[kk];
ci1[kk] = ci0[kk] - ci1[kk];
ci0[kk] = fi1;
}
}
/* Radix-4 iteration subroutine */
void R4TX(int nthpo, DCOMPLEX *c0, DCOMPLEX *c1,
DCOMPLEX *c2, DCOMPLEX *c3)
{
int k,kk;
double *cr0, *ci0, *cr1, *ci1, *cr2, *ci2, *cr3, *ci3;
double r1, r2, r3, r4, i1, i2, i3, i4;
cr0 = &(c0[0].re);
cr1 = &(c1[0].re);
cr2 = &(c2[0].re);
cr3 = &(c3[0].re);
ci0 = &(c0[0].im);
ci1 = &(c1[0].im);
ci2 = &(c2[0].im);
ci3 = &(c3[0].im);
for(k=1;k<=nthpo;k+=4)
{
kk = (k-1)*2; /* real and imag parts alternate */
r1 = cr0[kk] + cr2[kk];
r2 = cr0[kk] - cr2[kk];
r3 = cr1[kk] + cr3[kk];
r4 = cr1[kk] - cr3[kk];
i1 = ci0[kk] + ci2[kk];
i2 = ci0[kk] - ci2[kk];
i3 = ci1[kk] + ci3[kk];
i4 = ci1[kk] - ci3[kk];
cr0[kk] = r1 + r3;
ci0[kk] = i1 + i3;
cr1[kk] = r1 - r3;
ci1[kk] = i1 - i3;
cr2[kk] = r2 - i4;
ci2[kk] = i2 + r4;
cr3[kk] = r2 + i4;
ci3[kk] = i2 - r4;
}
}
/* Radix-8 iteration subroutine */
void R8TX(int nxtlt, int nthpo, int length,
DCOMPLEX *cc0, DCOMPLEX *cc1, DCOMPLEX *cc2, DCOMPLEX *cc3,
DCOMPLEX *cc4, DCOMPLEX *cc5, DCOMPLEX *cc6, DCOMPLEX *cc7)
{
int j,k,kk;
double scale, arg, tr, ti;
double c1, c2, c3, c4, c5, c6, c7;
double s1, s2, s3, s4, s5, s6, s7;
double ar0, ar1, ar2, ar3, ar4, ar5, ar6, ar7;
double ai0, ai1, ai2, ai3, ai4, ai5, ai6, ai7;
double br0, br1, br2, br3, br4, br5, br6, br7;
double bi0, bi1, bi2, bi3, bi4, bi5, bi6, bi7;