Commit 7da2cebb authored by Simone Butzer's avatar Simone Butzer

add additional input parameters

parent 67b9c79d
ifos3d
snapmerge
\ No newline at end of file
*.*
\ No newline at end of file
*.*
\ No newline at end of file
*.out
\ No newline at end of file
...@@ -10,9 +10,9 @@ ...@@ -10,9 +10,9 @@
#----------------------------------------------------------------------------------- #-----------------------------------------------------------------------------------
# #
#-------------- Domain Decomposition ----------------------------- #-------------- Domain Decomposition -----------------------------
number_of_processors_in_x-direction_(NPROCX) = 8 number_of_processors_in_x-direction_(NPROCX) = 2
number_of_processors_in_y-direction_(NPROCY) = 8 number_of_processors_in_y-direction_(NPROCY) = 2
number_of_processors_in_z-direction_(NPROCZ) = 8 number_of_processors_in_z-direction_(NPROCZ) = 2
# #
#-------------------- 3-D Grid ----------------------------------- #-------------------- 3-D Grid -----------------------------------
number_of_gridpoints_in_x-direction_(NX) = 160 number_of_gridpoints_in_x-direction_(NX) = 160
...@@ -161,6 +161,10 @@ minimum/maximum_iteration_number_(>0)_(ITMIN,ITMAX) = 1 , 80 ...@@ -161,6 +161,10 @@ minimum/maximum_iteration_number_(>0)_(ITMIN,ITMAX) = 1 , 80
filtering_(yes=1/no=0)_(FILT) = 1 filtering_(yes=1/no=0)_(FILT) = 1
maximum_number_frequencies_per_iteration_(NFMAX) = 5 maximum_number_frequencies_per_iteration_(NFMAX) = 5
number_of_timestep_per_period_used_for_inversion_(TAST) = 100 number_of_timestep_per_period_used_for_inversion_(TAST) = 100
average_model_parameter_(VP0,VS0,RHO0) = 6200.0, 3600.0, 2800.0
# velocities in m/s, density in kg/m³
parameter_class_weighting_factors_for_vp/vs/rho_(WEIGHT) = 1.0, 1.0, 0.0
# choose from 1.0 (full update) to 0.0 (no update)
# #
#------------------------Steplength estimation---------------------------------------- #------------------------Steplength estimation----------------------------------------
number_of_shots_used_for_steplength_estimation_(NSHOTS_STEP) = 4 number_of_shots_used_for_steplength_estimation_(NSHOTS_STEP) = 4
......
*.*
\ No newline at end of file
*.*
\ No newline at end of file
*.*
\ No newline at end of file
*.o
\ No newline at end of file
...@@ -8,7 +8,7 @@ ...@@ -8,7 +8,7 @@
# source code for model generation (elastic) # source code for model generation (elastic)
MODEL_SCR = hh_toy_florian.c MODEL_SCR = hh_toy.c
# Compiler (LAM: CC=hcc, CRAY T3E: CC=cc, SHARCNet: mpicc) # Compiler (LAM: CC=hcc, CRAY T3E: CC=cc, SHARCNet: mpicc)
......
...@@ -44,9 +44,9 @@ void exchange_par(void){ ...@@ -44,9 +44,9 @@ void exchange_par(void){
extern int METHOD; extern int METHOD;
extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS,LBFGS,EXTOBS; extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS,LBFGS,EXTOBS;
/*extern float F_INV;*/ /*extern float F_INV;*/
extern float TESTSTEP, WATER_HESS[3]; extern float TESTSTEP, WATER_HESS[3], WEIGHT[3], VP0, VS0, RHO0;
int idum[55]; int idum[55];
float fdum[36]; float fdum[42];
if (MYID == 0){ if (MYID == 0){
...@@ -81,10 +81,16 @@ void exchange_par(void){ ...@@ -81,10 +81,16 @@ void exchange_par(void){
fdum[29] = ALPHA; fdum[29] = ALPHA;
fdum[30] = BETA; fdum[30] = BETA;
fdum[31] = VPPML; fdum[31] = VPPML;
fdum[32] = TESTSTEP; fdum[32] = VP0;
fdum[33] = WATER_HESS[0]; fdum[33] = VS0;
fdum[34] = WATER_HESS[1]; fdum[34] = RHO0;
fdum[35] = WATER_HESS[2]; fdum[35] = WEIGHT[0];
fdum[36] = WEIGHT[1];
fdum[37] = WEIGHT[2];
fdum[38] = TESTSTEP;
fdum[39] = WATER_HESS[0];
fdum[40] = WATER_HESS[1];
fdum[41] = WATER_HESS[2];
idum[0] = FDORDER; idum[0] = FDORDER;
idum[1] = NPROCX; idum[1] = NPROCX;
...@@ -152,7 +158,7 @@ void exchange_par(void){ ...@@ -152,7 +158,7 @@ void exchange_par(void){
MPI_Barrier(MPI_COMM_WORLD); MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast(&idum,55,MPI_INT,0,MPI_COMM_WORLD); MPI_Bcast(&idum,55,MPI_INT,0,MPI_COMM_WORLD);
MPI_Bcast(&fdum,36,MPI_FLOAT,0,MPI_COMM_WORLD); MPI_Bcast(&fdum,42,MPI_FLOAT,0,MPI_COMM_WORLD);
MPI_Bcast(&SOURCE_FILE,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD); MPI_Bcast(&SOURCE_FILE,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD);
MPI_Bcast(&SIGNAL_FILE,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD); MPI_Bcast(&SIGNAL_FILE,STRING_SIZE,MPI_CHAR,0,MPI_COMM_WORLD);
...@@ -201,10 +207,16 @@ void exchange_par(void){ ...@@ -201,10 +207,16 @@ void exchange_par(void){
ALPHA = fdum[29]; ALPHA = fdum[29];
BETA = fdum[30]; BETA = fdum[30];
VPPML = fdum[31]; VPPML = fdum[31];
TESTSTEP=fdum[32]; VP0 = fdum[32];
WATER_HESS[0]=fdum[33]; VS0 = fdum[33];
WATER_HESS[1]= fdum[34]; RHO0 = fdum[34];
WATER_HESS[2]= fdum[35]; WEIGHT[0] = fdum[35];
WEIGHT[1] = fdum[36];
WEIGHT[2] = fdum[37];
TESTSTEP=fdum[38];
WATER_HESS[0]=fdum[39];
WATER_HESS[1]= fdum[40];
WATER_HESS[2]= fdum[41];
FDORDER = idum[0]; FDORDER = idum[0];
NPROCX = idum[1]; NPROCX = idum[1];
......
...@@ -57,4 +57,5 @@ float F_INV=0.0, TESTSTEP=0.0; ...@@ -57,4 +57,5 @@ float F_INV=0.0, TESTSTEP=0.0;
char MOD_OUT_FILE[STRING_SIZE], HESS_FILE[STRING_SIZE], GRAD_FILE[STRING_SIZE], SEIS_OBS_FILE[STRING_SIZE]; char MOD_OUT_FILE[STRING_SIZE], HESS_FILE[STRING_SIZE], GRAD_FILE[STRING_SIZE], SEIS_OBS_FILE[STRING_SIZE];
int DAMPTYPE, ITMIN, ITMAX, FILT, NFMAX, NSHOTS_STEP, TAST,EXTOBS; int DAMPTYPE, ITMIN, ITMAX, FILT, NFMAX, NSHOTS_STEP, TAST,EXTOBS;
int HESS, READ_HESS, REC_HESS, LBFGS; int HESS, READ_HESS, REC_HESS, LBFGS;
float WATER_HESS[3]={0.0, 0.0, 0.0}; float WATER_HESS[3]={0.0, 0.0, 0.0}, WEIGHT[3]={0.0, 0.0, 0.0};
\ No newline at end of file float VP0, VS0, RHO0;
\ No newline at end of file
...@@ -27,14 +27,14 @@ ...@@ -27,14 +27,14 @@
void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***grad1, float ***grad2,float ***grad3, float ***hess1, float ***hess2, float ***hess3, float finv, int iteration){ void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***grad1, float ***grad2,float ***grad3, float ***hess1, float ***hess2, float ***hess3, float finv, int iteration){
extern float WATER_HESS[3]; extern float WATER_HESS[3], VP0, VS0, RHO0;
extern int NXG, NYG, NZG, FW, NX, NY, NZ; extern int NXG, NYG, NZG, FW, NX, NY, NZ;
extern FILE *FP; extern FILE *FP;
int LBFGS=0; int LBFGS=0;
float wl[3],buf[3]; float wl[3],buf[3];
int i,j,k,buf1; int i,j,k,buf1;
float vp0=1240.0, vs0=730.0, rho0=1820.0;
buf[0]=0.0;buf[1]=0.0;buf[2]=0.0; buf[0]=0.0;buf[1]=0.0;buf[2]=0.0;
wl[0]=0.0;wl[1]=0.0;wl[2]=0.0; wl[0]=0.0;wl[1]=0.0;wl[2]=0.0;
...@@ -44,9 +44,9 @@ void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***g ...@@ -44,9 +44,9 @@ void hess_apply(int nx1, int nx2, int ny1, int ny2, int nz1, int nz2, float ***g
for (j=1;j<=ny2;j++){ for (j=1;j<=ny2;j++){
for (i=1;i<=nx2;i++){ for (i=1;i<=nx2;i++){
for (k=1;k<=nz2;k++){ for (k=1;k<=nz2;k++){
grad1[j][i][k]=grad1[j][i][k]*vp0; grad1[j][i][k]=grad1[j][i][k]*VP0;
grad2[j][i][k]=grad2[j][i][k]*vs0; grad2[j][i][k]=grad2[j][i][k]*VS0;
grad3[j][i][k]=grad3[j][i][k]*rho0; grad3[j][i][k]=grad3[j][i][k]*RHO0;
/*hess1[j][i][k]=hess1[j][i][k]*pow(vp0,2.0); /*hess1[j][i][k]=hess1[j][i][k]*pow(vp0,2.0);
hess2[j][i][k]=hess2[j][i][k]*pow(vs0,2.0); hess2[j][i][k]=hess2[j][i][k]*pow(vs0,2.0);
hess3[j][i][k]=hess3[j][i][k]*pow(rho0,2.0);*/ hess3[j][i][k]=hess3[j][i][k]*pow(rho0,2.0);*/
......
...@@ -27,20 +27,24 @@ void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float ...@@ -27,20 +27,24 @@ void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float
extern FILE *FP; extern FILE *FP;
extern int NX, NY, NZ, LBFGS; extern int NX, NY, NZ, LBFGS;
extern float VP0, VS0, RHO0, WEIGHT[3];
int j,i,k,l,l1; int j,i,k,l,l1;
float vpnew,vsnew,rhonew; float vpnew,vsnew,rhonew;
float max[3],buf[3],max1[3],dummy[3]; float max[3],buf[3],max1[3],dummy[3];
float vp,vs; float vp,vs;
int bfgsnum=5,w; int bfgsnum=5,w;
float scale1=1.0, scale2=1.0; float scale1=0.0, scale2=0.0, scale3=0.0;
float vp0=6200.0, vs0=3600.0, rho0=2800.0;
buf[0]=0.0;buf[1]=0.0;buf[2]=0.0; buf[0]=0.0;buf[1]=0.0;buf[2]=0.0;
max[0]=0.0;max[1]=0.0;max[2]=0.0; max[0]=0.0;max[1]=0.0;max[2]=0.0;
dummy[0]=0.0; dummy[1]=0.0; dummy[2]=0.0; dummy[0]=0.0; dummy[1]=0.0; dummy[2]=0.0;
max1[0]=0.0;max1[1]=0.0;max1[2]=0.0; max1[0]=0.0;max1[1]=0.0;max1[2]=0.0;
scale1=WEIGHT[0];
scale2=WEIGHT[1];
scale3=WEIGHT[2];
/*if(iteration<20)scale1=0.7; /*if(iteration<20)scale1=0.7;
if(iteration>40)scale2=0.7; if(iteration>40)scale2=0.7;
if(iteration>60)scale2=0.4;*/ if(iteration>60)scale2=0.4;*/
...@@ -126,20 +130,20 @@ void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float ...@@ -126,20 +130,20 @@ void modelupdate(int nx, int ny, int nz, float ***gradvp, float ***gradvs, float
/*update model*/ /*update model*/
vpnew=0.0; vpnew=0.0;
/*vpnew=sqrt(pi[j][i][k]/rho[j][i][k])+max[0]*step*gradconvp[j][i][k];*/ /*vpnew=sqrt(pi[j][i][k]/rho[j][i][k])+max[0]*step*gradconvp[j][i][k];*/
vpnew=sqrt(pi[j][i][k]/rho[j][i][k])+step*scale1*gradvp[j][i][k]*vp0; vpnew=sqrt(pi[j][i][k]/rho[j][i][k])+step*scale1*gradvp[j][i][k]*VP0;
vsnew=0.0; vsnew=0.0;
/*vsnew=sqrt(u[j][i][k]/rho[j][i][k])+max[1]*step*gradconvs[j][i][k];*/ /*vsnew=sqrt(u[j][i][k]/rho[j][i][k])+max[1]*step*gradconvs[j][i][k];*/
vsnew=sqrt(u[j][i][k]/rho[j][i][k])+step*scale2*gradvs[j][i][k]*vs0; vsnew=sqrt(u[j][i][k]/rho[j][i][k])+step*scale2*gradvs[j][i][k]*VS0;
rhonew=0.0; rhonew=0.0;
rhonew=rho[j][i][k]+0.0*step*gradrho[j][i][k]*rho0; rhonew=rho[j][i][k]+scale3*step*gradrho[j][i][k]*RHO0;
if(LBFGS){ if(LBFGS){
/*save normalised model differences for LBFGS*/ /*save normalised model differences for LBFGS*/
l++; l1++; l++; l1++;
bfgsmod1[w][l]=(vpnew-sqrt(pi[j][i][k]/rho[j][i][k]))/vp0; bfgsmod1[w][l]=(vpnew-sqrt(pi[j][i][k]/rho[j][i][k]))/VP0;
bfgsmod1[w][l1]=(vsnew-sqrt(u[j][i][k]/rho[j][i][k]))/vs0; bfgsmod1[w][l1]=(vsnew-sqrt(u[j][i][k]/rho[j][i][k]))/VS0;
/*bfgsmod3[w][l]=(rhonew-rho[j][i][k])/rho0;*/ /*bfgsmod3[w][l]=(rhonew-rho[j][i][k])/rho0;*/
} }
......
...@@ -48,7 +48,7 @@ int read_par(FILE *fp_in){ ...@@ -48,7 +48,7 @@ int read_par(FILE *fp_in){
extern int METHOD; extern int METHOD;
extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS,EXTOBS,LBFGS; extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS,EXTOBS,LBFGS;
/*extern float F_INV;*/ /*extern float F_INV;*/
extern float TESTSTEP, WATER_HESS[3]; extern float TESTSTEP, WATER_HESS[3], WEIGHT[3], VP0, VS0, RHO0;
/* definition of local variables */ /* definition of local variables */
char s[256], cline[256]=""; char s[256], cline[256]="";
...@@ -389,8 +389,7 @@ int read_par(FILE *fp_in){ ...@@ -389,8 +389,7 @@ int read_par(FILE *fp_in){
case 69 : case 69 :
sscanf(cline,"%s =%i , %i",s,&ITMIN,&ITMAX); sscanf(cline,"%s =%i , %i",s,&ITMIN,&ITMAX);
if(!METHOD) ITMAX=1; if(!METHOD) ITMAX=1;
break; break;
case 70 : case 70 :
sscanf(cline,"%s =%i",s,&FILT); sscanf(cline,"%s =%i",s,&FILT);
break; break;
...@@ -400,31 +399,38 @@ int read_par(FILE *fp_in){ ...@@ -400,31 +399,38 @@ int read_par(FILE *fp_in){
case 72 : case 72 :
sscanf(cline,"%s =%i",s,&TAST); sscanf(cline,"%s =%i",s,&TAST);
break; break;
case 73 : case 73:
sscanf(cline,"%s =%i",s,&NSHOTS_STEP); sscanf(cline,"%s =%f, %f, %f",s,&VP0,&VS0,&RHO0);
break; break;
case 74 : case 74 :
sscanf(cline,"%s =%f, %f, %f",s,&WEIGHT[0],&WEIGHT[1],&WEIGHT[2]);
break;
case 75 :
sscanf(cline,"%s =%i",s,&NSHOTS_STEP);
break;
case 76 :
sscanf(cline,"%s =%f",s,&TESTSTEP); sscanf(cline,"%s =%f",s,&TESTSTEP);
break; break;
case 75 : case 77 :
sscanf(cline,"%s =%i",s,&DAMPTYPE); sscanf(cline,"%s =%i",s,&DAMPTYPE);
break; break;
case 76 : case 78 :
sscanf(cline,"%s =%i",s,&HESS); sscanf(cline,"%s =%i",s,&HESS);
break; break;
case 77 : case 79 :
sscanf(cline,"%s =%i",s,&READ_HESS); sscanf(cline,"%s =%i",s,&READ_HESS);
break; break;
case 78 : case 80 :
sscanf(cline,"%s =%i",s,&REC_HESS); sscanf(cline,"%s =%i",s,&REC_HESS);
break; break;
case 79 : case 81 :
sscanf(cline,"%s =%f, %f, %f",s,&WATER_HESS[0],&WATER_HESS[1],&WATER_HESS[2]); sscanf(cline,"%s =%f, %f, %f",s,&WATER_HESS[0],&WATER_HESS[1],&WATER_HESS[2]);
break; break;
case 80 : case 82 :
sscanf(cline,"%s =%i",s,&LBFGS); sscanf(cline,"%s =%i",s,&LBFGS);
break; break;
case 81 : case 83 :
nvarin=sscanf(cline,"%s =%i ,%i ,%i",s,&ASCIIEBCDIC,&LITTLEBIG,&IEEEIBM); nvarin=sscanf(cline,"%s =%i ,%i ,%i",s,&ASCIIEBCDIC,&LITTLEBIG,&IEEEIBM);
switch(nvarin){ switch(nvarin){
case 0: ; case 0: ;
...@@ -439,7 +445,7 @@ int read_par(FILE *fp_in){ ...@@ -439,7 +445,7 @@ int read_par(FILE *fp_in){
} }
} }
LOG=0; LOG=0;
if(lineno<81) fprintf(stderr," Warning: only %d non-commentary lines of input parameters read (expected 81).\n",lineno); if(lineno<83) fprintf(stderr," Warning: only %d non-commentary lines of input parameters read (expected 83).\n",lineno);
/* else if (lineno>71) fprintf(stderr," Warning: %d non-commentary lines of input parameters read \n \t(expected and interpreted: 67).\n",lineno); */ /* else if (lineno>71) fprintf(stderr," Warning: %d non-commentary lines of input parameters read \n \t(expected and interpreted: 67).\n",lineno); */
/* else fprintf(stderr," %d non-commentary lines of input parameters read.\n",lineno); */ /* else fprintf(stderr," %d non-commentary lines of input parameters read.\n",lineno); */
fclose(fp_in); fclose(fp_in);
......
...@@ -26,8 +26,8 @@ ...@@ -26,8 +26,8 @@
void readhess(int nx, int ny, int nz, float *** hess1, float *** hess2, float ***hess3){ void readhess(int nx, int ny, int nz, float *** hess1, float *** hess2, float ***hess3){
extern int NX, NY, NZ, NXG, NYG, NZG, POS[4]; extern int NX, NY, NZ, NXG, NYG, NZG, POS[4];
extern float VP0, VS0, RHO0;
extern FILE *FP; extern FILE *FP;
float vp0=1240.0, vs0=730.0, rho0=1820.0;
/* local variables */ /* local variables */
float rhov, vp, vs; /*qp, qs;*/ float rhov, vp, vs; /*qp, qs;*/
...@@ -78,9 +78,9 @@ float vp0=1240.0, vs0=730.0, rho0=1820.0; ...@@ -78,9 +78,9 @@ float vp0=1240.0, vs0=730.0, rho0=1820.0;
hess3[jj][ii][kk]=rhov;*/ hess3[jj][ii][kk]=rhov;*/
hess1[jj][ii][kk]=vp*pow(vp0,2); hess1[jj][ii][kk]=vp*pow(VP0,2);
hess2[jj][ii][kk]=vs*pow(vs0,2); hess2[jj][ii][kk]=vs*pow(VS0,2);
hess3[jj][ii][kk]=rhov*pow(rho0,2); hess3[jj][ii][kk]=rhov*pow(RHO0,2);
} }
......
...@@ -42,7 +42,7 @@ void writepar(FILE *fp, int ns){ ...@@ -42,7 +42,7 @@ void writepar(FILE *fp, int ns){
extern int METHOD; extern int METHOD;
extern int NP, NPROCX, NPROCY, NPROCZ, MYID; extern int NP, NPROCX, NPROCY, NPROCZ, MYID;
extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS, LBFGS,EXTOBS; extern int ITMIN, ITMAX, FILT, NFMAX, TAST, NSHOTS_STEP, DAMPTYPE, HESS, READ_HESS, REC_HESS, LBFGS,EXTOBS;
extern float TESTSTEP,WATER_HESS[3]; extern float TESTSTEP,WATER_HESS[3], WEIGHT[3], VP0, VS0, RHO0;
/* definition of local variables */ /* definition of local variables */
char th1[3], file_ext[8]; char th1[3], file_ext[8];
char th2[3]; char th2[3];
...@@ -431,6 +431,8 @@ fprintf(fp," \n minimum/maximum_iteration_number: %d,%d\n",ITMIN,ITMAX); ...@@ -431,6 +431,8 @@ fprintf(fp," \n minimum/maximum_iteration_number: %d,%d\n",ITMIN,ITMAX);
fprintf(fp," \n filtering: %d\n",FILT); fprintf(fp," \n filtering: %d\n",FILT);
fprintf(fp," \n maximum_number_frequencies_per_iteration: %d\n",NFMAX); fprintf(fp," \n maximum_number_frequencies_per_iteration: %d\n",NFMAX);
fprintf(fp," \n number_of_timestep_per_wavelength_used_for_inversion: %d\n",TAST); fprintf(fp," \n number_of_timestep_per_wavelength_used_for_inversion: %d\n",TAST);
fprintf(fp," \n average_model_parameter VP0=%5.2f m/s, VS0=%5.2f m/s, RHO0=%5.2f kg/m^3\n",VP0, VS0, RHO0);
fprintf(fp," \n parameter_class_weighting_factors_for_vp: %5.2f, vs: %5.2f, rho: %5.2f\n",WEIGHT[0], WEIGHT[1], WEIGHT[2]);
fprintf(fp," \n------------------------Steplength estimation----------------------------------------\n"); fprintf(fp," \n------------------------Steplength estimation----------------------------------------\n");
fprintf(fp,"\n number_of_shots_used_for_steplength_estimation: %d\n",NSHOTS_STEP); fprintf(fp,"\n number_of_shots_used_for_steplength_estimation: %d\n",NSHOTS_STEP);
......
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