/*----------------------------------------------------------------------------------------- * 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 . -----------------------------------------------------------------------------------------*/ /*------------------------------------------------------------------------ * updating particle velocities at gridpoints [nx1...nx2][ny1...ny2] * by a staggered grid finite difference scheme of FDORDER accuracy in space * and second order accuracy in time * * ----------------------------------------------------------------------*/ #include "fd.h" void update_v_PML(int nx1, int nx2, int ny1, int ny2, int nt, float ** vx, float ** vxp1, float ** vxm1, float ** vy, float ** vyp1, float ** vym1, float ** uttx, float ** utty,float ** sxx, float ** syy, float ** sxy, float **rip, float **rjp, float ** srcpos_loc, float ** signals, float ** signals1, int nsrc, float ** absorb_coeff, float *hc, int infoout,int sw, float * K_x, float * a_x, float * b_x, float * K_x_half, float * a_x_half, float * b_x_half, float * K_y, float * a_y, float * b_y, float * K_y_half, float * a_y_half, float * b_y_half, float ** psi_sxx_x, float ** psi_syy_y, float ** psi_sxy_y, float ** psi_sxy_x){ int i, j,l,fdoh,m, h, h1; float dtdh, azi_rad; float vxtmp, vytmp; float sxx_x, syy_y, sxy_y, sxy_x; extern float DT, DH; double time1 = 0.0, time2; extern int MYID, SOURCE_TYPE, ADJOINT_TYPE, FDORDER; extern int FDORDER, PARAMETERIZATION; extern int FREE_SURF, BOUNDARY, FW; extern int NPROCX, NPROCY, POS[3]; extern FILE *FP; extern int VELOCITY; fdoh = FDORDER/2; dtdh = DT*DT/DH; /* drad = PI/180.0; angle = 135.0; */ if (infoout && (MYID==0)){ time1=MPI_Wtime(); fprintf(FP,"\n **Message from update_v (printed by PE %d):\n",MYID); fprintf(FP," Updating particle velocities ..."); } /* ------------------------------------------------------------ * Important! * rip and rjp are reciprocal values of averaged densities * ------------------------------------------------------------ */ switch (FDORDER){ case 2: for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ sxx_x = hc[1]*(sxx[j][i+1]-sxx[j][i]); sxy_x = hc[1]*(sxy[j][i]-sxy[j][i-1]); sxy_y = hc[1]*(sxy[j][i]-sxy[j-1][i]); syy_y = hc[1]*(syy[j+1][i]-syy[j][i]); /* left boundary */ if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){ psi_sxx_x[j][i] = b_x_half[i] * psi_sxx_x[j][i] + a_x_half[i] * sxx_x; sxx_x = sxx_x / K_x_half[i] + psi_sxx_x[j][i]; psi_sxy_x[j][i] = b_x[i] * psi_sxy_x[j][i] + a_x[i] * sxy_x; sxy_x = sxy_x / K_x[i] + psi_sxy_x[j][i]; } /* right boundary */ if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){ h1 = (i-nx2+2*FW); h=i; psi_sxx_x[j][h1] = b_x_half[h1] * psi_sxx_x[j][h1] + a_x_half[h1] * sxx_x; sxx_x = sxx_x / K_x_half[h1] + psi_sxx_x[j][h1]; psi_sxy_x[j][h1] = b_x[h1] * psi_sxy_x[j][h1] + a_x[h1] * sxy_x; sxy_x = sxy_x / K_x[h1] + psi_sxy_x[j][h1]; } /* top boundary */ if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){ psi_syy_y[j][i] = b_y_half[j] * psi_syy_y[j][i] + a_y_half[j] * syy_y; syy_y = syy_y / K_y_half[j] + psi_syy_y[j][i]; psi_sxy_y[j][i] = b_y[j] * psi_sxy_y[j][i] + a_y[j] * sxy_y; sxy_y = sxy_y / K_y[j] + psi_sxy_y[j][i]; } /* bottom boundary */ if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){ h1 = (j-ny2+2*FW); h = j; psi_syy_y[h1][i] = b_y_half[h1] * psi_syy_y[h1][i] + a_y_half[h1] * syy_y; syy_y = syy_y / K_y_half[h1] + psi_syy_y[h1][i]; psi_sxy_y[h1][i] = b_y[h1] * psi_sxy_y[h1][i] + a_y[h1] * sxy_y; sxy_y = sxy_y / K_y[h1] + psi_sxy_y[h1][i]; } if(sw==0){ if (VELOCITY==0){ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH; } else{ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH; } } if(sw==1){ if (VELOCITY==0){ vxp1[j][i] += vx[j][i]*DT; vyp1[j][i] += vy[j][i]*DT; } else{ vxp1[j][i] = vx[j][i]; vyp1[j][i] = vy[j][i]; } } vx[j][i] += DT*rip[j][i]*(sxx_x+sxy_y)/DH; vy[j][i] += DT*rjp[j][i]*(sxy_x+syy_y)/DH; } } break; case 4: for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ sxx_x = hc[1]*(sxx[j][i+1]-sxx[j][i]) + hc[2]*(sxx[j][i+2]-sxx[j][i-1]); sxy_x = hc[1]*(sxy[j][i]-sxy[j][i-1]) + hc[2]*(sxy[j][i+1]-sxy[j][i-2]); sxy_y = hc[1]*(sxy[j][i]-sxy[j-1][i]) + hc[2]*(sxy[j+1][i]-sxy[j-2][i]); syy_y = hc[1]*(syy[j+1][i]-syy[j][i]) + hc[2]*(syy[j+2][i]-syy[j-1][i]); /* left boundary */ if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){ psi_sxx_x[j][i] = b_x_half[i] * psi_sxx_x[j][i] + a_x_half[i] * sxx_x; sxx_x = sxx_x / K_x_half[i] + psi_sxx_x[j][i]; psi_sxy_x[j][i] = b_x[i] * psi_sxy_x[j][i] + a_x[i] * sxy_x; sxy_x = sxy_x / K_x[i] + psi_sxy_x[j][i]; } /* right boundary */ if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){ h1 = (i-nx2+2*FW); h=i; /*psi_sxx_x[j][i] = b_x_half[i] * psi_sxx_x[j][i] + a_x_half[i] * sxx_x; sxx_x = sxx_x / K_x_half[i] + psi_sxx_x[j][i];*/ psi_sxx_x[j][h1] = b_x_half[h1] * psi_sxx_x[j][h1] + a_x_half[h1] * sxx_x; sxx_x = sxx_x / K_x_half[h1] + psi_sxx_x[j][h1]; psi_sxy_x[j][h1] = b_x[h1] * psi_sxy_x[j][h1] + a_x[h1] * sxy_x; sxy_x = sxy_x / K_x[h1] + psi_sxy_x[j][h1]; } /* top boundary */ if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){ psi_syy_y[j][i] = b_y_half[j] * psi_syy_y[j][i] + a_y_half[j] * syy_y; syy_y = syy_y / K_y_half[j] + psi_syy_y[j][i]; psi_sxy_y[j][i] = b_y[j] * psi_sxy_y[j][i] + a_y[j] * sxy_y; sxy_y = sxy_y / K_y[j] + psi_sxy_y[j][i]; } /* bottom boundary */ if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){ h1 = (j-ny2+2*FW); h = j; /*psi_syy_y[j][i] = b_y_half[j] * psi_syy_y[j][i] + a_y_half[j] * syy_y; syy_y = syy_y / K_y_half[j] + psi_syy_y[j][i];*/ psi_syy_y[h1][i] = b_y_half[h1] * psi_syy_y[h1][i] + a_y_half[h1] * syy_y; syy_y = syy_y / K_y_half[h1] + psi_syy_y[h1][i]; psi_sxy_y[h1][i] = b_y[h1] * psi_sxy_y[h1][i] + a_y[h1] * sxy_y; sxy_y = sxy_y / K_y[h1] + psi_sxy_y[h1][i]; } if(sw==0){ if (VELOCITY==0){ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH;} else{ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH;} } if(sw==1){ if (VELOCITY==0){ vxp1[j][i] += vx[j][i]*DT; vyp1[j][i] += vy[j][i]*DT;} else{ vxp1[j][i] = vx[j][i]; vyp1[j][i] = vy[j][i];} } vx[j][i] += DT*rip[j][i]*(sxx_x+sxy_y)/DH; vy[j][i] += DT*rjp[j][i]*(sxy_x+syy_y)/DH; } } break; case 6: for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ sxx_x = hc[1]*(sxx[j][i+1]-sxx[j][i]) + hc[2]*(sxx[j][i+2]-sxx[j][i-1]) + hc[3]*(sxx[j][i+3]-sxx[j][i-2]); sxy_x = hc[1]*(sxy[j][i]-sxy[j][i-1]) + hc[2]*(sxy[j][i+1]-sxy[j][i-2]) + hc[3]*(sxy[j][i+2]-sxy[j][i-3]); sxy_y = hc[1]*(sxy[j][i]-sxy[j-1][i]) + hc[2]*(sxy[j+1][i]-sxy[j-2][i]) + hc[3]*(sxy[j+2][i]-sxy[j-3][i]); syy_y = hc[1]*(syy[j+1][i]-syy[j][i]) + hc[2]*(syy[j+2][i]-syy[j-1][i]) + hc[3]*(syy[j+3][i]-syy[j-2][i]); /* left boundary */ if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){ psi_sxx_x[j][i] = b_x_half[i] * psi_sxx_x[j][i] + a_x_half[i] * sxx_x; sxx_x = sxx_x / K_x_half[i] + psi_sxx_x[j][i]; psi_sxy_x[j][i] = b_x[i] * psi_sxy_x[j][i] + a_x[i] * sxy_x; sxy_x = sxy_x / K_x[i] + psi_sxy_x[j][i]; } /* right boundary */ if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){ h1 = (i-nx2+2*FW); h=i; /*psi_sxx_x[j][i] = b_x_half[i] * psi_sxx_x[j][i] + a_x_half[i] * sxx_x; sxx_x = sxx_x / K_x_half[i] + psi_sxx_x[j][i];*/ psi_sxx_x[j][h1] = b_x_half[h1] * psi_sxx_x[j][h1] + a_x_half[h1] * sxx_x; sxx_x = sxx_x / K_x_half[h1] + psi_sxx_x[j][h1]; psi_sxy_x[j][h1] = b_x[h1] * psi_sxy_x[j][h1] + a_x[h1] * sxy_x; sxy_x = sxy_x / K_x[h1] + psi_sxy_x[j][h1]; } /* top boundary */ if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){ psi_syy_y[j][i] = b_y_half[j] * psi_syy_y[j][i] + a_y_half[j] * syy_y; syy_y = syy_y / K_y_half[j] + psi_syy_y[j][i]; psi_sxy_y[j][i] = b_y[j] * psi_sxy_y[j][i] + a_y[j] * sxy_y; sxy_y = sxy_y / K_y[j] + psi_sxy_y[j][i]; } /* bottom boundary */ if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){ h1 = (j-ny2+2*FW); h = j; /*psi_syy_y[j][i] = b_y_half[j] * psi_syy_y[j][i] + a_y_half[j] * syy_y; syy_y = syy_y / K_y_half[j] + psi_syy_y[j][i];*/ psi_syy_y[h1][i] = b_y_half[h1] * psi_syy_y[h1][i] + a_y_half[h1] * syy_y; syy_y = syy_y / K_y_half[h1] + psi_syy_y[h1][i]; psi_sxy_y[h1][i] = b_y[h1] * psi_sxy_y[h1][i] + a_y[h1] * sxy_y; sxy_y = sxy_y / K_y[h1] + psi_sxy_y[h1][i]; } if(sw==0){ if (VELOCITY==0){ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH;} else{ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH;} } if(sw==1){ if (VELOCITY==0){ vxp1[j][i] += vx[j][i]*DT; vyp1[j][i] += vy[j][i]*DT;} else{ vxp1[j][i] = vx[j][i]; vyp1[j][i] = vy[j][i];} } vx[j][i] += DT*rip[j][i]*(sxx_x+sxy_y)/DH; vy[j][i] += DT*rjp[j][i]*(sxy_x+syy_y)/DH; } } break; case 8: for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ sxx_x = hc[1]*(sxx[j][i+1]-sxx[j][i]) + hc[2]*(sxx[j][i+2]-sxx[j][i-1]) + hc[3]*(sxx[j][i+3]-sxx[j][i-2]) + hc[4]*(sxx[j][i+4]-sxx[j][i-3]); sxy_x = hc[1]*(sxy[j][i]-sxy[j][i-1]) + hc[2]*(sxy[j][i+1]-sxy[j][i-2]) + hc[3]*(sxy[j][i+2]-sxy[j][i-3]) + hc[4]*(sxy[j][i+3]-sxy[j][i-4]); sxy_y = hc[1]*(sxy[j][i]-sxy[j-1][i]) + hc[2]*(sxy[j+1][i]-sxy[j-2][i]) + hc[3]*(sxy[j+2][i]-sxy[j-3][i]) + hc[4]*(sxy[j+3][i]-sxy[j-4][i]); syy_y = hc[1]*(syy[j+1][i]-syy[j][i]) + hc[2]*(syy[j+2][i]-syy[j-1][i]) + hc[3]*(syy[j+3][i]-syy[j-2][i]) + hc[4]*(syy[j+4][i]-syy[j-3][i]); /* left boundary */ if((!BOUNDARY) && (POS[1]==0) && (i<=FW)){ psi_sxx_x[j][i] = b_x_half[i] * psi_sxx_x[j][i] + a_x_half[i] * sxx_x; sxx_x = sxx_x / K_x_half[i] + psi_sxx_x[j][i]; psi_sxy_x[j][i] = b_x[i] * psi_sxy_x[j][i] + a_x[i] * sxy_x; sxy_x = sxy_x / K_x[i] + psi_sxy_x[j][i]; } /* right boundary */ if((!BOUNDARY) && (POS[1]==NPROCX-1) && (i>=nx2-FW+1)){ h1 = (i-nx2+2*FW); h=i; /*psi_sxx_x[j][i] = b_x_half[i] * psi_sxx_x[j][i] + a_x_half[i] * sxx_x; sxx_x = sxx_x / K_x_half[i] + psi_sxx_x[j][i];*/ psi_sxx_x[j][h1] = b_x_half[h1] * psi_sxx_x[j][h1] + a_x_half[h1] * sxx_x; sxx_x = sxx_x / K_x_half[h1] + psi_sxx_x[j][h1]; psi_sxy_x[j][h1] = b_x[h1] * psi_sxy_x[j][h1] + a_x[h1] * sxy_x; sxy_x = sxy_x / K_x[h1] + psi_sxy_x[j][h1]; } /* top boundary */ if((POS[2]==0) && (!(FREE_SURF)) && (j<=FW)){ psi_syy_y[j][i] = b_y_half[j] * psi_syy_y[j][i] + a_y_half[j] * syy_y; syy_y = syy_y / K_y_half[j] + psi_syy_y[j][i]; psi_sxy_y[j][i] = b_y[j] * psi_sxy_y[j][i] + a_y[j] * sxy_y; sxy_y = sxy_y / K_y[j] + psi_sxy_y[j][i]; } /* bottom boundary */ if((POS[2]==NPROCY-1) && (j>=ny2-FW+1)){ h1 = (j-ny2+2*FW); h = j; /*psi_syy_y[j][i] = b_y_half[j] * psi_syy_y[j][i] + a_y_half[j] * syy_y; syy_y = syy_y / K_y_half[j] + psi_syy_y[j][i];*/ psi_syy_y[h1][i] = b_y_half[h1] * psi_syy_y[h1][i] + a_y_half[h1] * syy_y; syy_y = syy_y / K_y_half[h1] + psi_syy_y[h1][i]; psi_sxy_y[h1][i] = b_y[h1] * psi_sxy_y[h1][i] + a_y[h1] * sxy_y; sxy_y = sxy_y / K_y[h1] + psi_sxy_y[h1][i]; } if(sw==0){ if (VELOCITY==0){ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH;} else{ vxp1[j][i] = rip[j][i]*(sxx_x+sxy_y)/DH; vyp1[j][i] = rjp[j][i]*(sxy_x+syy_y)/DH;} } if(sw==1){ if (VELOCITY==0){ vxp1[j][i] += vx[j][i]*DT; vyp1[j][i] += vy[j][i]*DT;} else{ vxp1[j][i] = vx[j][i]; vyp1[j][i] = vy[j][i];} } vx[j][i] += DT*rip[j][i]*(sxx_x+sxy_y)/DH; vy[j][i] += DT*rjp[j][i]*(sxy_x+syy_y)/DH; } } break; case 10: for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ vx[j][i]+= ( hc[1]*(sxx[j][i+1]-sxx[j][i] ) + hc[2]*(sxx[j][i+2]-sxx[j][i-1]) + hc[3]*(sxx[j][i+3]-sxx[j][i-2]) + hc[4]*(sxx[j][i+4]-sxx[j][i-3]) + hc[5]*(sxx[j][i+5]-sxx[j][i-4]) + hc[1]*(sxy[j][i] -sxy[j-1][i]) + hc[2]*(sxy[j+1][i]-sxy[j-2][i]) + hc[3]*(sxy[j+2][i]-sxy[j-3][i]) + hc[4]*(sxy[j+3][i]-sxy[j-4][i]) + hc[5]*(sxy[j+4][i]-sxy[j-5][i]) )*dtdh*rip[j][i]; vy[j][i]+= ( hc[1]*(syy[j+1][i]-syy[j][i] ) + hc[2]*(syy[j+2][i]-syy[j-1][i]) + hc[3]*(syy[j+3][i]-syy[j-2][i]) + hc[4]*(syy[j+4][i]-syy[j-3][i]) + hc[5]*(syy[j+5][i]-syy[j-4][i]) + hc[1]*(sxy[j][i] -sxy[j][i-1]) + hc[2]*(sxy[j][i+1]-sxy[j][i-2]) + hc[3]*(sxy[j][i+2]-sxy[j][i-3]) + hc[4]*(sxy[j][i+3]-sxy[j][i-4]) + hc[5]*(sxy[j][i+4]-sxy[j][i-5]) )*dtdh*rjp[j][i]; } } break; case 12: for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ vx[j][i]+= ( hc[1]*(sxx[j][i+1]-sxx[j][i] ) + hc[2]*(sxx[j][i+2]-sxx[j][i-1]) + hc[3]*(sxx[j][i+3]-sxx[j][i-2]) + hc[4]*(sxx[j][i+4]-sxx[j][i-3]) + hc[5]*(sxx[j][i+5]-sxx[j][i-4]) + hc[6]*(sxx[j][i+6]-sxx[j][i-5]) + hc[1]*(sxy[j][i] -sxy[j-1][i]) + hc[2]*(sxy[j+1][i]-sxy[j-2][i]) + hc[3]*(sxy[j+2][i]-sxy[j-3][i]) + hc[4]*(sxy[j+3][i]-sxy[j-4][i]) + hc[5]*(sxy[j+4][i]-sxy[j-5][i]) + hc[6]*(sxy[j+5][i]-sxy[j-6][i]) )*dtdh*rip[j][i]; vy[j][i]+= ( hc[1]*(syy[j+1][i]-syy[j][i] ) + hc[2]*(syy[j+2][i]-syy[j-1][i]) + hc[3]*(syy[j+3][i]-syy[j-2][i]) + hc[4]*(syy[j+4][i]-syy[j-3][i]) + hc[5]*(syy[j+5][i]-syy[j-4][i]) + hc[6]*(syy[j+6][i]-syy[j-5][i]) + hc[1]*(sxy[j][i] -sxy[j][i-1]) + hc[2]*(sxy[j][i+1]-sxy[j][i-2]) + hc[3]*(sxy[j][i+2]-sxy[j][i-3]) + hc[4]*(sxy[j][i+3]-sxy[j][i-4]) + hc[5]*(sxy[j][i+4]-sxy[j][i-5]) + hc[6]*(sxy[j][i+5]-sxy[j][i-6]) )*dtdh*rjp[j][i]; } } break; default: for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ vxtmp = 0; vytmp = 0; for (m=1; m<=fdoh; m++) { vxtmp += hc[m]*( sxx[j][i+m] - sxx[j][i-m+1] ) + hc[m]*( sxy[j+m-1][i] - sxy[j-m][i] ); vytmp += hc[m]*( syy[j+m][i] - syy[j-m+1][i] ) + hc[m]*( sxy[j][i+m-1] - sxy[j][i-m] ); } vx[j][i] += vxtmp*dtdh*rip[j][i]; vy[j][i] += vytmp*dtdh*rjp[j][i]; } } break; } /* end of switch(FDORDER) */ /* Forward Modelling (sw==0) */ if(sw==0){ for (l=1;l<=nsrc;l++) { i=(int)srcpos_loc[1][l]; j=(int)srcpos_loc[2][l]; azi_rad=srcpos_loc[7][l]*PI/180; SOURCE_TYPE=(int)srcpos_loc[8][l]; if(SOURCE_TYPE==2){vx[j][i] += (DT*rip[j][i]*signals1[l][nt])/(DH*DH);} /* single force in x */ if(SOURCE_TYPE==3){vy[j][i] += (DT*rjp[j][i]*signals1[l][nt])/(DH*DH);} /* single force in y */ if(SOURCE_TYPE==4){vx[j][i] += (DT*rip[j][i]*sin(azi_rad) * signals[l][nt])/(DH*DH); /* rotated force in x */ vy[j][i] += (DT*rjp[j][i]*cos(azi_rad) * signals[l][nt])/(DH*DH);} /* rotated force in y */ }} /* Backpropagation (sw==1) */ if(sw==1){ for (l=1;l<=nsrc;l++) { i=(int)srcpos_loc[1][l]; j=(int)srcpos_loc[2][l]; if(ADJOINT_TYPE==1){vx[j][i] += signals[l][nt]; /* single force in x */ vy[j][i] += signals1[l][nt];} /* + single force in y */ if(ADJOINT_TYPE==2){vy[j][i] += signals1[l][nt];} /* single force in y */ if(ADJOINT_TYPE==3){vx[j][i] += signals[l][nt];} /* single force in x */ }} /*if (FW>0.0) for (j=ny1;j<=ny2;j++){ for (i=nx1;i<=nx2;i++){ vx[j][i]*=absorb_coeff[j][i]; vy[j][i]*=absorb_coeff[j][i]; sxy[j][i]*=absorb_coeff[j][i]; sxx[j][i]*=absorb_coeff[j][i]; syy[j][i]*=absorb_coeff[j][i]; }}*/ if (infoout && (MYID==0)){ time2=MPI_Wtime(); fprintf(FP," finished (real time: %4.2f s).\n",time2-time1); } }