Commit 09fcd4b8 authored by Florian Wittkamp's avatar Florian Wittkamp

Merge branch 'Joint_development' into development

parents 2f88f362 47f69677
......@@ -2337,7 +2337,6 @@ int main(int argc, char **argv){
hin++;
}
/*fclose(FP2);*/
if((EPRECOND==1)&&(EPRECOND_ITER==iter||(EPRECOND_ITER==0))){
if (WAVETYPE==1 || WAVETYPE==3) eprecond(Wr,pvx,pvy);
......@@ -2551,7 +2550,6 @@ int main(int argc, char **argv){
if(INVMAT1==1){
/* calculate density gradient */
waveconv_rho_shot_z[j][i] = ( (pu[j][i] * pu[j][i] * waveconv_mu_z[j][i]) + waveconv_rho_s_z[j][i]);
}
......@@ -2570,15 +2568,18 @@ int main(int argc, char **argv){
/* calculate and apply energy preconditioning */
/* -------------------------------------------- */
if((EPRECOND==1)||(EPRECOND==3)){
/* calculate energy weights */
if(EPRECOND_ITER==iter||(EPRECOND_ITER==0)) {
fprintf(FP,"\n Calculating approx. Hessian for shot %i. EPRECOND=%i, EPSILON_WE=%f",ishot,EPRECOND,EPSILON_WE);
if(WAVETYPE==1 || WAVETYPE==3) {
eprecond1(We,Ws,Wr);
eprecond1(We,Ws,Wr,EPSILON_WE);
if(EPRECOND_PER_SHOT) We_max=global_maximum(We);
}
if(WAVETYPE==2 || WAVETYPE==3) {
eprecond1(We_SH,Ws_SH,Wr_SH);
eprecond1(We_SH,Ws_SH,Wr_SH,EPSILON_WE_SH);
if(EPRECOND_PER_SHOT) We_max_SH=global_maximum(We_SH);
}
......
/* This function is copied from DENISE Black Edition from D. Koehn */
#include "fd.h"
void eprecond1(float ** We, float ** Ws, float ** Wr){
void eprecond1(float ** We, float ** Ws, float ** Wr, float epsilon){
extern int NX, NY, IDX, IDY, DTINV, EPRECOND, VERBOSE;
extern int POS[3], NXG;
extern float DH;
int i, j, k, l, ii, jj;
float maxWetmp, maxWe, x, y, xmin, xmax;
extern float EPSILON_WE;
xmin = DH;
xmax = NXG*DH;
......@@ -56,11 +55,12 @@ void eprecond1(float ** We, float ** Ws, float ** Wr){
/* estimate maximum of We */
MPI_Allreduce(&maxWetmp,&maxWe,1,MPI_FLOAT,MPI_MAX,MPI_COMM_WORLD);
/* regularize energy weighting to avoid divison by zero */
for (j=1;j<=NY;j=j+IDY){
for (i=1;i<=NX;i=i+IDX){
We[j][i] = We[j][i] + (EPSILON_WE*maxWe);
We[j][i] = We[j][i] + (epsilon*maxWe);
}
}
......
......@@ -99,7 +99,7 @@ void exchange_par(void){
extern int EPRECOND;
extern int EPRECOND_ITER;
extern float EPSILON_WE;
extern float EPSILON_WE,EPSILON_WE_SH;
extern int EPRECOND_PER_SHOT;
extern int LBFGS_SURFACE;
......@@ -113,16 +113,18 @@ void exchange_par(void){
extern float WOLFE_C1_SL;
extern float WOLFE_C2_SL;
/* definition of local variables */
/* definition of local variables */
/* NPAR is set in fd.h and must be set to the max. number of elements in idum/fdum +1 (because vector index starts at 0) */
int idum[NPAR];
float fdum[NPAR];
if (MYID == 0){
/***********/
/* Float */
/***********/
fdum[1] = DH;
fdum[2] = TIME;
fdum[3] = DT;
......@@ -201,10 +203,14 @@ void exchange_par(void){
fdum[61] = JOINT_INVERSION_PSV_SH_ALPHA_RHO;
fdum[62]=EPSILON_WE;
fdum[63]=EPSILON_WE_SH;
fdum[63]=WOLFE_C1_SL;
fdum[64]=WOLFE_C2_SL;
fdum[64]=WOLFE_C1_SL;
fdum[65]=WOLFE_C2_SL;
/***********/
/* Integer */
/***********/
idum[1] = NPROCX;
idum[2] = NPROCY;
......@@ -391,6 +397,9 @@ void exchange_par(void){
MPI_Barrier(MPI_COMM_WORLD);
/***********/
/* Float */
/***********/
DH=fdum[1];
TIME=fdum[2];
DT=fdum[3];
......@@ -468,9 +477,14 @@ void exchange_par(void){
JOINT_INVERSION_PSV_SH_ALPHA_RHO = fdum[61];
EPSILON_WE=fdum[62];
EPSILON_WE_SH=fdum[63];
WOLFE_C1_SL=fdum[64];
WOLFE_C2_SL=fdum[65];
WOLFE_C1_SL=fdum[63];
WOLFE_C2_SL=fdum[64];
/***********/
/* Integer */
/***********/
NPROCX = idum[1];
NPROCY = idum[2];
......
......@@ -520,7 +520,7 @@ void apply_workflow(float ** workflow,int workflow_lines,char workflow_header[ST
void eprecond(float ** W, float ** vx, float ** vy);
void eprecond_SH(float ** W, float ** vz);
void eprecond1(float ** We, float ** Ws, float ** Wr);
void eprecond1(float ** We, float ** Ws, float ** Wr, float epsilon);
/* Matrix Operations */
float average_matrix(float ** matrix);
......
......@@ -85,7 +85,7 @@ float EPS_SCALE, SCALEFAC;
float PRO;
int EPRECOND;
float EPSILON_WE;
float EPSILON_WE,EPSILON_WE_SH;
int EPRECOND_ITER;
int EPRECOND_PER_SHOT;
......
......@@ -25,9 +25,6 @@
* ----------------------------------------------------------------------*/
#include "fd.h"
/* Declare local functions */
float ** joint_inversion_grad ( float ** gradiant_1,float ** gradiant_2, float alpha, int joint_type) {
/*
* Input parameters:
......@@ -57,11 +54,31 @@ float ** joint_inversion_grad ( float ** gradiant_1,float ** gradiant_2, float a
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
gradiant_1[j][i]=((1-alpha)*gradiant_1[j][i]/max1+alpha*gradiant_2[j][i]/max2);
if (isnan(gradiant_1[j][i])) {
gradiant_1[j][i]=0.0;
}
}
}
joint_gradiant=gradiant_1;
break;
case 2:
/* joint_type=2: Arithmetic mean without normalization */
for (j=1;j<=NY;j++){
for (i=1;i<=NX;i++){
gradiant_1[j][i]=((1-alpha)*gradiant_1[j][i]+alpha*gradiant_2[j][i]);
if (isnan(gradiant_1[j][i])) {
gradiant_1[j][i]=0.0;
}
}
}
joint_gradiant=gradiant_1;
break;
}
return joint_gradiant;
};
......
......@@ -106,7 +106,7 @@ void read_par_json(FILE *fp, char *fileinp){
extern int EPRECOND;
extern int EPRECOND_ITER;
extern float EPSILON_WE;
extern float EPSILON_WE,EPSILON_WE_SH;
extern int EPRECOND_PER_SHOT;
extern int LBFGS_SURFACE;
......@@ -491,6 +491,12 @@ void read_par_json(FILE *fp, char *fileinp){
}
if (get_float_from_objectlist("EPSILON_WE",number_readobjects,&EPSILON_WE,varname_list, value_list))
err("Variable EPSILON_WE could not be retrieved from the json input file!");
if (get_float_from_objectlist("EPSILON_WE_SH",number_readobjects,&EPSILON_WE_SH,varname_list, value_list)) {
EPSILON_WE_SH=EPSILON_WE;
fprintf(fp,"Variable EPSILON_WE_SH is set to EPSILON_WE=%f.\n",EPSILON_WE_SH);
}
}
if (get_int_from_objectlist("TESTSHOT_START",number_readobjects,&TESTSHOT_START,varname_list, value_list))
......
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