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

Merge branch 'feature/streamer' into develop

parents a19a1cf8 a3da04c0
......@@ -324,6 +324,8 @@ These receiver coordinates in REC\_FILE are shifted by REFREC[1], REFREC[2] into
If READREC=0 the receiver locations must be specified in the parameter file. In this case, it is assumed that the receivers are located along a straight line. The first receiver position is defined by (XREC1, YREC1), and the last receiver position by (XREC1, YREC1). The spacing between receivers is NGEOPH grid points.
READREC=2 is made for a moving streamer aquisition geometry. The receiver positions for shot 1 must be specified in the parameter file (as described for READREC=0). For all other shots the receiver positions are moved according to the movement of the following source. This implementation allows a memory saving usage of moving geometries. The receivers must be in the same depth and all receivers must be located in the model for all shots.
Receivers are always located on full grid indices, i.e. a receiver that is located between two grid points will be shifted by the FD program to the closest next grid point. It is not yet possible
to output seismograms for arbitrary receiver locations since this would require a certain wavefield interpolation.
......
This diff is collapsed.
......@@ -70,6 +70,7 @@ IFOS2D= \
IFOS2D.c \
stf.c \
window_cos.c \
alloc_sections.c \
calc_mat_change_test.c \
calc_res.c \
calc_misfit.c \
......
/*-----------------------------------------------------------------------------------------
* 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>.
-----------------------------------------------------------------------------------------*/
/* Computation of local source coordinates
*
*
*/
#include "fd.h"
void alloc_sections(int ntr,int ns,float ***sectionvx,float ***sectionvy,float ***sectionvz,float ***sectionp,float ***sectionpnp1,float ***sectionpn,float ***sectioncurl,float ***sectiondiv,
float ***sectionpdata,float ***sectionpdiff,float ***sectionpdiffold,float ***sectionvxdata,float ***sectionvxdiff,float ***sectionvxdiffold,float ***sectionvydata,
float ***sectionvydiff,float ***sectionvydiffold,float ***sectionvzdata,float ***sectionvzdiff,float ***sectionvzdiffold)
{
extern int SEISMO, WAVETYPE;
extern FILE *FP;
switch (SEISMO){
case 1 : /* particle velocities only */
switch (WAVETYPE) {
case 1:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
*sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
*sectionvz=matrix(1,ntr,1,ns);
break;
}
break;
case 2 : /* pressure only */
*sectionp=matrix(1,ntr,1,ns);
*sectionpnp1=matrix(1,ntr,1,ns);
*sectionpn=matrix(1,ntr,1,ns);
break;
case 3 : /* curl and div only */
*sectioncurl=matrix(1,ntr,1,ns);
*sectiondiv=matrix(1,ntr,1,ns);
break;
case 4 : /* everything */
switch (WAVETYPE) {
case 1:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
*sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
*sectionvz=matrix(1,ntr,1,ns);
break;
}
*sectioncurl=matrix(1,ntr,1,ns);
*sectiondiv=matrix(1,ntr,1,ns);
*sectionp=matrix(1,ntr,1,ns);
break;
case 5 : /* everything except curl and div*/
switch (WAVETYPE) {
case 1:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
break;
case 2:
*sectionvz=matrix(1,ntr,1,ns);
break;
case 3:
*sectionvx=matrix(1,ntr,1,ns);
*sectionvy=matrix(1,ntr,1,ns);
*sectionvz=matrix(1,ntr,1,ns);
break;
}
*sectionp=matrix(1,ntr,1,ns);
break;
}
*sectionpdata=matrix(1,ntr,1,ns);
*sectionpdiff=matrix(1,ntr,1,ns);
*sectionpdiffold=matrix(1,ntr,1,ns);
switch (WAVETYPE) {
case 1:
*sectionvxdata=matrix(1,ntr,1,ns);
*sectionvxdiff=matrix(1,ntr,1,ns);
*sectionvxdiffold=matrix(1,ntr,1,ns);
*sectionvydata=matrix(1,ntr,1,ns);
*sectionvydiff=matrix(1,ntr,1,ns);
*sectionvydiffold=matrix(1,ntr,1,ns);
break;
case 2:
*sectionvzdata=matrix(1,ntr,1,ns);
*sectionvzdiff=matrix(1,ntr,1,ns);
*sectionvzdiffold=matrix(1,ntr,1,ns);
break;
case 3:
*sectionvxdata=matrix(1,ntr,1,ns);
*sectionvxdiff=matrix(1,ntr,1,ns);
*sectionvxdiffold=matrix(1,ntr,1,ns);
*sectionvydata=matrix(1,ntr,1,ns);
*sectionvydiff=matrix(1,ntr,1,ns);
*sectionvydiffold=matrix(1,ntr,1,ns);
*sectionvzdata=matrix(1,ntr,1,ns);
*sectionvzdiff=matrix(1,ntr,1,ns);
*sectionvzdiffold=matrix(1,ntr,1,ns);
break;
}
}
void dealloc_sections(int ntr,int ns,int **recpos_loc,float **sectionvx,float **sectionvy,float **sectionvz,float **sectionp,float **sectionpnp1,float **sectionpn,float **sectioncurl,float **sectiondiv,
float **sectionpdata,float **sectionpdiff,float **sectionpdiffold,float **sectionvxdata,float **sectionvxdiff,float **sectionvxdiffold,float **sectionvydata,
float **sectionvydiff,float **sectionvydiffold,float **sectionvzdata,float **sectionvzdiff,float **sectionvzdiffold)
{
extern int SEISMO, WAVETYPE;
extern FILE *FP;
free_imatrix(recpos_loc,1,3,1,ntr);
switch (SEISMO){
case 1 : /* particle velocities only */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvx,1,ntr,1,ns);
free_matrix(sectionvy,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvz,1,ntr,1,ns);
}
break;
case 2 : /* pressure only */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionp,1,ntr,1,ns);
free_matrix(sectionpn,1,ntr,1,ns);
free_matrix(sectionpnp1,1,ntr,1,ns);
}
break;
case 3 : /* curl and div only */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectioncurl,1,ntr,1,ns);
free_matrix(sectiondiv,1,ntr,1,ns);
}
break;
case 4 : /* everything */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvx,1,ntr,1,ns);
free_matrix(sectionvy,1,ntr,1,ns);
free_matrix(sectionp,1,ntr,1,ns);
free_matrix(sectioncurl,1,ntr,1,ns);
free_matrix(sectiondiv,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvz,1,ntr,1,ns);
}
break;
case 5 : /* everything except curl and div */
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvx,1,ntr,1,ns);
free_matrix(sectionvy,1,ntr,1,ns);
free_matrix(sectionp,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvz,1,ntr,1,ns);
}
break;
}
if (WAVETYPE==1 || WAVETYPE==3) {
free_matrix(sectionvxdata,1,ntr,1,ns);
free_matrix(sectionvxdiff,1,ntr,1,ns);
free_matrix(sectionvydata,1,ntr,1,ns);
free_matrix(sectionvydiff,1,ntr,1,ns);
free_matrix(sectionvydiffold,1,ntr,1,ns);
free_matrix(sectionvxdiffold,1,ntr,1,ns);
free_matrix(sectionpdata,1,ntr,1,ns);
free_matrix(sectionpdiff,1,ntr,1,ns);
free_matrix(sectionpdiffold,1,ntr,1,ns);
}
if (WAVETYPE==2 || WAVETYPE==3) {
free_matrix(sectionvzdata,1,ntr,1,ns);
free_matrix(sectionvzdiff,1,ntr,1,ns);
free_matrix(sectionvzdiffold,1,ntr,1,ns);
}
}
\ No newline at end of file
......@@ -277,7 +277,7 @@ void checkfd(FILE *fp, float ** prho, float ** ppi, float ** pu, float ** ptaus,
therefore we determine the minimum/maximum position in y-direction by the ZREC1 variable and vice versa.
this has to be considered for the receiver line coordinates specified in both the input file and separate source/receiver files*/
if (READREC==0) {
if (READREC==0 || READREC==2) {
if (XREC1>XREC2) {
srec_maxx=XREC1;
srec_minx=XREC2;
......
......@@ -47,6 +47,10 @@ void spat_filt(float ** waveconv, int iter, int sws);
float norm(float ** waveconv, int iter, int sws);
void alloc_sections(int ntr,int ns,float ***sectionvx,float ***sectionvy,float ***sectionvz,float ***sectionp,float ***sectionpnp1,float ***sectionpn,float ***sectioncurl,float ***sectiondiv,
float ***sectionpdata,float ***sectionpdiff,float ***sectionpdiffold,float ***sectionvxdata,float ***sectionvxdiff,float ***sectionvxdiffold,float ***sectionvydata,
float ***sectionvydiff,float ***sectionvydiffold,float ***sectionvzdata,float ***sectionvzdiff,float ***sectionvzdiffold);
void av_mat(float ** pi, float ** u,
float ** ppijm, float ** puip, float ** pujm);
......@@ -95,6 +99,10 @@ void count_killed_traces(int ntr, int swstestshot, int ntr_glob, int **recpos_lo
void create_trkill_table(int ** killtable, int ntr_glob, int **recpos, int nsrc_glob, float **srcpos, int ishot, float kill_offset_lower, float kill_offset_upper);
void dealloc_sections(int ntr,int ns,int **recpos_loc,float **sectionvx,float **sectionvy,float **sectionvz,float **sectionp,float **sectionpnp1,float **sectionpn,float **sectioncurl,float **sectiondiv,
float **sectionpdata,float **sectionpdiff,float **sectionpdiffold,float **sectionvxdata,float **sectionvxdiff,float **sectionvxdiffold,float **sectionvydata,
float **sectionvydiff,float **sectionvydiffold,float **sectionvzdata,float **sectionvzdiff,float **sectionvzdiffold);
float exchange_L2(float L2, int sw, int bcast_l2);
void exchange_rsg(float ** vx, float ** vy, float ** vz,
......@@ -226,7 +234,7 @@ void readmod_elastic(float ** rho, float ** pi, float ** u);
void readmod_elastic_es(float ** rho, float ** pi, float ** u, float ** matmod, int is);
int **receiver(FILE *fp, int *ntr);
int **receiver(int* ntr, float** srcpos, int shotno);
void save_checkpoint(int nx1, int nx2, int ny1, int ny2,
float ** vx, float ** vy, float ** sxx, float ** syy, float ** sxy);
......
......@@ -354,7 +354,7 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("READREC",number_readobjects,&READREC,varname_list, value_list))
declare_error("Variable READREC could not be retrieved from the json input file!");
else {
if (READREC==0) {
if (READREC==0 || READREC==2) {
if (get_float_from_objectlist("XREC1",number_readobjects,&XREC1,varname_list, value_list))
declare_error("Variable XREC1 could not be retrieved from the json input file!");
if (get_float_from_objectlist("XREC2",number_readobjects,&XREC2,varname_list, value_list))
......@@ -366,7 +366,7 @@ void read_par_json(FILE *fp, char *fileinp){
if (get_int_from_objectlist("NGEOPH",number_readobjects,&NGEOPH,varname_list, value_list))
declare_error("Variable NGEOPH could not be retrieved from the json input file!");
}
else {
if (READREC>0) {
if (get_string_from_objectlist("REC_FILE",number_readobjects,REC_FILE,varname_list, value_list))
declare_error("Variable REC_FILE could not be retrieved from the json input file!");
}
......@@ -1052,7 +1052,7 @@ void read_par_json(FILE *fp, char *fileinp){
}*/
/* receiver file */
if (READREC)
if (READREC==1)
{
if (access(REC_FILE,0) != 0)
{
......
......@@ -24,42 +24,44 @@
#include "fd.h"
int **receiver(FILE *fp, int *ntr){
int **receiver(int *ntr, float **srcpos, int shotno){
/* declaration of extern variables */
extern char REC_FILE[STRING_SIZE];
extern float DH, REFREC[4], REC_ARRAY_DEPTH, REC_ARRAY_DIST, FW;
extern int READREC, NGEOPH, DRX, REC_ARRAY, NX;
extern int MYID;
extern int READREC, NGEOPH, DRX, REC_ARRAY, NXG, NYG;
extern int MYID, VERBOSE,TRKILL;
extern FILE *FP;
int **recpos1, **recpos = NULL, nxrec=0, nyrec=0, nzrec=0;
int itr=1, itr1=0, itr2=0, recflag=0, c, ifw, n, i, j;
int nxrec1, nxrec2, nyrec1, nyrec2, nzrec1=1, nzrec2=1;
int **recpos1, **recpos = NULL, nxrec=0, nyrec=0, nzrec=0, recdist;
int itr=1, itr1=0, itr2=0, recflag=0, c, ifw, n, i, j, ntr1=0;
int nxrec1, nxrec2, nxrec3, nyrec1, nyrec2, nzrec1=1, nzrec2=1, offset, shotdist;
extern float XREC1, YREC1, XREC2, YREC2;
float xrec, yrec;
FILE *fpr;
FILE *fpr,*f;
char filename[STRING_SIZE];
if (MYID==0)
{
if (READREC){ /* read receiver positions from file */
fprintf(fp,"\n Reading receiver positions from file: \n\t%s\n",REC_FILE);
if (READREC==1){ /* read receiver positions from file */
fprintf(FP,"\n Reading receiver positions from file: \n\t%s\n",REC_FILE);
fpr=fopen(REC_FILE,"r");
if (fpr==NULL) declare_error(" Receiver file could not be opened !");
*ntr=0;
while ((c=fgetc(fpr)) != EOF)
if (c=='\n') ++(*ntr);
if (c=='\n') ++(ntr1);
rewind(fpr);
recpos1=imatrix(1,3,1,*ntr);
for (itr=1;itr<=*ntr;itr++){
recpos1=imatrix(1,3,1,ntr1);
for (itr=1;itr<=ntr1;itr++){
fscanf(fpr,"%f%f\n",&xrec, &yrec);
recpos1[1][itr]=iround((xrec+REFREC[1])/DH);
recpos1[2][itr]=iround((yrec+REFREC[2])/DH);
recpos1[3][itr]=iround((0.0+REFREC[3])/DH);
}
fclose(fpr);
fprintf(fp," Message from function receiver (written by PE %d):\n",MYID);/***/
fprintf(fp," Number of receiver positions found: %i\n",*ntr);
fprintf(FP," Message from function receiver (written by PE %d):\n",MYID);/***/
fprintf(FP," Number of receiver positions found: %i\n",ntr1);
/* recpos1[1][itr] X in m
* recpos1[2][itr] Y in m
......@@ -68,15 +70,15 @@ int **receiver(FILE *fp, int *ntr){
/* check if more than one receiver is located
at the same gridpoint */
for (itr=1;itr<=(*ntr-1);itr++)
for (itr1=itr+1;itr1<=*ntr;itr1++)
for (itr=1;itr<=(ntr1-1);itr++)
for (itr1=itr+1;itr1<=ntr1;itr1++)
if ((recpos1[1][itr]==recpos1[1][itr1])
&& (recpos1[2][itr]==recpos1[2][itr1])
&& (recpos1[3][itr]==recpos1[3][itr1]))
recpos1[1][itr1]=-(++recflag);
recpos=imatrix(1,3,1,*ntr-recflag);
for (itr=1;itr<=*ntr;itr++)
recpos=imatrix(1,3,1,ntr1-recflag);
for (itr=1;itr<=ntr1;itr++)
if (recpos1[1][itr]>0){
recpos[1][++itr2]=recpos1[1][itr];
recpos[2][itr2]=recpos1[2][itr];
......@@ -85,36 +87,73 @@ int **receiver(FILE *fp, int *ntr){
*ntr=itr2;
if (recflag>0){
fprintf(fp,"\n\n");
fprintf(fp," Warning:\n");
fprintf(fp," Several receivers located at the same gridpoint !\n");
fprintf(fp," Number of receivers reduced to %i\n", *ntr);
fprintf(fp,"\n\n");
fprintf(FP,"\n\n");
fprintf(FP," Warning:\n");
fprintf(FP," Several receivers located at the same gridpoint !\n");
fprintf(FP," Number of receivers reduced to %i\n", *ntr);
fprintf(FP,"\n\n");
}
free_imatrix(recpos1,1,3,1,*ntr);
free_imatrix(recpos1,1,3,1,ntr1);
}
else if (REC_ARRAY>0){
ifw=iround(FW/DH); /* frame width in gridpoints */
*ntr=(1+(NX-2*ifw)/DRX)*REC_ARRAY;
recpos=imatrix(1,3,1,*ntr);
itr=0;
for (n=0;n<=REC_ARRAY-1;n++){
j=iround((REC_ARRAY_DEPTH+REC_ARRAY_DIST*(float)n)/DH);
for (i=ifw;i<=NX-ifw;i+=DRX){
itr++;
recpos[1][itr]=i;
recpos[2][itr]=j;
else if (READREC==2){
fprintf(FP,"\n Moving streamer acquisition geometry choosen. \n");
/* decide if XREC1 or XREC2 is the nearest receiver to shot 1 */
if (abs(XREC1-srcpos[1][1])<abs(XREC2-srcpos[1][1])){
nxrec1=iround(XREC1/DH); /* (nxrec1,nyrec1) and (nxrec2,nyrec2) are */
nyrec1=iround(YREC1/DH); /* the positions of the first and last receiver*/
nxrec2=iround(XREC2/DH); /* in gridpoints */
nyrec2=iround(YREC2/DH);
offset=iround((XREC1-srcpos[1][1])/DH);
}
else {
nxrec1=iround(XREC2/DH); /* (nxrec1,nyrec1) and (nxrec2,nyrec2) are */
nyrec1=iround(YREC2/DH); /* the positions of the first and last receiver*/
nxrec2=iround(XREC1/DH); /* in gridpoints */
nyrec2=iround(YREC1/DH);
offset=iround((XREC2-srcpos[1][1])/DH);
}
if (nyrec1 != nyrec2){
fprintf(FP,"\n\n");
fprintf(FP," Warning:\n");
fprintf(FP," Receiver positions are not in the same depth !\n");
fprintf(FP," Depth of all receivers is seth to %i m\n",YREC1);
fprintf(FP,"\n\n");
}
if (offset<0) recdist=-NGEOPH;
else recdist=NGEOPH;
*ntr=iround((nxrec2-nxrec1)/recdist)+1;
if (shotno>0) nxrec1=nxrec1+((srcpos[1][shotno]-srcpos[1][1])/DH);
recpos=imatrix(1,3,1,*ntr);
n=0;
for (itr=1;itr<=*ntr;itr++){
nxrec=nxrec1+(itr-1)*recdist;
if (nxrec>FW && nxrec<(NXG-FW) && nyrec1>0 && nyrec1<(NYG-FW)){
recpos[1][itr]=nxrec;
recpos[2][itr]=nyrec1;
} else declare_error(" Receiver positions of current shot are out of model boundaries !");
}
if (shotno>0){
/* write receiver positions to file */
sprintf(filename,"%s.shot%i",REC_FILE,shotno);
f=fopen(filename,"wb");
for (i=1;i<=*ntr;i++) {
fprintf(f,"%f \t %f \n",recpos[1][i]*DH,recpos[2][i]*DH);
}
fclose(f);
}
if (VERBOSE==1) {
fprintf(FP," Message from function receiver (written by PE %d):\n",MYID);
fprintf(FP," Number of receiver positions found: %i\n",*ntr);
fprintf(FP," First receiver position for shot %i: %f (x), %f (y)\n",shotno,recpos[1][1]*DH,recpos[2][1]*DH);
fprintf(FP," Last receiver position for shot %i: %f (x), %f (y)\n",shotno,recpos[1][*ntr]*DH,recpos[2][*ntr]*DH);
}
}
else{ /* straight horizontal or vertical
line of receivers */
......
......@@ -184,7 +184,7 @@ void write_par(FILE *fp){
if (SEISMO){
fprintf(fp," ------------------------- RECEIVER --------------------------\n");
if (READREC){
if (READREC==1){
fprintf(fp," reading receiver positions from file \n");
fprintf(fp,"\t%s\n\n",REC_FILE);
fprintf(fp," reference_point_for_receiver_coordinate_system:\n");
......
Markdown is supported
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