Commit 5020ca0c authored by Florian Wittkamp's avatar Florian Wittkamp

Merge branch 'develop' into Seitosh/import

parents 86dd3cb5 939652b7
stages:
- build
gcc:
stage: build
script:
- export cc=gcc && ${cc} --version
- export OMPI_MPICC=${cc} && mpicc --version
- cd par/ && make install
\ No newline at end of file
# What is IFOS2D? # What is IFOS2D?
**IFOS2D** (**I**nversion of **F**ull **O**bserved **S**eismograms) is a 2-D elastic full waveform inversion code. **IFOS2D** (**I**nversion of **F**ull **O**bserved **S**eismograms) is a 2-D elastic full waveform inversion code.
The inversion problem is solved by a conjugate gradient method and the gradients are computed in the time domain by the adjoint state method. The inversion problem is solved by a conjugate gradient method and the gradients are computed in the time domain by the adjoint state method.
The forward modeling is done by a time domain Finite-Difference scheme. The forward modeling is done by a time domain Finite-Difference scheme.
IFOS2D is the reverse (inverse) of our 2-D Finite-Difference forward solver [**SOFI2D**](https://git.scc.kit.edu/GPIAG-Software/SOFI2D). IFOS2D is the reverse (inverse) of our 2-D Finite-Difference forward solver [**SOFI2D**](https://git.scc.kit.edu/GPIAG-Software/SOFI2D).
The [**manual**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/wikis/home) is included in the download archive or can be downloaded [here](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/wikis/home). The [**manual**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/wikis/home) is included in the download archive or can be downloaded [here](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/wikis/home).
# Download and Newsletter # Download and Newsletter
Release: [![build status](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/badges/master/build.svg)](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/commits/master) Beta: [![build status](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/badges/develop/build.svg)](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/commits/develop)
You can download the [**latest Release 2.0.2**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tags/Release_2.0.2) or the current [**Beta-Version**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tree/develop). You can download the [**latest Release 2.0.2**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tags/Release_2.0.2) or the current [**Beta-Version**](https://git.scc.kit.edu/GPIAG-Software/IFOS2D/tree/develop).
To receive news and updates please [register](http://www.gpi.kit.edu/Software-FWI.php) on the email list [IFOS@lists.kit.edu](http://www.gpi.kit.edu/Software-FWI.php). To receive news and updates please [register](http://www.gpi.kit.edu/Software-FWI.php) on the email list [IFOS@lists.kit.edu](http://www.gpi.kit.edu/Software-FWI.php).
Please use this list also to ask questions or to report problems or bugs. Please use this list also to ask questions or to report problems or bugs.
\ No newline at end of file
...@@ -322,7 +322,9 @@ The locations of the receivers may either be specified in a separate file REC\_F ...@@ -322,7 +322,9 @@ The locations of the receivers may either be specified in a separate file REC\_F
\end{verbatim}}} \end{verbatim}}}
These receiver coordinates in REC\_FILE are shifted by REFREC[1], REFREC[2] into the x- and y-direction, respectively. This allows for completely moving the receiver spread without modifying REC\_FILE. This may be useful for the simulation of moving profiles in reflection seismics. These receiver coordinates in REC\_FILE are shifted by REFREC[1], REFREC[2] into the x- and y-direction, respectively. This allows for completely moving the receiver spread without modifying REC\_FILE. This may be useful for the simulation of moving profiles in reflection seismics.
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. 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 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. to output seismograms for arbitrary receiver locations since this would require a certain wavefield interpolation.
......
/*-----------------------------------------------------------------------------------------
* Copyright (C) 2005 <name of author>
*
* This file is part of DENISE.
*
* DENISE 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.
*
* DENISE 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 DENISE. See file COPYING and/or <http://www.gnu.org/licenses/gpl-2.0.html>.
-----------------------------------------------------------------------------------------*/
/* $Id: hh_el.c,v 2.3 2007/08/21 13:16:19 tbohlen Exp $*/
/*
* Model defined by flnode file. In this version a constant Q-model with Q_p=Q_s is used. It
* is defined via the relaxation frequencies and the TAU value given in the input file.
*/
#include "fd.h"
void model(float ** rho, float ** pi, float ** u, float ** taus, float ** taup, float * eta){
/*--------------------------------------------------------------------------*/
/* extern variables */
extern int NX, NY, NXG, NYG, POS[3], MYID, L;
extern float DH, DT, *FL, TAU;
extern char MFILE[STRING_SIZE];
extern char INV_MODELFILE[STRING_SIZE];
/* local variables */
float muv, piv, vp, vs, rhov, ts, tp, *pts;
int i, j, ii, jj, l;
char modfile[STRING_SIZE];
FILE *flfile;
int nodes;
char cline[256];
float *fldepth, *flrho, *flvp, *flvs;
/*---------------------------------------
Trench "Ettlinger Linie" units in m
----------------------------------------*/
/* trench=1 - elliptic shaped trench , trench=2 - V shaped trench */
int trench=2;
float width=10;
float depth=3;
float pos=NXG/2*DH;
// Parameters of trench "Ettlinger Linie" velocities in km/s and density in kg/m^3
float evs=0.0917;
float evp=0.2569;
float erho=1.5;
/* fopen FL-node file*/
flfile=fopen("model_true/flnodes.toy_example","r");
if (flfile==NULL) declare_error(" FL-file could not be opened !");
/* nodes specified in FL-node File*/
nodes=7;
/*-----------No changes below that line needed ---------------------------*/
/* 2 linear lines as border for the v-shape */
float m1=2*depth/width;
float m2=-m1;
float c1=depth*(1-(2*pos/width));
float c2=depth*(1+(2*pos/width));
fldepth=vector(1,nodes);
flrho=vector(1,nodes);
flvp=vector(1,nodes);
flvs=vector(1,nodes);
/* vector for maxwellbodies */
pts=vector(1,L);
for (l=1;l<=L;l++) {
pts[l]=1.0/(2.0*PI*FL[l]);
eta[l]=DT/pts[l];
}
/*read FL nodes from File*/
for (l=1;l<=nodes;l++){
fgets(cline,255,flfile);
if (cline[0]!='#'){
sscanf(cline,"%f%f%f%f",&fldepth[l], &flrho[l], &flvp[l], &flvs[l]);
}
else l=l-1;
}
if(MYID==0){
printf(" ------------------------------------------------------------------ \n\n");
printf(" Information of FL nodes: \n\n");
printf(" \t depth \t rho \t vp \t vs \n\n");
for (l=1;l<=nodes;l++){
printf(" \t %f \t %f \t %f \t %f\n\n",fldepth[l],flrho[l],flvp[l],flvs[l]);
}
printf(" ------------------------------------------------------------------ \n\n");
}
/*-----------------------------------------------------------------------*/
/*-----------------------------------------------------------------------*/
/* loop over global grid */
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++){
vp=0.0; vs=0.0; rhov=0.0;
vp=(DH*(j-1)-fldepth[l])*(flvp[l+1]-flvp[l])/(fldepth[l+1]-fldepth[l])+flvp[l];
vp=vp*1000.0;
vs=(DH*(j-1)-fldepth[l])*(flvs[l+1]-flvs[l])/(fldepth[l+1]-fldepth[l])+flvs[l];
vs=vs*1000.0;
rhov=(DH*(j-1)-fldepth[l])*(flrho[l+1]-flrho[l])/(fldepth[l+1]-fldepth[l])+flrho[l];
rhov=rhov*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* 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;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
}
}
for (j=(int)(fldepth[nodes]/DH)+1;j<=NYG;j++){
vp=0.0; vs=0.0; rhov=0.0;
vp=flvp[nodes]*1000.0; vs=flvs[nodes]*1000.0; rhov=flrho[nodes]*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* 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;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
/*Trench Ettlinger Linie*/
if (trench==1) {
for (j=1;j<=(int) depth/DH+1;j++){
if ( ( ((pos-i*DH) * (pos-i*DH))/(width/2*width/2) + ((j*DH)*(j*DH))/(depth*depth) ) <=1){
vp=evp*1000.0;
vs=evs*1000.0;
rhov=erho*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* 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;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
}
}
/*Trench V-shape Ettlinger Linie*/
if (trench==2) {
for (j=1;j<=(int) depth/DH+1;j++){
if ( ((i*DH)*m1+c1 >=0) && j*DH <= (i*DH)*m1+c1 && j*DH <= (i*DH)*m2+c2){
vp=evp*1000.0;
vs=evs*1000.0;
rhov=erho*1000.0;
muv=vs;
piv=vp;
ts=TAU;
tp=TAU;
/* 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;
u[jj][ii]=muv;
rho[jj][ii]=rhov;
pi[jj][ii]=piv;
taus[jj][ii]=ts;
taup[jj][ii]=tp;
}
}
}
}
}
sprintf(modfile,"%s_rho_it0.bin",INV_MODELFILE);
writemod(modfile,rho,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_rho_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_vs_it0.bin",INV_MODELFILE);
writemod(modfile,u,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_vs_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_vp_it0.bin",INV_MODELFILE);
writemod(modfile,pi,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_vp_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_taus_it0.bin",INV_MODELFILE);
writemod(modfile,taus,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_taus_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
sprintf(modfile,"%s_taup_it0.bin",INV_MODELFILE);
writemod(modfile,taup,3);
MPI_Barrier(MPI_COMM_WORLD);
if (MYID==0) mergemod(modfile,3);
MPI_Barrier(MPI_COMM_WORLD);
sprintf(modfile,"%s_taup_it0.bin.%i%i",INV_MODELFILE,POS[1],POS[2]);
remove(modfile);
free_vector(fldepth,1,nodes);
free_vector(flrho,1,nodes);
free_vector(flvp,1,nodes);
free_vector(flvs,1,nodes);
free_vector(pts,1,L);
}
...@@ -7,6 +7,14 @@ ...@@ -7,6 +7,14 @@
# compiling all libraries and IFOS # compiling all libraries and IFOS
make install MODEL=../genmod/toy_example_true.c make install MODEL=../genmod/toy_example_true.c
#-------------------------------------------------------------------------------------
# Insert the trench "Ettlinger Linie" into the toy_example model
# in the file ../genmod/toy_example_true_trench.c the size and position can be adjusted
# Either an elliptic trench or a v-shaped trench can be choosen
#make install MODEL=../genmod/toy_example_true_trench.c
#---------------------------------------------------------------------------------------
# starting IFOS for forward modeling # starting IFOS for forward modeling
# (depending on the MPI package you maybe have to adjust the following programme call, # (depending on the MPI package you maybe have to adjust the following programme call,
# e.g. when using openMPI you do not have to use the 'lamboot' command) # e.g. when using openMPI you do not have to use the 'lamboot' command)
......
...@@ -53,7 +53,7 @@ int main(int argc, char **argv){ ...@@ -53,7 +53,7 @@ int main(int argc, char **argv){
int sum_killed_traces=0, sum_killed_traces_testshots=0, killed_traces=0, killed_traces_testshots=0; int sum_killed_traces=0, sum_killed_traces_testshots=0, killed_traces=0, killed_traces_testshots=0;
int *ptr_killed_traces=&killed_traces, *ptr_killed_traces_testshots=&killed_traces_testshots; int *ptr_killed_traces=&killed_traces, *ptr_killed_traces_testshots=&killed_traces_testshots;
float energy, energy_sum, energy_all_shots, energy_sum_all_shots; float energy, energy_sum, energy_all_shots, energy_sum_all_shots = 0.0;
float energy_SH, energy_sum_SH, energy_all_shots_SH, energy_sum_all_shots_SH; float energy_SH, energy_sum_SH, energy_all_shots_SH, energy_sum_all_shots_SH;
float L2_SH, L2sum_SH, L2_all_shots_SH, L2sum_all_shots_SH; float L2_SH, L2sum_SH, L2_all_shots_SH, L2sum_all_shots_SH;
...@@ -157,9 +157,15 @@ int main(int argc, char **argv){ ...@@ -157,9 +157,15 @@ int main(int argc, char **argv){
int nfrq=0; int nfrq=0;
int FREQ_NR=1; int FREQ_NR=1;
float JOINT_EQUAL_PSV=0.0, JOINT_EQUAL_SH=0.0;
float JOINT_EQUAL_PSV_all=0.0, JOINT_EQUAL_SH_all=0.0;
int JOINT_EQUAL_new_max=1;
FILE *fprec, *FPL2; FILE *fprec, *FPL2;
FILE *FPL2_JOINT;
char L2_joint_log[STRING_SIZE];
/* General parameters */ /* General parameters */
int nt_out; int nt_out;
...@@ -237,10 +243,14 @@ int main(int argc, char **argv){ ...@@ -237,10 +243,14 @@ int main(int argc, char **argv){
/* In the following, NX and NY denote size of the local grid ! */ /* In the following, NX and NY denote size of the local grid ! */
NX = IENDX; NX = IENDX;
NY = IENDY; NY = IENDY;
/* Reading source positions from SOURCE_FILE */
srcpos=sources(&nsrc);
nsrc_glob=nsrc;
ishot=0;
if (SEISMO){ if (SEISMO){
recpos=receiver(FP, &ntr); recpos=receiver(&ntr, srcpos, ishot);
recswitch = ivector(1,ntr); recswitch = ivector(1,ntr);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch); recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr; ntr_glob=ntr;
...@@ -705,105 +715,14 @@ int main(int argc, char **argv){ ...@@ -705,105 +715,14 @@ int main(int argc, char **argv){
} }
if (ntr>0){ if (ntr>0){
switch (SEISMO){ alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
case 1 : /* particle velocities only */ &sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
switch (WAVETYPE) { &sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
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;
}
} }
/* Memory for seismic data */ /* Memory for seismic data */
sectionread=matrix(1,ntr_glob,1,ns); sectionread=matrix(1,ntr_glob,1,ns);
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;
}
/* Memory for inversion for source time function */ /* Memory for inversion for source time function */
if((INV_STF==1)||(TIME_FILT==1) || (TIME_FILT==2)){ if((INV_STF==1)||(TIME_FILT==1) || (TIME_FILT==2)){
sectionp_conv=matrix(1,ntr_glob,1,NT); sectionp_conv=matrix(1,ntr_glob,1,NT);
...@@ -851,10 +770,6 @@ int main(int argc, char **argv){ ...@@ -851,10 +770,6 @@ int main(int argc, char **argv){
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
/* Reading source positions from SOURCE_FILE */
srcpos=sources(&nsrc);
nsrc_glob=nsrc;
if(FORWARD_ONLY==0&&USE_WORKFLOW){ if(FORWARD_ONLY==0&&USE_WORKFLOW){
read_workflow(FILE_WORKFLOW,&workflow, &workflow_lines,workflow_header); read_workflow(FILE_WORKFLOW,&workflow, &workflow_lines,workflow_header);
} }
...@@ -1072,27 +987,41 @@ int main(int argc, char **argv){ ...@@ -1072,27 +987,41 @@ int main(int argc, char **argv){
C_rho = rho_avg*rho_avg; C_rho = rho_avg*rho_avg;
} }
/* Seperate PSV and SH logging in case of a joint inversion */
if(WAVETYPE==3){
sprintf(L2_joint_log,"%s_JOINT",MISFIT_LOG_FILE);
}
/* Open Log File for L2 norm */ /* Open Log File for L2 norm */
if(FORWARD_ONLY!=1){ if(!FORWARD_ONLY && MYID==0){
if(MYID==0){
if(iter==1){ if(iter==1){
FPL2=fopen(MISFIT_LOG_FILE,"w");
/* Write header for misfit log file */ FPL2=fopen(MISFIT_LOG_FILE,"w");
if(GRAD_METHOD==1&&VERBOSE) {
if (TIME_FILT==0){ /* Write header for misfit log file */
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \n");} if(GRAD_METHOD==1&&VERBOSE) {
else{ if (TIME_FILT==0){
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t F_LOW_PASS \n"); fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \n");}
} else{
fprintf(FPL2,"opteps_vp \t epst1[1] \t epst1[2] \t epst1[3] \t L2t[1] \t L2t[2] \t L2t[3] \t L2t[4] \t F_LOW_PASS \n");
} }
} }
if(iter>1){FPL2=fopen(MISFIT_LOG_FILE,"a");}
if(WAVETYPE==3) FPL2_JOINT=fopen(L2_joint_log,"w");
} else {
FPL2=fopen(MISFIT_LOG_FILE,"a");
if(WAVETYPE==3) FPL2_JOINT=fopen(L2_joint_log,"a");
} }
} }
/* initialization of L2 calculation */ /* initialization of L2 calculation */
L2=0.0; L2=0.0;
Lcount=0;
energy=0.0; energy=0.0;
L2_all_shots=0.0; L2_all_shots=0.0;
energy_all_shots=0.0; energy_all_shots=0.0;
...@@ -1156,6 +1085,28 @@ int main(int argc, char **argv){ ...@@ -1156,6 +1085,28 @@ int main(int argc, char **argv){
/*------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------------*/
for (ishot=1;ishot<=nshots;ishot+=SHOTINC){ for (ishot=1;ishot<=nshots;ishot+=SHOTINC){
if (SEISMO && READREC==2){
if (ntr>0) {
dealloc_sections(ntr,ns,recpos_loc,sectionvx,sectionvy,sectionvz,sectionp,sectionpnp1,sectionpn,sectioncurl,sectiondiv,
sectionpdata,sectionpdiff,sectionpdiffold,sectionvxdata,sectionvxdiff,sectionvxdiffold,sectionvydata,
sectionvydiff,sectionvydiffold,sectionvzdata,sectionvzdiff,sectionvzdiffold);
}
free_imatrix(recpos,1,3,1,ntr_glob);
recpos=receiver(&ntr, srcpos, ishot);
recpos_loc = splitrec(recpos,&ntr_loc, ntr, recswitch);
ntr_glob=ntr;
ntr=ntr_loc;
if (ntr>0){
alloc_sections(ntr,ns,&sectionvx,&sectionvy,&sectionvz,&sectionp,&sectionpnp1,&sectionpn,&sectioncurl,&sectiondiv,
&sectionpdata,&sectionpdiff,&sectionpdiffold,&sectionvxdata,&sectionvxdiff,&sectionvxdiffold,&sectionvydata,
&sectionvydiff,&sectionvydiffold,&sectionvzdata,&sectionvzdiff,&sectionvzdiffold);
}
if (ntr) group_id = 1;
else group_id = 0;
MPI_Comm_split(MPI_COMM_WORLD, group_id, MYID, &MPI_COMM_NTR);
MPI_Comm_rank(MPI_COMM_NTR,