...
 
Commits (21)
......@@ -350,6 +350,23 @@ Default values are:
If SEISMO$>$0 seismograms recorded at the receiver positions are written to the corresponding output files. The sampling rate of the seismograms is NDT*DT seconds. In case of a small time step interval and a high number of time steps, it might be useful to choose a high NDT in order to avoid a unnecessary detailed sampling of the seismograms and consequently large files of seismogram data. Possible output formats of the seismograms are SU, ASCII and BINARY. It is recommended to use SU format for saving the seismograms. The main advantage of this format is that the time step interval (NDT*DT) and the acquisition geometry (shot and receiver locations) are stored in the corresponding SU header words. Also additional header words like offset are set by IFOS2D. This format thus facilitates a further visualization and processing of the synthetic seismograms. Note, however, that SU cannot handle sampling rates smaller than 1.0e-6 seconds and the number of samples is limited to about 32.000. In such cases, you should increase the sampling rate by increasing NDT. If this is impossible (for example because the Nyquist criterion is violated) you must choose a different output format (ASCII or binary). File endings will be added to SEIS\_FILE automatically.
\section{Wavefield separation}
\label{wavefield_sep}
{\color{blue}{\begin{verbatim}
"Wavefield separation" : "comment",
"WAVESEP" : "1",
"VEL" : "1500",
"DENS" : "1000",
"ANGTAPI" : "0",
"ANGTAPO" : "70",
\end{verbatim}}}
{\color{red}{\begin{verbatim}
Default values are:
WAVESEP=0
\end{verbatim}}}
If WAVESEP$=$1 a wavefield separation is performed after every forward modeling. The recorded wavefield is separated into an upgoing and downgoing wavefield by using the recorded p-data and vz-data (\cite{amundsen:93}). Therefore, SEISMO$=$4 or SEISMO$=$5 are mandatory. For the following inversion only the upgpoing wavefield is used. For the calculation of the filter the velocity (VEL in m/s) and the density (DENS in kg/m$^3$) are needed, as well as the angles of the taper of the filter (ANGTAPI, ANGTAPO in degree with $0^\circ$ $<=$ ANGTAPI, ANGTAPO $<=$ $90^\circ$ and ANGTAPI $<=$ ANGTAPO). ANGTAPO describes the maximum used dipping angle. ANGTAPI defines the start of the taper function of the filter.
\section{Q-approximation}
\label{q_approx}
......
......@@ -3079,3 +3079,13 @@ author = {Choi, Y. and Alkhalifah, T.},
pages = {748--758},
year ={2012}
}
@article{amundsen:93,
title={Wavenumber-based filtering of marine point-source data},
author={Amundsen, Lasse},
journal=Geo,
volume={58},
number={9},
pages={1335--1348},
year={1993},
publisher={Society of Exploration Geophysicists}
}
......@@ -58,6 +58,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",
......
......@@ -58,6 +58,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",
......
......@@ -156,6 +156,8 @@ int main(int argc, char **argv){
float *F_LOW_PASS_EXT=NULL;
int nfrq=0;
int FREQ_NR=1;
char pf[STRING_SIZE];
float JOINT_EQUAL_PSV=0.0, JOINT_EQUAL_SH=0.0;
float JOINT_EQUAL_PSV_all=0.0, JOINT_EQUAL_SH_all=0.0;
......@@ -814,6 +816,18 @@ int main(int argc, char **argv){
MPI_Barrier(MPI_COMM_WORLD);
/* MPI split for processors used for pup parallelisation */
int myid_pup, pupid=0, no_myid;
MPI_Comm MPI_COMM_PUP;
no_myid=9; /* must be power of 2 +1 */
if (MYID<=no_myid-1) pupid = 1;
else pupid = 0;
MPI_Comm_split(MPI_COMM_WORLD, pupid, MYID, &MPI_COMM_PUP);
MPI_Comm_rank(MPI_COMM_PUP, &myid_pup);
/* end of MPI split for processors with ntr>0 */
SHOTINC=1;
RECINC=1;
......@@ -1376,7 +1390,29 @@ 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 (VERBOSE==1 && 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 for forward modelling of STF*/
if (WAVESEP==1) {
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup,no_myid);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished"); /* save pup*/
if (VERBOSE==1 && MYID==0) {
fprintf(FP,"Write PUP data to file...\n ");
sprintf(pf,"%s_pup.su.stf.shot%i",SEIS_FILE,ishot);
outseis_glob(FP,fopen(pf,"w"), 0, fulldata_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
}
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
h=1;
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectionp[h][j]=fulldata_p[recpos_loc[3][i]][j];
}
h++;
}
}
if (VERBOSE==0 && 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);
break;
} /* end of switch (SEISMO) */
......@@ -1907,7 +1943,36 @@ 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 && VERBOSE==1) 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 start wavefield separation...");
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup,no_myid);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished"); /* save pup*/
if (VERBOSE==1 && MYID==0) {
fprintf(FP,"Write PUP data to file...\n ");
sprintf(pf,"%s_pup.su.shot%i",SEIS_FILE,ishot);
outseis_glob(FP,fopen(pf,"w"), 0, fulldata_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
}
if ((VERBOSE==1) && (MYID == 0)) fprintf(FP,"wavefield separation finished.");
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
if (FORWARD_ONLY==0) {
if(ntr>0) {
h=1;
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectionp[h][j]=fulldata_p[recpos_loc[3][i]][j];
}
h++;
}
}
}
}
if (MYID==0 && VERBOSE==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);
break;
} /* end of switch (SEISMO) */
......@@ -3505,6 +3570,29 @@ int main(int argc, char **argv){
}
}
/* wavefield separation */
if (WAVESEP==1) {
catseis(sectionvy, fulldata_vy, recswitch, ntr_glob, MPI_COMM_WORLD);
catseis(sectionp, fulldata_p, recswitch, ntr_glob, MPI_COMM_WORLD);
if (MYID == 0) fprintf(FP,"\n start wavefield separation...");
if (pupid) pup(fulldata_p, fulldata_vy, ntr_glob,recpos,ishot,ns,iter,MPI_COMM_PUP,myid_pup,no_myid);
if ((VERBOSE==1)&&(MYID == 0)) fprintf(FP,"\n wavefield separation finished"); /* save pup*/
if (VERBOSE==1 && MYID==0) {
fprintf(FP,"Write PUP data to file...\n ");
sprintf(pf,"%s_pup.su.step.shot%i",SEIS_FILE,ishot);
outseis_glob(FP,fopen(pf,"w"), 0, fulldata_p,recpos,recpos_loc,ntr_glob,srcpos,1,ns,SEIS_FORMAT,ishot,0);
}
if ((VERBOSE==1) && (MYID == 0)) fprintf(FP,"\n wavefield separation finished");
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast( &fulldata_p[1][1], ntr_glob*ns, MPI_FLOAT, 0, MPI_COMM_WORLD);
h=1;
for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){
sectionp[h][j]=fulldata_p[recpos_loc[3][i]][j];
}
h++;
}
}
/*-------------------------------------------------------------------------------*/
/*------------------ end loop over timesteps (forward model) step length -------*/
......
......@@ -181,7 +181,9 @@ IFOS2D= \
readmod_viscac.c \
time_window_glob.c \
create_trkill_table.c \
filter_frequencies.c
filter_frequencies.c \
pup.c \
fft.c
# -------------
# Targes
......
......@@ -119,6 +119,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 */
......@@ -220,7 +223,13 @@ void exchange_par(void){
fdum[68]=TRKILL_OFFSET_LOWER;
fdum[69]=TRKILL_OFFSET_UPPER;
fdum[70]=LBFGS_SCALE_GRADIENTS;
fdum[74]=LBFGS_SCALE_GRADIENTS;
fdum[70]=VEL;
fdum[71]=DENS;
fdum[72]=ANGTAPI;
fdum[73]=ANGTAPO;
/***********/
/* Integer */
......@@ -375,9 +384,10 @@ void exchange_par(void){
idum[114]=TRKILL_OFFSET;
idum[115]=TRKILL_STF_OFFSET;
idum[116]=TRKILL_STF_OFFSET_INVERT;
idum[117]=WAVESEP;
idum[117]=JOINT_EQUAL_WEIGHTING;
idum[118]=STF_FULL;
idum[118]=JOINT_EQUAL_WEIGHTING;
idum[119]=STF_FULL;
} /** if (MYID == 0) **/
MPI_Barrier(MPI_COMM_WORLD);
......@@ -503,8 +513,14 @@ 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];
LBFGS_SCALE_GRADIENTS=fdum[70];
LBFGS_SCALE_GRADIENTS=fdum[74];
/***********/
/* Integer */
......@@ -662,9 +678,10 @@ void exchange_par(void){
TRKILL_OFFSET=idum[114];
TRKILL_STF_OFFSET=idum[115];
TRKILL_STF_OFFSET_INVERT=idum[116];
WAVESEP=idum[117];
JOINT_EQUAL_WEIGHTING=idum[117];
STF_FULL=idum[118];
JOINT_EQUAL_WEIGHTING=idum[118];
STF_FULL=idum[119];
if ( MYID!=0 && L>0 ) {
FL=vector(1,L);
}
......
......@@ -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,9 @@
#define REQUEST_COUNT 4
#define WORKFLOW_MAX_VAR 13
/****************************/
/* declaration of functions */
/****************************/
void window_cos(float **win, int npad, int nsrc, float it1, float it2, float it3, float it4);
......@@ -131,7 +133,9 @@ void FFT_filt(float ** data, float freqshift, int ntr, int ns,int method);
/*short FFT(short int dir,long m,float *x,float *y);*/
void FFT(int isign, unsigned long nlog, float *re, float *im); /* NR version*/
// void FFT(int isign, unsigned long nlog, float *re, float *im); /* NR version*/
void fft(float *array, float *arrayim, int npad, int dir); /* 1d fft*/
float *filter_frequencies(int *nfrq);
......@@ -209,6 +213,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_vy, int ntr_glob, int** recpos, int ishot, int ns, int iter, MPI_Comm newcomm_nodentr, int myid_pup, int no_myid );
float *rd_sour(int *nts,FILE* fp_source);
......
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* 1D-FFT
* ----------------------------------------------------------------------*/
#include "fd.h"
#include <fftw3.h>
void fft(float *array, float *arrayim, int npad, int dir) {
/* declaration of local variables */
int i,j;
/* declaration of variables for FFTW3*/
fftw_complex *in,*out;
fftw_plan p;
in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * npad);
if (dir==1) {
/* write data into complex vector */
for (j=0;j<npad;j++){
in[j][0]=array[j+1];
in[j][1]=arrayim[j+1];
}
p = fftw_plan_dft_1d(npad, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
fftw_execute(p);
/* scramble data (reorder seismogram to put frequency 0 in the middle) */
for (i=0,j=npad/2;j<npad;j++,i++){
array[j+1]=out[i][0];
array[i+1]=out[j][0];
arrayim[j+1]=out[i][1];
arrayim[i+1]=out[j][1];
}
} else {
/* unscramble data */
for (i=0,j=npad/2;j<npad;j++,i++){
in[i][0]=array[j+1];
in[j][0]=array[i+1];
in[i][1]=arrayim[j+1];
in[j][1]=arrayim[i+1];
}
p = fftw_plan_dft_1d(npad, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
fftw_execute(p);
/* write output into data matrix */
for(j=0;j<npad;j++){
array[j+1] = out[j][0];
arrayim[j+1] = out[j][1];
}
}
fftw_free(in);
fftw_free(out);
fftw_destroy_plan(p);
}
/*-----------------------------------------------------------------------------------------
* 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) {
COMPLEX *C;
int j, i, k;
C = cplxvector(1, NYG*NXG);
if (dir==1) {
k = 1;
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
C[k].re = array[j][i];
C[k].im = 0.0;
k++;
}
forward_fft2f(C, NYG, NXG);
k = 1;
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
array[j][i] = C[k].re;
arrayim[j][i] = C[k].im;
k++;
}
} else {
k = 1;
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
C[k].re = array[j][i];
C[k].im = arrayim[j][i];
k++;
}
inverse_fft2f(C, NYG, NXG);
k = 1;
for (i=1; i<=NXG; i++)
for (j=1; j<=NYG; j++) {
array[j][i] = C[k].re;
arrayim[j][i] = 0.0;
k++;
}
}
free_cplxvector(C, 1, NYG*NXG);
}
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* 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 from Standard C to ANSI C by Stefan Gustavson
*
* Rows and columns must each be an integral power of 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,rows,cols) if (1){\
register int j,k;\
for (j = row, k = 0 ; k < cols ; j+=rows, k++){\
stageBuff[k].re = buffer[j].re;\
stageBuff[k].im = buffer[j].im;\
}\
} else {}
#define StoreRow(buffer,row,rows,cols) if (1){\
register int j,k;\
for (j = row, k = 0 ; k < cols ; j+=rows, 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 = col*rows,k = 0 ; k < rows ; j++, k++){\
stageBuff1[k].re = buffer[j].re;\
stageBuff1[k].im = buffer[j].im;\
}\
} else {}
#define StoreCol(buffer,col,rows,cols) if (1){\
register int j,k;\
for (j = col*rows,k = 0 ; k < rows ; j++, k++){\
buffer[j].re = stageBuff1[k].re;\
buffer[j].im = stageBuff1[k].im;\
}\
} else {}
/* Do something useful with an error message */
#define handle_error(msg) fprintf(FP,msg)
DCOMPLEX *stageBuff; /* buffer to hold a row or column at a time */
DCOMPLEX *stageBuff1; /* 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 */
extern FILE *FP;
/* 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); }
}
/* Allocate space for stageBuff */
int allocateBuffer1(int size) {
stageBuff1 = (DCOMPLEX *) malloc(size*sizeof(DCOMPLEX));
if (stageBuff1==(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);
free((char *)stageBuff1);
}
/* 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;
double *cr0, *cr1, *cr2, *cr3, *cr4, *cr5, *cr6, *cr7;
double *ci0, *ci1, *ci2, *ci3, *ci4, *ci5, *ci6, *ci7;
cr0 = &(cc0[0].re);
cr1 = &(cc1[0].re);
cr2 = &(cc2[0].re);
cr3 = &(cc3[0].re);
cr4 = &(cc4[0].re);
cr5 = &(cc5[0].re);
cr6 = &(cc6[0].re);
cr7 = &(cc7[0].re);
ci0 = &(cc0[0].im);
ci1 = &(cc1[0].im);
ci2 = &(cc2[0].im);
ci3 = &(cc3[0].im);
ci4 = &(cc4[0].im);
ci5 = &(cc5[0].im);
ci6 = &(cc6[0].im);
ci7 = &(cc7[0].im);
scale = TWOPI/length;
for (j=1; j<=nxtlt; j++) {
arg = (j-1)*scale;
c1 = cos(arg);
s1 = sin(arg);
c2 = c1*c1 - s1*s1;
s2 = c1*s1 + c1*s1;
c3 = c1*c2 - s1*s2;
s3 = c2*s1 + s2*c1;
c4 = c2*c2 - s2*s2;
s4 = c2*s2 + c2*s2;
c5 = c2*c3 - s2*s3;
s5 = c3*s2 + s3*c2;
c6 = c3*c3 - s3*s3;
s6 = c3*s3 + c3*s3;
c7 = c3*c4 - s3*s4;
s7 = c4*s3 + s4*c3;
for (k=j; k<=nthpo; k+=length) {
kk = (k-1)*2; /* index by twos; re & im alternate */
ar0 = cr0[kk] + cr4[kk];
ar1 = cr1[kk] + cr5[kk];
ar2 = cr2[kk] + cr6[kk];
ar3 = cr3[kk] + cr7[kk];
ar4 = cr0[kk] - cr4[kk];
ar5 = cr1[kk] - cr5[kk];
ar6 = cr2[kk] - cr6[kk];
ar7 = cr3[kk] - cr7[kk];
ai0 = ci0[kk] + ci4[kk];
ai1 = ci1[kk] + ci5[kk];
ai2 = ci2[kk] + ci6[kk];
ai3 = ci3[kk] + ci7[kk];
ai4 = ci0[kk] - ci4[kk];
ai5 = ci1[kk] - ci5[kk];
ai6 = ci2[kk] - ci6[kk];
ai7 = ci3[kk] - ci7[kk];
br0 = ar0 + ar2;
br1 = ar1 + ar3;
br2 = ar0 - ar2;
br3 = ar1 - ar3;
br4 = ar4 - ai6;
br5 = ar5 - ai7;
br6 = ar4 + ai6;
br7 = ar5 + ai7;
bi0 = ai0 + ai2;
bi1 = ai1 + ai3;
bi2 = ai0 - ai2;
bi3 = ai1 - ai3;
bi4 = ai4 + ar6;
bi5 = ai5 + ar7;
bi6 = ai4 - ar6;
bi7 = ai5 - ar7;
cr0[kk] = br0 + br1;
ci0[kk] = bi0 + bi1;
if (j>1) {
cr1[kk] = c4*(br0-br1) - s4*(bi0-bi1);
cr2[kk] = c2*(br2-bi3) - s2*(bi2+br3);
cr3[kk] = c6*(br2+bi3) - s6*(bi2-br3);
ci1[kk] = c4*(bi0-bi1) + s4*(br0-br1);
ci2[kk] = c2*(bi2+br3) + s2*(br2-bi3);
ci3[kk] = c6*(bi2-br3) + s6*(br2+bi3);
tr = IRT2*(br5-bi5);
ti = IRT2*(br5+bi5);
cr4[kk] = c1*(br4+tr) - s1*(bi4+ti);
ci4[kk] = c1*(bi4+ti) + s1*(br4+tr);
cr5[kk] = c5*(br4-tr) - s5*(bi4-ti);
ci5[kk] = c5*(bi4-ti) + s5*(br4-tr);
tr = -IRT2*(br7+bi7);
ti = IRT2*(br7-bi7);
cr6[kk] = c3*(br6+tr) - s3*(bi6+ti);
ci6[kk] = c3*(bi6+ti) + s3*(br6+tr);
cr7[kk] = c7*(br6-tr) - s7*(bi6-ti);
ci7[kk] = c7*(bi6-ti) + s7*(br6-tr);
} else {
cr1[kk] = br0 - br1;
cr2[kk] = br2 - bi3;
cr3[kk] = br2 + bi3;
ci1[kk] = bi0 - bi1;
ci2[kk] = bi2 + br3;
ci3[kk] = bi2 - br3;
tr = IRT2*(br5-bi5);
ti = IRT2*(br5+bi5);
cr4[kk] = br4 + tr;
ci4[kk] = bi4 + ti;
cr5[kk] = br4 - tr;
ci5[kk] = bi4 - ti;
tr = -IRT2*(br7+bi7);
ti = IRT2*(br7-bi7);
cr6[kk] = br6 + tr;
ci6[kk] = bi6 + ti;
cr7[kk] = br6 - tr;
ci7[kk] = bi6 - ti;
}
}
}
}
/*
* FFT842 (Name kept from the original Fortran version)
* This routine replaces the input DCOMPLEX vector by its
* finite discrete complex fourier transform if in==FFT_FORWARD.
* It replaces the input DCOMPLEX vector by its finite discrete
* complex inverse fourier transform if in==FFT_INVERSE.
*
* The implementation is a radix-2 FFT, but with faster shortcuts for
* radix-4 and radix-8. It performs as many radix-8 iterations as
* possible, and then finishes with a radix-2 or -4 iteration if needed.
*/
void FFT842(int direction, int n, DCOMPLEX *b)
/* direction: FFT_FORWARD or FFT_INVERSE
* n: length of vector
* *b: input vector */
{
double fn, r, fi;
int L[16],L1,L2,L3,L4,L5,L6,L7,L8,L9,L10,L11,L12,L13,L14,L15;
/* int j0,j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14;*/
int j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14;
int i, j, ij, ji, ij1, ji1;
/* int nn, n2pow, n8pow, nthpo, ipass, nxtlt, length;*/
int n2pow, n8pow, nthpo, ipass, nxtlt, length;
n2pow = fastlog2(n);
nthpo = n;
fn = 1.0 / (double)nthpo; /* Scaling factor for inverse transform */
if (direction==FFT_FORWARD)
/* Conjugate the input */
for (i=0; i<n; i++) {
b[i].im = -b[i].im;
}
if (direction==FFT_INVERSE)
/* Scramble the inputs */
for (i=0,j=n/2; j<n; i++,j++) {
r = b[j].re;
fi = b[j].im;
b[j].re = b[i].re;
b[j].im = b[i].im;
b[i].re = r;
b[i].im = fi;
}
n8pow = n2pow/3;
if (n8pow) {
/* Radix 8 iterations */
for (ipass=1; ipass<=n8pow; ipass++) {
nxtlt = 0x1 << (n2pow - 3*ipass);
length = 8*nxtlt;
R8TX(nxtlt, nthpo, length,
b, b+nxtlt, b+2*nxtlt, b+3*nxtlt,
b+4*nxtlt, b+5*nxtlt, b+6*nxtlt, b+7*nxtlt);
}
}
if (n2pow%3 == 1) {
/* A final radix 2 iteration is needed */
R2TX(nthpo, b, b+1);
}
if (n2pow%3 == 2) {
/* A final radix 4 iteration is needed */
R4TX(nthpo, b, b+1, b+2, b+3);
}
for (j=1; j<=15; j++) {
L[j] = 1;
if (j-n2pow <= 0) { L[j] = 0x1 << (n2pow + 1 - j); }
}
L15=L[1];
L14=L[2];
L13=L[3];
L12=L[4];
L11=L[5];
L10=L[6];
L9=L[7];
L8=L[8];
L7=L[9];
L6=L[10];
L5=L[11];
L4=L[12];
L3=L[13];
L2=L[14];
L1=L[15];
ij = 1;
for (j1=1; j1<=L1; j1++)
for (j2=j1; j2<=L2; j2+=L1)
for (j3=j2; j3<=L3; j3+=L2)
for (j4=j3; j4<=L4; j4+=L3)
for (j5=j4; j5<=L5; j5+=L4)
for (j6=j5; j6<=L6; j6+=L5)
for (j7=j6; j7<=L7; j7+=L6)
for (j8=j7; j8<=L8; j8+=L7)
for (j9=j8; j9<=L9; j9+=L8)
for (j10=j9; j10<=L10; j10+=L9)
for (j11=j10; j11<=L11; j11+=L10)
for (j12=j11; j12<=L12; j12+=L11)
for (j13=j12; j13<=L13; j13+=L12)
for (j14=j13; j14<=L14; j14+=L13)
for (ji=j14; ji<=L15; ji+=L14) {
ij1 = ij-1;
ji1 = ji-1;
if (ij-ji<0) {
r = b[ij1].re;
b[ij1].re = b[ji1].re;
b[ji1].re = r;
fi = b[ij1].im;
b[ij1].im = b[ji1].im;
b[ji1].im = fi;
}
ij++;
}
if (direction==FFT_FORWARD) /* Take conjugates & unscramble outputs */
for (i=0,j=n/2; j<n; i++,j++) {
r = b[j].re;
fi = b[j].im;
b[j].re = b[i].re;
b[j].im = -b[i].im;
b[i].re = r;
b[i].im = -fi;
}
if (direction==FFT_INVERSE) /* Scale outputs */
for (i=0; i<nthpo; i++) {
b[i].re *= fn;
b[i].im *= fn;
}
}
/*
* Compute 2D DFT on float data:
* forward if direction==FFT_FORWARD,
* inverse if direction==FFT_INVERSE.
*/
int fft2f(COMPLEX *array, int rows, int cols, int direction) {
int i;
if (!power_of_2(rows) || !power_of_2(cols)) {
handle_error("\n ****************************************************");
handle_error("\n fft: input array must have dimensions a power of 2! \n");
handle_error(" ****************************************************\n");
return (ERROR);
}
/* Allocate 1D buffer */
bigBuff = array;
// maxsize = rows>cols ? rows : cols;
allocateBuffer(cols);
allocateBuffer1(rows);
// if(errflag != NO_ERROR) return(errflag);
/* Compute transform row by row */
if (cols>1)
for (i=0; i<rows; i++) {
LoadRow(bigBuff,i,rows,cols);
FFT842(direction,cols,stageBuff);
StoreRow(bigBuff,i,rows,cols);
}
/* Compute transform column by column */
if (rows>1)
for (i=0; i<cols; i++) {
LoadCol(bigBuff,i,rows,cols);
FFT842(direction,rows,stageBuff1);
StoreCol(bigBuff,i,rows,cols);
}
freeBuffer();
return (NO_ERROR);
}
// /* Finally, the entry points we announce in kube_fft.h */
/* Perform forward 2D transform on a COMPLEX array. */
int forward_fft2f(COMPLEX *array, int rows, int cols) {
return (fft2f(array, rows, cols, FFT_FORWARD));
}
/* Perform inverse 2D transform on a COMPLEX array. */
int inverse_fft2f(COMPLEX *array, int rows, int cols) {
return (fft2f(array, rows, cols, FFT_INVERSE));
}
......@@ -148,3 +148,6 @@ int JOINT_EQUAL_WEIGHTING;
float JOINT_INVERSION_PSV_SH_ALPHA_VS;
float JOINT_INVERSION_PSV_SH_ALPHA_RHO;
int SNAPSHOT_START,SNAPSHOT_END,SNAPSHOT_INCR;
int WAVESEP;
float VEL, DENS, ANGTAPI, ANGTAPO;
\ No newline at end of file
......@@ -17,7 +17,7 @@
-----------------------------------------------------------------------------------------*/
/*------------------------------------------------------------------------
* Write seismograms to disk
* Read seismograms from disk
* ----------------------------------------------------------------------*/
#include "fd.h"
#include "segy.h"
......@@ -27,6 +27,7 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
/* declaration of extern variables */
extern int NDT, MYID;
extern char DATA_DIR[STRING_SIZE];
extern char SEIS_FILE[STRING_SIZE];
extern float TIME, DH, DT, REFREC[4];
char data[STRING_SIZE];
FILE *fpdata;
......@@ -39,7 +40,19 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
sprintf(data,"%s_vy.su.shot%d",DATA_DIR,comp);
}
/* sws 3 -- 6 not used */
if(sws==3){ /* open forward modeled seismic data pup*/
sprintf(data,"%s_pup.su.shot%d",SEIS_FILE,comp);
}
if(sws==4){ /* open forward modeled seismic data pup used for step length calculation */
sprintf(data,"%s_pup.su.step.shot%d",SEIS_FILE,comp);
}
if(sws==5){ /* open seismic data pup */
sprintf(data,"%s_pup.su.shot%d",DATA_DIR,comp);
}
/* SWS==6 not used */
if(sws==7){ /* open convolved seismic data vx */
sprintf(data,"%s_vx.su.conv.shot%d",DATA_DIR,comp);
......@@ -57,7 +70,6 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
sprintf(data,"%s_vz.su.shot%d",DATA_DIR,comp);
}
fpdata = fopen(data,"r");
if (fpdata==NULL) {
if(MYID==0) printf(" Was not able to read %s",data);
......
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
/*
*-------------------------------------------------------------
*
* perform wavefield separation in FK domain and save it to file
*
*
*-------------------------------------------------------------
*/
#include "fd.h"
void pup(float **data_p, float **data_vy, int ntr_glob, int **recpos, int ishot, int ns_glob, int iter, MPI_Comm newcomm_nodentr,int myid_pup,int no_myid) {
extern float DT,DH, VEL, DENS, ANGTAPI, ANGTAPO;
extern char SEIS_FILE[STRING_SIZE];
extern int VERBOSE;
extern FILE *FP;
int i, j, k, nx_pot, ny_pot,pvz;
double t0,t1,t2,t3,t4,t5,t6,t7,t8;
int ntr_start,ntr_end,ns_start,ns_end, ntr,ns;
float **filtwk,*w, *kx;
float nyq_k, nyq_f, dk, df, angi, ango, ka, dtr;
FILE *file;
float *col,*col_im,*row,*row_im;
float **data_pim,**data_im1,**data_vyim,**data_pup,**data_pupim,**data_pt,**data_t1,**data_vyt,**data_pup1;
char pf[STRING_SIZE], pf1[STRING_SIZE], pf2[STRING_SIZE];
t1=MPI_Wtime();
/* double the dimensions to avoid artefacts and recompute the dimensions to become values of 2^m */
nx_pot = (int)(pow(2,(ceil(log(ntr_glob*2)/log(2)))));
ny_pot = (int)(pow(2,(ceil(log(ns_glob*2)/log(2)))));
fprintf(FP,"\n ny_pot: %i; nx_pot: %i \n",ny_pot,nx_pot);
ntr_start=nx_pot/(no_myid)*(myid_pup)+1;
ntr_end=nx_pot/(no_myid)*(myid_pup+1);
ns_start=ny_pot/(no_myid)*(myid_pup)+1;
ns_end=ny_pot/(no_myid)*(myid_pup+1);
ns=ns_end-ns_start+1;
ntr=ntr_end-ntr_start+1;
printf("MYID: %i; ntr_start: %i; ntr_end: %i; ns_start: %i; ns_end: %i; ntr: %i, ns: %i \n",myid_pup,ntr_start,ntr_end,ns_start,ns_end,ntr,ns);
data_im1=matrix(1,ny_pot,1,nx_pot);
data_t1=matrix(1,ny_pot,1,nx_pot);
data_pup1=matrix(1,ny_pot,1,nx_pot);
data_pt=matrix(1,ns,1,nx_pot);
data_pim=matrix(1,ns,1,nx_pot);
data_vyt=matrix(1,ns,1,nx_pot);
data_vyim=matrix(1,ns,1,nx_pot);
data_pup=matrix(1,ns,1,nx_pot);
data_pupim=matrix(1,ns,1,nx_pot);
col = vector(1,ny_pot);
col_im = vector(1,ny_pot);
row = vector(1,nx_pot);
row_im = vector(1,nx_pot);
dtr=fabs((recpos[1][ntr_glob]-recpos[1][1])/(ntr_glob-1)*DH);
if (VERBOSE==1) {
fprintf(FP,"\n Receiver sampling was calculatet to dtr: %f m \n ",dtr);
}
filtwk=matrix(1,ns,1,nx_pot);
w=vector(1,ny_pot);
kx=vector(1,nx_pot/2);
nyq_k = 1/(2*dtr); /* Nyquist of data in first dimension */
nyq_f = 1/(2*DT); /* Nyquist of data in time dimension */
dk = 1/(nx_pot*dtr); /* Wavenumber increment */
df = 1/(ny_pot*DT); /* Frequency increment */
for (j=1; j<=ny_pot; j++) {
w[j]=(-nyq_f+(j-1)*df)*2*PI; /* vector of angular frequencies*/
}
for (j=1; j<=nx_pot/2; j++) {
kx[j]=((j-1)*dk)*2*PI; /* vector of angular wavenumbers*/
}
angi = ANGTAPI * PI/180.0; /* angle at which taper begin */
ango = ANGTAPO * PI/180.0; /* angle at which taper ends */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* start calculating filter */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* calculate Filter vor wavefield separation for first quadrant (after Kluever, 2008) */
if (VERBOSE==1 && myid_pup==0) {
fprintf(FP,"Calculate Filter for wavefield separation... \n ");
}
for (k=ns_start,j=1; k<=ns_end; k++,j++) {
ka = fabs(w[k])/VEL;
for (i=1; i<=nx_pot/2; i++) {
if (kx[i] < ka*sin(angi)) {
filtwk[j][i+nx_pot/2] = (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
filtwk[j][nx_pot/2-(i-1)] = filtwk[j][i+nx_pot/2];
} else {
if (kx[i] < ka*sin(ango) && kx[i] > ka*sin(angi)) {
filtwk[j][i+nx_pot/2] = 0.5*(1.0+cos(PI/((ka*sin(ango))-(ka*sin(angi))) * (kx[i]-(ka*sin(angi))))) * (DENS*fabs(w[k]))/sqrt((ka*ka)-(kx[i]*kx[i]));
filtwk[j][nx_pot/2-(i-1)] = filtwk[j][i+nx_pot/2];
} else {
filtwk[j][i+nx_pot/2] = 0.0;
filtwk[j][nx_pot/2-(i-1)] = 0.0;
}
}
}
}
t0=MPI_Wtime();
fprintf(FP,"\ntime for calculating filter: %f\n",t0-t1);
/* write filter to file*/
if (VERBOSE==1 && iter==1 && ishot==1) {
sprintf(pf1,"%s_spectrum.filter_shot%i.it%i.%i.%i.%i.bin",SEIS_FILE,ishot, iter, ns, nx_pot,myid_pup);
file = fopen(pf1,"wb");
for (i=1; i<=nx_pot; i++)
for (j=1; j<=ns; j++) {
fwrite(&filtwk[j][i],sizeof(float),1,file);
}
fclose(file);
}
t4=MPI_Wtime();
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* start fft for pressure and vz */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* columns */
for (pvz=1; pvz<=2; pvz++) {
if (VERBOSE==1 && pvz==1) fprintf(FP,"\n starting fft of pressure wavefield...\n");
else if (VERBOSE==1 && pvz==2) fprintf(FP,"\n starting fft of vertical velocity wavefield...\n");
for (k=ntr_start; k<=ntr_end; k++) {
if (k<=ntr_glob) {
for (i=1; i<=ns_glob; i++) {
if (pvz==1){
col[i]=data_p[k][i];
} else if (pvz==2)
col[i]=data_vy[k][i];
}
fft(col, col_im, ny_pot,1);
for (i=1; i<=ny_pot; i++) {
data_t1[i][k]=col[i];
data_im1[i][k]=col_im[i];
col[i]=0;
col_im[i]=0;
}
}
}
MPI_Allreduce(MPI_IN_PLACE, &data_t1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
MPI_Allreduce(MPI_IN_PLACE, &data_im1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
/* rows */
for (i=ns_start,j=1; i<=ns_end; i++,j++) {
for (k=1; k<=nx_pot; k++) {
row[k]=data_t1[i][k];
row_im[k]=data_im1[i][k];
}
fft(row, row_im, nx_pot,1);
for (k=1; k<=nx_pot; k++) {
if (pvz==1) {
data_pt[j][k]=row[k];
data_pim[j][k]=row_im[k];
} else if (pvz==2) {
data_vyt[j][k]=row[k];
data_vyim[j][k]=row_im[k];
}
}
}
for (i=1; i<=ny_pot; i++) {
for (k=1; k<=nx_pot; k++) {
data_t1[i][k]=0;
data_im1[i][k]=0;
}
}
// if (VERBOSE==1 && iter==1 && ishot==1 && pvz==1) {
// sprintf(pf1,"%s_spectrum.p_shot%i.it%i.%i.%i.%i.bin",SEIS_FILE,ishot, iter, ns, nx_pot,myid_pup);
// file = fopen(pf1,"wb");
// for (i=1; i<=nx_pot; i++)
// for (j=1; j<=ns; j++) {
// fwrite(&data_pt[j][i],sizeof(float),1,file);
// }
// fclose(file);
// } else if (VERBOSE==1 && iter==1 && ishot==1 && pvz==2) {
// sprintf(pf1,"%s_spectrum.vz_shot%i.it%i.%i.%i.%i.bin",SEIS_FILE,ishot, iter, ns, nx_pot,myid_pup);
// file = fopen(pf1,"wb");
// for (i=1; i<=nx_pot; i++)
// for (j=1; j<=ns; j++) {
// fwrite(&data_vyt[j][i],sizeof(float),1,file);
// }
// fclose(file);
// }
}
for (k=1; k<=ntr_glob; k++) {
for (i=1; i<=ns_glob; i++) {
data_p[k][i]=0;
}
}
t5=MPI_Wtime();
fprintf(FP,"\ntime for parallised ffts: %f\n",t5-t4);
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* perform wavefield separation */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
if (VERBOSE==1) fprintf(FP,"\n starting wavefield separation...\n");
for (i=1; i<=ns; i++) {
for (k=1; k<=nx_pot; k++) {
data_pup[i][k]=0.5*(data_pt[i][k]-(filtwk[i][k]*data_vyt[i][k]));
data_pupim[i][k]=0.5*(data_pim[i][k]-(filtwk[i][k]*data_vyim[i][k]));
}
}
// if (VERBOSE==1 && iter==1 && ishot==1) {
// sprintf(pf1,"%s_spectrum.pup_shot%i.it%i.%i.%i.%i.bin",SEIS_FILE,ishot, iter, ns, nx_pot,myid_pup);
// file = fopen(pf1,"wb");
// for (i=1; i<=nx_pot; i++)
// for (j=1; j<=ns; j++) {
// fwrite(&data_pup[j][i],sizeof(float),1,file);
// }
// fclose(file);
// }
t6=MPI_Wtime();
fprintf(FP,"\ntime for wavefield separation: %f\n",t6-t5);
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
/* start inverse fft for pup */
/*%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
if (VERBOSE==1) fprintf(FP,"\n starting inverse fft of pup wavefield...\n");
/* rows */
for (i=ns_start,j=1; i<=ns_end; i++,j++) {
for (k=1; k<=nx_pot; k++) {
row[k]=data_pup[j][k];
row_im[k]=data_pupim[j][k];
}
fft(row, row_im, nx_pot,-1);
for (k=1; k<=nx_pot; k++) {
data_t1[i][k]=row[k];
data_im1[i][k]=row_im[k];
}
}
MPI_Allreduce(MPI_IN_PLACE, &data_t1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
MPI_Allreduce(MPI_IN_PLACE, &data_im1[1][1], nx_pot*ny_pot, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
for (k=ntr_start; k<=ntr_end; k++) {
for (i=1; i<=ny_pot; i++) {
col[i]=data_t1[i][k];
col_im[i]=data_im1[i][k];
}
fft(col, col_im, ny_pot,-1);
if (k<=ntr_glob){
for (i=1; i<=ns_glob; i++) {
data_p[k][i]=col[i];
}
}
}
MPI_Allreduce(MPI_IN_PLACE, &data_p[1][1], ntr_glob*ns_glob, MPI_FLOAT, MPI_SUM, newcomm_nodentr);
if (VERBOSE==1) fprintf(FP,"\n wavefield separation finished.\n");
t8=MPI_Wtime();
fprintf(FP,"\ntime for reverse fft: %f\n",t8-t6);
fprintf(FP,"\ntime for pup: %f\n",t8-t1);
free_matrix(filtwk,1,ns,1,nx_pot);
free_vector(w,1,ny_pot);
free_vector(kx,1,ny_pot);
free_matrix(data_im1, 1,ny_pot,1,nx_pot);
free_matrix(data_t1, 1,ny_pot,1,nx_pot);
free_matrix(data_pup1, 1,ny_pot,1,nx_pot);
free_matrix(data_pt, 1,ns,1,nx_pot);
free_matrix(data_pim, 1,ns,1,nx_pot);
free_matrix(data_vyt, 1,ns,1,nx_pot);
free_matrix(data_vyim, 1,ns,1,nx_pot);
free_matrix(data_pup, 1,ns,1,nx_pot);
free_matrix(data_pupim, 1,ns,1,nx_pot);
free_vector(col,1,ny_pot);
free_vector(col_im,1,ny_pot);
free_vector(row,1,nx_pot);
free_vector(row_im,1,nx_pot);
}
\ No newline at end of file
......@@ -127,6 +127,9 @@ void read_par_json(FILE *fp, char *fileinp){
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;
extern int STF_FULL;
/* definition of local variables */
......@@ -341,6 +344,8 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("IDY",number_readobjects,&IDY,varname_list, value_list)){
IDY=1;
fprintf(fp,"Variable IDY is set to default value %d.\n",IDY);}
/*=================================
section seismogramm parameters
......@@ -397,6 +402,7 @@ void read_par_json(FILE *fp, char *fileinp){
declare_error("Variable REFRECY could not be retrieved from the json input file!");
}
}
/*=================================
section general model and log parameters
......@@ -978,6 +984,45 @@ void read_par_json(FILE *fp, char *fileinp){
}
/*=================================
section wavefield separation
=================================*/
if (get_int_from_objectlist("WAVESEP",number_readobjects,&WAVESEP,varname_list, value_list)){
WAVESEP=0;
fprintf(fp,"Variable WAVESEP is set to default value %i.\n",WAVESEP);}
else {
if (SEISMO < 4 && WAVESEP == 1) {
WAVESEP=0;
fprintf(fp,"Variable SEISMO has to be set to 4 or 5 to use wavefield separation. WAVESEP is set to %i now.\n",WAVESEP);
}
if (ADJOINT_TYPE != 4 && WAVESEP == 1 && FORWARD_ONLY == 0) {
WAVESEP=0;
fprintf(fp,"Variable ADJOINT_TYPE has to be set to 4 to use wavefield separation. WAVESEP is set to %i now.\n",WAVESEP);
}
if (WAVETYPE == 2 && WAVESEP == 1) {
WAVESEP=0;
fprintf(fp,"Variable WAVETYPE has to be set to 1 or 3 to use wavefield separation. WAVESEP is set to %i now.\n",WAVESEP);
}
if (WAVESEP==1) {
if (get_float_from_objectlist("VEL",number_readobjects,&VEL,varname_list, value_list)){
VEL=1500;
fprintf(fp,"Variable VEL is set to default value %f.\n",VEL);}
if (get_float_from_objectlist("DENS",number_readobjects,&DENS,varname_list, value_list)){
DENS=1000;
fprintf(fp,"Variable DENS is set to default value %f.\n",DENS);}
if (get_float_from_objectlist("ANGTAPI",number_readobjects,&ANGTAPI,varname_list, value_list)){
ANGTAPI=0;
fprintf(fp,"Variable ANGTAPI is set to default value %f.\n",ANGTAPI);}
if (get_float_from_objectlist("ANGTAPO",number_readobjects,&ANGTAPO,varname_list, value_list)){
ANGTAPO=70;
fprintf(fp,"Variable ANGTAPO is set to default value %f.\n",ANGTAPO);}
}
}
fprintf(fp,"\nEnd of setting default values\n");
fprintf(fp,"=====================================\n\n");
......
......@@ -27,7 +27,7 @@ void seismo_ssg(int lsamp, int ntr, int **recpos, float **sectionvx,
float **sectionvy,float **sectionvz, float **sectionp, float **sectioncurl, float **sectiondiv,
float **vx, float **vy, float **vz, float **sxx, float **syy, float **sp, float **pi, float **u, float *hc){
extern int NDT, SEISMO, FDORDER, ACOUSTIC,WAVETYPE;
extern int NDT, SEISMO, FDORDER, ACOUSTIC<