Skip to content
GitLab
Menu
Projects
Groups
Snippets
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Sign in
Toggle navigation
Menu
Open sidebar
GPIAG-Software
IFOS2D
Commits
0803b755
Commit
0803b755
authored
Apr 05, 2016
by
niklas.thiel
Browse files
included VERBOSE, cleaned up
parent
bcca82b5
Changes
4
Show whitespace changes
Inline
Side-by-side
src/IFOS2D.c
View file @
0803b755
...
@@ -1428,7 +1428,7 @@ int main(int argc, char **argv){
...
@@ -1428,7 +1428,7 @@ int main(int argc, char **argv){
if
(
WAVESEP
==
1
)
{
if
(
WAVESEP
==
1
)
{
if
(
MYID
==
0
)
fprintf
(
FP
,
"
\n
start wavefield separation..."
);
if
(
MYID
==
0
)
fprintf
(
FP
,
"
\n
start wavefield separation..."
);
if
(
MYID
==
0
)
pup
(
fulldata_p
,
fulldata_vy
,
FP
,
ntr_glob
,
recpos
,
recpos_loc
,
ntr_glob
,
srcpos
,
ishot
,
ns
,
iter
,
1
);
if
(
MYID
==
0
)
pup
(
fulldata_p
,
fulldata_vy
,
FP
,
ntr_glob
,
recpos
,
recpos_loc
,
ntr_glob
,
srcpos
,
ishot
,
ns
,
iter
,
1
);
if
(
MYID
==
0
)
fprintf
(
FP
,
"
\n
wavefield separation finished"
);
if
((
VERBOSE
)
&&
(
MYID
==
0
)
)
fprintf
(
FP
,
"
\n
wavefield separation finished"
);
MPI_Barrier
(
MPI_COMM_WORLD
);
MPI_Barrier
(
MPI_COMM_WORLD
);
inseis
(
fprec
,
ishot
,
sectionread
,
ntr_glob
,
ns
,
3
,
iter
);
inseis
(
fprec
,
ishot
,
sectionread
,
ntr_glob
,
ns
,
3
,
iter
);
h
=
1
;
h
=
1
;
...
@@ -1463,7 +1463,7 @@ int main(int argc, char **argv){
...
@@ -1463,7 +1463,7 @@ int main(int argc, char **argv){
if
(
ADJOINT_TYPE
==
4
){
if
(
ADJOINT_TYPE
==
4
){
if
(
WAVESEP
==
1
)
{
inseis
(
fprec
,
ishot
,
sectionp_obs
,
ntr_glob
,
ns
,
5
,
iter
);
if
(
WAVESEP
==
1
)
{
inseis
(
fprec
,
ishot
,
sectionp_obs
,
ntr_glob
,
ns
,
5
,
iter
);
}
else
{
inseis
(
fprec
,
ishot
,
sectionp_obs
,
ntr_glob
,
ns
,
9
,
iter
);}
}
else
{
inseis
(
fprec
,
ishot
,
sectionp_obs
,
ntr_glob
,
ns
,
9
,
iter
);}
timedomain_filt
(
sectionp_obs
,
F
C
,
ORDER
,
ntr_glob
,
ns
,
1
);
timedomain_filt
(
sectionp_obs
,
F
_LOW_PASS
,
ORDER
,
ntr_glob
,
ns
,
1
);
}
}
}
}
...
@@ -1974,7 +1974,7 @@ int main(int argc, char **argv){
...
@@ -1974,7 +1974,7 @@ int main(int argc, char **argv){
if
(
WAVESEP
==
1
)
{
if
(
WAVESEP
==
1
)
{
if
(
MYID
==
0
)
fprintf
(
FP
,
"
\n
start wavefield separation..."
);
if
(
MYID
==
0
)
fprintf
(
FP
,
"
\n
start wavefield separation..."
);
if
(
MYID
==
0
)
pup
(
fulldata_p
,
fulldata_vy
,
FP
,
ntr_glob
,
recpos
,
recpos_loc
,
ntr_glob
,
srcpos
,
ishot
,
ns
,
iter
,
1
);
if
(
MYID
==
0
)
pup
(
fulldata_p
,
fulldata_vy
,
FP
,
ntr_glob
,
recpos
,
recpos_loc
,
ntr_glob
,
srcpos
,
ishot
,
ns
,
iter
,
1
);
if
(
MYID
==
0
)
fprintf
(
FP
,
"wavefield separation finished."
);
if
((
VERBOSE
)
&&
(
MYID
==
0
)
)
fprintf
(
FP
,
"wavefield separation finished."
);
if
(
FORWARD_ONLY
==
0
)
{
if
(
FORWARD_ONLY
==
0
)
{
MPI_Barrier
(
MPI_COMM_WORLD
);
MPI_Barrier
(
MPI_COMM_WORLD
);
inseis
(
fprec
,
ishot
,
sectionread
,
ntr_glob
,
ns
,
3
,
iter
);
inseis
(
fprec
,
ishot
,
sectionread
,
ntr_glob
,
ns
,
3
,
iter
);
...
@@ -3515,7 +3515,7 @@ int main(int argc, char **argv){
...
@@ -3515,7 +3515,7 @@ int main(int argc, char **argv){
catseis
(
sectionp
,
fulldata_p
,
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
(
MYID
==
0
)
fprintf
(
FP
,
"
\n
start wavefield separation..."
);
if
(
MYID
==
0
)
pup
(
fulldata_p
,
fulldata_vy
,
FP
,
ntr_glob
,
recpos
,
recpos_loc
,
ntr_glob
,
srcpos
,
ishot
,
ns
,
iter
,
2
);
if
(
MYID
==
0
)
pup
(
fulldata_p
,
fulldata_vy
,
FP
,
ntr_glob
,
recpos
,
recpos_loc
,
ntr_glob
,
srcpos
,
ishot
,
ns
,
iter
,
2
);
if
(
MYID
==
0
)
fprintf
(
FP
,
"
\n
wavefield separation finished"
);
if
((
VERBOSE
)
&&
(
MYID
==
0
)
)
fprintf
(
FP
,
"
\n
wavefield separation finished"
);
MPI_Barrier
(
MPI_COMM_WORLD
);
MPI_Barrier
(
MPI_COMM_WORLD
);
inseis
(
fprec
,
ishot
,
sectionread
,
ntr_glob
,
ns
,
4
,
iter
);
inseis
(
fprec
,
ishot
,
sectionread
,
ntr_glob
,
ns
,
4
,
iter
);
h
=
1
;
h
=
1
;
...
...
src/fft2d.c
View file @
0803b755
/*
/*-----------------------------------------------------------------------------------------
------------------------------------------------------------------------------
* Copyright (C) 2016 For the list of authors, see file AUTHORS.
*
* Copyright (C) 2013, see file AUTHORS for list of authors/contributors
*
*
* This file is part of
PROTEU
S.
* This file is part of
IFO
S.
*
*
*
PROTEU
S is free software: you can redistribute it and/or modify
*
IFO
S is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, version 2.0 of the License only.
* the Free Software Foundation, version 2.0 of the License only.
*
*
*
PROTEU
S is distributed in the hope that it will be useful,
*
IFO
S is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
* GNU General Public License for more details.
*
*
* You should have received a copy of the GNU General Public License
* You should have received a copy of the GNU General Public License
* along with PROTEUS. See file COPYING and/or
* along with IFOS. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
* <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
*
------------------------------------------------------------------------------
*/
/*------------------------------------------------------------------------
* Forward and inverse discrete 2D Fourier transforms
* ----------------------------------------------------------------------*/
/**************************************************************
/*
* kube-gustavson-fft.c
*
* Forward and inverse discrete 2D Fourier transforms.
*
* Originally a FORTRAN implementation published in:
* Originally a FORTRAN implementation published in:
* Programming for Digital Signal Processing, IEEE Press 1979,
* Programming for Digital Signal Processing, IEEE Press 1979,
* Section 1, by G. D. Bergland and M. T. Dolan.
* 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
* 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.
* 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"
#include "fd.h"
...
...
src/pup.c
View file @
0803b755
...
@@ -31,7 +31,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -31,7 +31,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
extern
float
DT
,
DH
;
extern
float
DT
,
DH
;
extern
float
VEL
,
DENS
,
ANGTAPI
,
ANGTAPO
;
extern
float
VEL
,
DENS
,
ANGTAPI
,
ANGTAPO
;
extern
char
SEIS_FILE
[
STRING_SIZE
];
extern
char
SEIS_FILE
[
STRING_SIZE
];
extern
int
SEIS_FORMAT
;
extern
int
SEIS_FORMAT
,
VERBOSE
;
extern
FILE
*
FP
;
extern
FILE
*
FP
;
int
i
,
j
,
k
,
nx_pot
,
ny_pot
;
int
i
,
j
,
k
,
nx_pot
,
ny_pot
;
...
@@ -44,7 +44,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -44,7 +44,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
/* calculate receiver sampling */
/* calculate receiver sampling */
dtr
=
(
recpos
[
1
][
ntr_glob
]
-
recpos
[
1
][
1
])
/
(
ntr_glob
-
1
)
*
DH
;
dtr
=
(
recpos
[
1
][
ntr_glob
]
-
recpos
[
1
][
1
])
/
(
ntr_glob
-
1
)
*
DH
;
fprintf
(
FP
,
"
\n
Receiver sampling was calculatet to dtr: %f m
\n
"
,
dtr
);
if
(
VERBOSE
)
fprintf
(
FP
,
"
\n
Receiver sampling was calculatet to dtr: %f m
\n
"
,
dtr
);
/* double the dimensions to avoid artefacts and recompute the dimensions to become values of 2^m */
/* 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
)))));
nx_pot
=
(
int
)(
pow
(
2
,(
ceil
(
log
(
ntr_glob
*
2
)
/
log
(
2
)))));
...
@@ -92,7 +92,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -92,7 +92,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
ango
=
ANGTAPO
*
PI
/
180
.
0
;
/* angle at which taper ends */
ango
=
ANGTAPO
*
PI
/
180
.
0
;
/* angle at which taper ends */
/* calculate Filter vor wavefield separation for first quadrant (after Kluever, 2008) */
/* calculate Filter vor wavefield separation for first quadrant (after Kluever, 2008) */
fprintf
(
FP
,
"Calculate Filter for wavefield separation...
\n
"
);
if
(
VERBOSE
)
fprintf
(
FP
,
"Calculate Filter for wavefield separation...
\n
"
);
for
(
k
=
1
;
k
<=
ny_pot
/
2
;
k
++
)
{
for
(
k
=
1
;
k
<=
ny_pot
/
2
;
k
++
)
{
ka
=
fabs
(
w
[
k
])
/
VEL
;
ka
=
fabs
(
w
[
k
])
/
VEL
;
for
(
i
=
1
;
i
<=
nx_pot
/
2
;
i
++
)
{
for
(
i
=
1
;
i
<=
nx_pot
/
2
;
i
++
)
{
...
@@ -141,7 +141,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -141,7 +141,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
}
}
}
}
fprintf
(
FP
,
"2D FFT of p-data and vz-data...
\n
"
);
if
(
VERBOSE
)
fprintf
(
FP
,
"2D FFT of p-data and vz-data...
\n
"
);
/* transformation into fk domain */
/* transformation into fk domain */
fft2
(
data_pt
,
data_pim
,
ny_pot
,
nx_pot
,
1
);
fft2
(
data_pt
,
data_pim
,
ny_pot
,
nx_pot
,
1
);
...
@@ -166,7 +166,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -166,7 +166,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
fclose
(
file
);
fclose
(
file
);
}
}
fprintf
(
FP
,
"Perform wavefield separation...
\n
"
);
if
(
VERBOSE
)
fprintf
(
FP
,
"Perform wavefield separation...
\n
"
);
/* perform wavefield separation */
/* perform wavefield separation */
for
(
k
=
1
;
k
<=
nx_pot
;
k
++
)
{
for
(
k
=
1
;
k
<=
nx_pot
;
k
++
)
{
...
@@ -184,12 +184,12 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -184,12 +184,12 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
fclose
(
file
);
fclose
(
file
);
}
}
fprintf
(
FP
,
"Reverse 2D FFT...
\n
"
);
if
(
VERBOSE
)
fprintf
(
FP
,
"Reverse 2D FFT...
\n
"
);
/* reverse 2d fft */
/* reverse 2d fft */
fft2
(
data_pup
,
data_pupim
,
ny_pot
,
nx_pot
,
-
1
);
fft2
(
data_pup
,
data_pupim
,
ny_pot
,
nx_pot
,
-
1
);
fprintf
(
FP
,
"Write PUP data to file...
\n
"
);
if
(
VERBOSE
)
fprintf
(
FP
,
"Write PUP data to file...
\n
"
);
/* write spectrum to file */
/* write spectrum to file */
if
(
ishot
==
1
&&
iter
==
1
)
{
if
(
ishot
==
1
&&
iter
==
1
)
{
...
@@ -200,7 +200,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -200,7 +200,7 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
fclose
(
file
);
fclose
(
file
);
}
}
fprintf
(
FP
,
"Write PUP-data into p-data array...
\n
"
);
if
(
VERBOSE
)
fprintf
(
FP
,
"Write PUP-data into p-data array...
\n
"
);
/* overwrite p data */
/* overwrite p data */
for
(
k
=
1
;
k
<=
ntr_glob
;
k
++
)
{
for
(
k
=
1
;
k
<=
ntr_glob
;
k
++
)
{
for
(
i
=
1
;
i
<=
ns
;
i
++
)
{
for
(
i
=
1
;
i
<=
ns
;
i
++
)
{
...
@@ -218,7 +218,6 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
...
@@ -218,7 +218,6 @@ void pup(float **data_p, float **data_vy, FILE *fp, int ntr_glob, int **recpos,
outseis_glob
(
fp
,
fopen
(
pf
,
"w"
),
0
,
data_p
,
recpos
,
recpos_loc
,
ntr
,
srcpos
,
1
,
ns
,
SEIS_FORMAT
,
ishot
,
1
);
outseis_glob
(
fp
,
fopen
(
pf
,
"w"
),
0
,
data_p
,
recpos
,
recpos_loc
,
ntr
,
srcpos
,
1
,
ns
,
SEIS_FORMAT
,
ishot
,
1
);
free_matrix
(
filtwk
,
1
,
ny_pot
/
2
,
1
,
nx_pot
/
2
);
free_matrix
(
filtwk
,
1
,
ny_pot
/
2
,
1
,
nx_pot
/
2
);
free_matrix
(
filtwk_full
,
1
,
ny_pot
,
1
,
nx_pot
);
free_matrix
(
filtwk_full
,
1
,
ny_pot
,
1
,
nx_pot
);
free_matrix
(
data_pim
,
1
,
ny_pot
,
1
,
nx_pot
);
free_matrix
(
data_pim
,
1
,
ny_pot
,
1
,
nx_pot
);
...
...
src/read_par_json.c
View file @
0803b755
...
@@ -408,7 +408,15 @@ void read_par_json(FILE *fp, char *fileinp){
...
@@ -408,7 +408,15 @@ void read_par_json(FILE *fp, char *fileinp){
else
{
else
{
if
(
SEISMO
<
4
&&
WAVESEP
==
1
)
{
if
(
SEISMO
<
4
&&
WAVESEP
==
1
)
{
WAVESEP
=
0
;
WAVESEP
=
0
;
fprintf
(
fp
,
"Variable SEISMO has to be set to 4 or 5 to use wavefield separation. WAVESEP is set now to %i.
\n
"
,
WAVESEP
);
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
)
{
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
(
WAVESEP
==
1
)
{
if
(
get_float_from_objectlist
(
"VEL"
,
number_readobjects
,
&
VEL
,
varname_list
,
value_list
)){
if
(
get_float_from_objectlist
(
"VEL"
,
number_readobjects
,
&
VEL
,
varname_list
,
value_list
)){
...
...
Write
Preview
Markdown
is supported
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment