Commit a98a59f2 authored by laura.gassner's avatar laura.gassner

fixed bugs in snapmerge.c and exchange_par.c, and made a few other minor adjustments.

parent e3375c7d
......@@ -22,34 +22,40 @@
*/
#include "fd.h"
#include "../src/fd.h"
void model_acoustic(float **rho, float **pi){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY, NXG, NYG, POS[3], MYID;
extern int NX, NY, NXG, NYG, POS[3], MYID, FDORDER;
extern float DH, DT;
extern char MFILE[STRING_SIZE];
extern char INV_MODELFILE[STRING_SIZE];
extern char TAPER_FILE_NAME[STRING_SIZE], TAPER_FILE_NAME_RHO[STRING_SIZE];
/* local variables */
float piv, vp, rhov;
int i, j, ii, jj, l;
float piv, vp, rhov, taperv, **taper;
int i, j, ii, jj, l, nd;
char modfile[STRING_SIZE];
FILE *flfile;
int nodes;
char cline[256];
float *fldepth, *flrho, *flvp;
float *fldepth, *flrho, *flvp, *fltaper;
nd = FDORDER/2 + 1;
/*read FL nodes from File*/
nodes=4;
nodes=5;
fldepth=vector(1,nodes);
flrho=vector(1,nodes);
flvp=vector(1,nodes);
fltaper=vector(1,nodes);
taper=matrix(-nd+1,NY+nd,-nd+1,NX+nd);
flfile=fopen("model_true/flnodes.toy_example_ac.start","r");
if (flfile==NULL) err(" FL-file could not be opened !");
......@@ -57,7 +63,7 @@ void model_acoustic(float **rho, float **pi){
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
sscanf(cline,"%f%f%f",&fldepth[l], &flrho[l], &flvp[l]);
sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &fltaper[l]);
}
else l=l-1;
......@@ -122,9 +128,53 @@ void model_acoustic(float **rho, float **pi){
}
}
/* loop over global grid - taper */
for (i=1;i<=NXG;i++){
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
if ((i==1) && (MYID==0)){
printf("depth: %f m: double node\n",fldepth[l]);}}
else{
for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
taperv=0.0;
taperv=(DH*(j-1)-fldepth[l])*(fltaper[l+1]-fltaper[l])/(fldepth[l+1]-fldepth[l])+fltaper[l];
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
taper[jj][ii]=taperv;
}
}
}
}
for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
taperv=0.0;
taperv=fltaper[nodes];
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
taper[jj][ii]=taperv;
}
}
}
free_vector(fldepth,1,nodes);
free_vector(flrho,1,nodes);
free_vector(flvp,1,nodes);
free_vector(fltaper,1,nodes);
/**************************************************/
......@@ -164,8 +214,9 @@ void model_acoustic(float **rho, float **pi){
for (l=1;l<nodes;l++){
if (fldepth[l]==fldepth[l+1]){
if ((i==1) && (MYID==0)){
printf("depth: %f m: double node\n",fldepth[l]);}}
else{
printf("depth: %f m: double node\n",fldepth[l]);
}
}else{
for (j=(int)(fldepth[l]/DH)+1;j<=(int)(fldepth[l+1]/DH);j++){
rhov=0.0;
......@@ -176,13 +227,11 @@ void model_acoustic(float **rho, float **pi){
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
rho[jj][ii]=rhov;
// pi[jj][ii]=piv;
}
}
}
......@@ -196,8 +245,7 @@ void model_acoustic(float **rho, float **pi){
/* only the PE which belongs to the current global gridpoint
is saving model parameters in his local arrays */
if ((POS[1]==((i-1)/NX)) &&
(POS[2]==((j-1)/NY))){
if ((POS[1]==((i-1)/NX)) && (POS[2]==((j-1)/NY))){
ii=i-POS[1]*NX;
jj=j-POS[2]*NY;
......@@ -227,4 +275,22 @@ void model_acoustic(float **rho, float **pi){
sprintf(modfile,"%s_vp_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s",TAPER_FILE_NAME);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s",TAPER_FILE_NAME_RHO);
writemod(modfile,taper,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s.%i%i",TAPER_FILE_NAME_RHO,POS[1],POS[2]);
remove(modfile);
free_matrix(taper,-nd+1,NY+nd,-nd+1,NX+nd);
}
......@@ -32,15 +32,11 @@
"QUELLTYP" : "3",
"QUELLTYP values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4" : "comment",
"SRCREC" : "1",
"SRCREC values : read source positions from SOURCE_FILE=1, PLANE_WAVE=2" : "comment",
"SOURCE_FILE" : "./source/sources.dat",
"RUN_MULTIPLE_SHOTS" : "1",
"PLANE_WAVE_DEPTH" : "0.0",
"PHI" : "0.0",
"TS" : "0.032",
"Acoustic Computation" : "comment",
"ACOUSTIC" : "1",
"Model" : "comment",
"READMOD" : "0",
......@@ -92,22 +88,10 @@
"SNAP_FORMAT" : "3",
"SNAP_FILE" : "./snap/waveform_forward",
"Receiver array" : "comment",
"REC_ARRAY" : "0",
"REC_ARRAY_DEPTH" : "70.0",
"REC_ARRAY_DIST" : "40.0",
"DRX" : "4",
"Monitoring the simulation" : "comment",
"LOG_FILE" : "log/2layer.log",
"LOG" : "1",
"Checkpoints" : "comment",
"CHECKPTREAD" : "0",
"CHECKPTWRITE" : "0",
"CHECKPTFILE" : "tmp/checkpoint_fdveps",
"General inversion parameters" : "comment",
"INVMAT1" : "1",
......
......@@ -32,15 +32,11 @@
"QUELLTYP" : "3",
"QUELLTYP values (point_source): explosive=1;force_in_x=2;force_in_y=3;rotated_force=4" : "comment",
"SRCREC" : "1",
"SRCREC values : read source positions from SOURCE_FILE=1, PLANE_WAVE=2" : "comment",
"SOURCE_FILE" : "./source/sources.dat",
"RUN_MULTIPLE_SHOTS" : "1",
"PLANE_WAVE_DEPTH" : "0.0",
"PHI" : "0.0",
"TS" : "0.032",
"Acoustic Computation" : "comment",
"ACOUSTIC" : "1",
"Model" : "comment",
"READMOD" : "0",
......@@ -92,29 +88,16 @@
"SNAP_FORMAT" : "3",
"SNAP_FILE" : "./snap/waveform_forward",
"Receiver array" : "comment",
"REC_ARRAY" : "0",
"REC_ARRAY_DEPTH" : "70.0",
"REC_ARRAY_DIST" : "40.0",
"DRX" : "4",
"Monitoring the simulation" : "comment",
"LOG_FILE" : "log/2layer.log",
"LOG" : "1",
"Checkpoints" : "comment",
"CHECKPTREAD" : "0",
"CHECKPTWRITE" : "0",
"CHECKPTFILE" : "tmp/checkpoint_fdveps",
"General inversion parameters" : "comment",
"ITERMAX" : "10",
"DATA_DIR" : "su/measured_data/DENISE_real",
"INVMAT1" : "1",
"INVMAT" : "0",
"INVTYPE" : "2",
"QUELLTYPB" : "1",
"MISFIT_LOG_FILE" : "L2_LOG.dat",
"VELOCITY" : "0",
......@@ -173,10 +156,6 @@
"Minimum number of iteration per frequency" : "comment",
"MIN_ITER" : "0",
"Cosine taper" : "comment",
"TAPER" : "0",
"TAPERLENGTH" : "10",
"Definition of gradient taper geometry" : "comment",
"SWS_TAPER_GRAD_VERT" : "0",
"SWS_TAPER_GRAD_HOR" : "0",
......
......@@ -76,7 +76,7 @@
"TAU" : "0.0966",
"General inversion parameters" : "comment",
"ITERMAX" : "1",
"ITERMAX" : "100",
"DATA_DIR" : "su/measured_data/toy_example_ac",
"INVMAT1" : "1",
"INVMAT" : "0",
......@@ -124,10 +124,9 @@
"SRTRADIUS" : "3.0",
"FILTSIZE" : "1",
"SWS_TAPER_FILE" : "1",
"TAPER_FILE_NAME" : "taper.bin",
"TAPER_FILE_NAME_U" : "taper_u.bin",
"TAPER_FILE_NAME_RHO" : "taper_rho.bin",
"TAPER_FILE_NAME" : "taper/taper.bin",
"TAPER_FILE_NAME_U" : "taper/taper_u.bin",
"TAPER_FILE_NAME_RHO" : "taper/taper_rho.bin",
"Upper and lower limits for model parameters" : "comment",
"VPUPPERLIM" : "5000",
......@@ -156,4 +155,16 @@
"FC_END" : "70.0",
"FC_INCR" : "10.0",
"ORDER" : "4",
"ZERO_PHASE" : "1",
"Trace killing" : "comment",
"TRKILL" : "0",
"TRKILL_FILE" : "./trace_kill/trace_kill.dat",
"Time windowing and damping" : "comment",
"TIMEWIN" : "0",
"PICKS_FILE" : "./picked_times/PickedTimes",
"TWLENGTH_PLUS" : "4.0",
"TWLENGTH_MINUS" : "0.01",
"GAMMA" : "100000",
}
#FLNODE
#converted by MOCOX V1.2 from final.mod.gremlin
# 0.0000 1 0 6371.0000
# depth rho vp
# depth rho vp taper
# in m in g/cm^3 in km/s
#4
0.0000 1.0200 1.4800
500.0000 1.0200 1.4800
500.0000 2.0010 1.5000
2000.0000 2.0010 2.5000
#5
0.0000 1.0200 1.4800 0
500.0000 1.0200 1.4800 0
500.0000 2.0010 1.5000 0
550.0000 2.0010 1.5000 1
2000.0000 2.0010 2.5000 1
......@@ -40,5 +40,7 @@ mpirun -np 4 nice -19 ../bin/denise in_and_out/toy_example/toy_example_ac_INV.js
make clean
# rm jacobian/toy_example/*.old.*.*
# rm model/toy_example/*.bin.*.*
rm jacobian/toy_example/*toy_example*.*.*.*
rm model/*.bin.*.*
rm model/toy_example/*.bin.*.*
rm taper/*.bin.*.*
\ No newline at end of file
......@@ -38,13 +38,15 @@ float intseis_s, intseis_sd;
float *picked_times=NULL;
float **intseis_section=NULL, **intseis_sectiondata=NULL;
float **intseis_sectiondata_envelope=NULL, **intseis_section_envelope=NULL, **intseis_section_hilbert=NULL, **dummy_1=NULL, **dummy_2=NULL;
if(LNORM==8){
intseis_section_envelope = matrix(1,ntr,1,ns);
intseis_sectiondata_envelope = matrix(1,ntr,1,ns);
intseis_section_hilbert = matrix(1,ntr,1,ns);
dummy_1 = matrix(1,ntr,1,ns);
dummy_2 = matrix(1,ntr,1,ns);
}
}
intseis_section = matrix(1,ntr,1,ns); /* declaration of variables for integration */
intseis_sectiondata = matrix(1,ntr,1,ns);
if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
......@@ -53,13 +55,12 @@ if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWI
umax=ntr*ns;
zero(&sectiondiff[1][1],umax);
/* TRACE KILLING */
int ** kill_tmp, *kill_vector; /* declaration of variables for trace killing */
char trace_kill_file[STRING_SIZE];
FILE *ftracekill;
if(TRKILL==1){
if(TRKILL==1){
/* sectiondiff will be set to zero */
umax=ntr*ns;
zero(&sectiondiff[1][1],umax);
......@@ -86,16 +87,13 @@ if(TRKILL==1){
}
} /* end if(TRKILL)*/
RMS=0.0;
Lcount=1;
/* Integration of measured and synthetic data */
for(i=1;i<=ntr;i++){
intseis_s=0;
intseis_sd=0;
if (VELOCITY==0){ /* only integtration if displacement seismograms are inverted */
for(j=1;j<=ns;j++){
intseis_s += section[i][j];
......@@ -113,25 +111,22 @@ for(i=1;i<=ntr;i++){
/* TIME WINDOWING */
if(TIMEWIN==1){
time_window(intseis_section, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(intseis_sectiondata, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(intseis_section, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(intseis_sectiondata, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
}
/* NORMALIZE TRACES */
if(NORMALIZE==1){
normalize_data(intseis_section,ntr,ns);
normalize_data(intseis_sectiondata,ntr,ns);
normalize_data(intseis_section,ntr,ns);
normalize_data(intseis_sectiondata,ntr,ns);
}
/* calculate RMS */
/*for(i=1;i<=ntr;i++){
for(j=1;j<=ns;j++){*/
/*RMS += (sectionpdata[i][j]-sectionvx[i][j]) * (sectionpdata[i][j]-sectionvx[i][j]);*/
/*RMS += (sectionpdata[i][j]) * (sectionpdata[i][j]);*/
/* Lcount++;
/* Lcount++;
}
} */
......@@ -139,10 +134,8 @@ normalize_data(intseis_sectiondata,ntr,ns);
RMS=1.0;
/* calculate weighted data residuals and reverse time direction */
/* calculate kind of "energy" */
/* calculate envelope and hilbert transform for LNORM==8 */
if (LNORM==8){
calc_envelope(intseis_sectiondata,intseis_sectiondata_envelope,ns,ntr);
......@@ -158,7 +151,6 @@ if (LNORM==8){
/* end of calculate envelope and hilbert transform for LNORM==8 */
for(i=1;i<=ntr;i++){
if((TRKILL==1)&&(kill_vector[i]==1))
continue;
......@@ -175,7 +167,6 @@ for(i=1;i<=ntr;i++){
abs_sectiondata=sqrt(abs_sectiondata);
abs_section=sqrt(abs_section);
}
/* calculate residual seismograms and norm */
if((TRKILL==1)&&(kill_vector[i]==1))
......@@ -185,10 +176,7 @@ for(i=1;i<=ntr;i++){
invtime=ns;
for(j=1;j<=ns;j++){
/*printf("%d \t %d \t %e \t %e \n",i,j,sectionpdata[i][j],sectionp[i][j]);*/
/* calculate L1 residuals */
if(LNORM==1){
if(((sectiondata[i][j]-section[i][j])/RMS)>0){signL1=1.0;}
......@@ -225,7 +213,6 @@ for(i=1;i<=ntr;i++){
/*sectionpdiff[i][invtime]=sectionp[i][j];*/
/* calculate norm for different LNORM */
if((LNORM==2) && (swstestshot==1)){
L2+=sectiondiff[i][invtime]*sectiondiff[i][invtime];
......@@ -245,7 +232,6 @@ for(i=1;i<=ntr;i++){
/*L2+=sectiondiff[i][invtime];*/
invtime--; /* reverse time direction */
}
}
......@@ -253,7 +239,6 @@ for(i=1;i<=ntr;i++){
l2=L2;
// taper(sectiondiff,ntr,ns);
/* printf("\n MYID = %i IN CALC_RES: L2 = %10.12f \n",MYID,l2); */
/* free memory for integrated seismograms */
......@@ -262,9 +247,12 @@ free_matrix(intseis_sectiondata,1,ntr,1,ns);
/* free memory for time windowing and trace killing */
if(TIMEWIN==1) free_vector(picked_times,1,ntr);
if(TRKILL==1){
free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
free_ivector(kill_vector,1,ntr);}
free_imatrix(kill_tmp,1,ntr_glob,1,nsrc_glob);
free_ivector(kill_vector,1,ntr);
}
if(LNORM==8){
free_matrix(intseis_sectiondata_envelope,1,ntr,1,ns);
free_matrix(intseis_section_envelope,1,ntr,1,ns);
......
This diff is collapsed.
......@@ -69,7 +69,7 @@ void exchange_par(void){
extern int STEPMAX;
extern float EPS_SCALE, SCALEFAC;
extern float PRO;
extern int TRKILL;
extern int TRKILL, TRKILL_STF;
extern char TRKILL_FILE[STRING_SIZE], TRKILL_FILE_STF[STRING_SIZE];
extern int TIMEWIN, NORMALIZE;
extern float TWLENGTH_PLUS, TWLENGTH_MINUS, GAMMA;
......@@ -192,7 +192,7 @@ void exchange_par(void){
idum[18] = SRCREC;
idum[19] = IDX;
idum[20] = IDY;
idum[21] = 0;
idum[21] = TRKILL_STF;
idum[22] = 0;
idum[23] = SNAP_FORMAT;
idum[24] = SEISMO;
......@@ -426,7 +426,7 @@ void exchange_par(void){
SRCREC = idum[18];
IDX = idum[19];
IDY = idum[20];
TRKILL_STF=idum[21];
SNAP_FORMAT = idum[23];
SEISMO = idum[24];
......
......@@ -70,10 +70,10 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
sprintf(data,"%s_p.su.shot%d",DATA_DIR,comp);
}
/*printf("%s\n",data);*/
fpdata = fopen(data,"r");
if (fpdata==NULL) err(" Seismograms for inversion were not found ");
/* declaration of local variables */
int i,j;
......@@ -93,9 +93,6 @@ void inseis(FILE *fp, int comp, float **section, int ntr, int ns, int sws, int
dump=tr.data[j];
section[tracl1][j+1]=dump;
}
}
fclose(fpdata);
}
......@@ -24,8 +24,7 @@
*/
#include "fd.h"
int **splitrec(int **recpos,int *ntr_loc, int ntr, int *recswitch)
{
int **splitrec(int **recpos,int *ntr_loc, int ntr, int *recswitch){
extern int IENDX, IENDY, POS[4], REC1, REC2;
......
......@@ -30,8 +30,8 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
/* declaration of global variables */
extern float DT, DH;
extern int SEIS_FORMAT, MYID, NT, QUELLART, TIME_FILT;
extern char SEIS_FILE_VY[STRING_SIZE], PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
extern int SEIS_FORMAT, MYID, NT, QUELLART, TIME_FILT, TIMEWIN;
extern char SEIS_FILE_VY[STRING_SIZE], SEIS_FILE_P[STRING_SIZE], PARA[STRING_SIZE], DATA_DIR[STRING_SIZE];
extern int TRKILL_STF, NORMALIZE;
extern char TRKILL_FILE_STF[STRING_SIZE];
extern char SIGNAL_FILE[STRING_SIZE];
......@@ -41,20 +41,25 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
char trace_kill_file[STRING_SIZE];
FILE *ftracekill;
float *picked_times=NULL;
/* --------------- declaration of variables --------------- */
unsigned int nrec, nsamp, i, npairs;
float dt;
float xr=0.0, yr=0.0;
float XS=0.0, YS=0.0;
char conv_y[STRING_SIZE], qw[STRING_SIZE], conv_y_tmp[STRING_SIZE];
char conv_y[STRING_SIZE], qw[STRING_SIZE], conv_y_tmp[STRING_SIZE], obs_y_tmp[STRING_SIZE], mod_y_tmp[STRING_SIZE];
/* variables for wavelet */
int nt, nts;
float tshift, amp=0.0, fc, tau, t, ts, ag;
float * wavelet, * stf_conv_wavelet, *psource=NULL;
wavelet=vector(1,ns);
stf_conv_wavelet=vector(1,ns);
if(TIMEWIN) picked_times = vector(1,ntr); /* declaration of variables for TIMEWIN */
printf("\n================================================================================================\n\n");
printf("\n ***** Inversion of Source Time Function - shot: %d - it: %d ***** \n\n",ishot,iter);
......@@ -94,10 +99,14 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
}
if(i==ntr_glob)printf(" ***** \n\n");
}
}
/* trace killing ends here */
if(TIMEWIN==1){
time_window(sectionvy, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
time_window(sectionvy_obs, picked_times, iter, ntr_glob,recpos_loc, ntr, ns, ishot);
}
/* NORMALIZE TRACES */
/*if(NORMALIZE==1){*/
normalize_data(sectionvy,ntr_glob,ns);
......@@ -113,6 +122,7 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
XS=srcpos[1][ishot];
YS=srcpos[2][ishot];
/* TF Software: see libstfinv */
struct CTriples data;
......@@ -146,20 +156,19 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
stf.sampling.dt=dt;
/*char para[]="fbd:tshift=0.0"; /* parameter string */
initstfinvengine(data, stf, PARA);
runstfinvengine();
/* END TF Software */
if (QUELLART==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r"));
if (QUELLART==7){
psource=vector(1,ns);
inseis_source_wavelet(psource,ns,ishot);}
psource=vector(1,ns);
if (QUELLART==3) psource=rd_sour(&nts,fopen(SIGNAL_FILE,"r"));
if (QUELLART==7){
inseis_source_wavelet(psource,ns,ishot);
}
/* calculating wavelet SIN**3 for convoling with STF */
tshift=srcpos[4][ishot];
......@@ -234,19 +243,23 @@ int ntr_glob,int ntr, float ** srcpos, int ishot, int ns, int iter, int nshots,
}/* end of for (nt=1;nt<=ns;nt++) */
/* convolving wavelet with STF */
conv_FD(wavelet,source_time_function,stf_conv_wavelet,ns);