Commit 3b01f88e authored by Tilman Steinweg's avatar Tilman Steinweg

Bugfix

To exchange the wavefield only a wrapping layer of FDORDER/2 Gridpoints is necessary.
The FDORDER/2+1 layer was probably implemented because it fixed indirectly the bug (inconsistency)
in the exchange_s/v functions.
parent 8cdd1868
......@@ -38,7 +38,7 @@ MPI_Request *req_send, MPI_Request *req_rec){
/* comunication initialisation for persistent communication */
fdo2 = 2*(FDORDER/2 + 1);
fdo2 = 2*(FDORDER/2);
/* buffer arrays are copied into local buffers using buffered send (bsend),
......
......@@ -32,14 +32,13 @@ void exchange_s(int nd, float ** sxx, float ** syy,
extern int NX, NY, POS[3], NPROCX, NPROCY, BOUNDARY, FDORDER;
extern int INDEX[5];
extern int INDEX[5],MYID;
extern const int TAG1,TAG2,TAG5,TAG6;
MPI_Status status;
int i, j, fdo, fdo2, n, l;
fdo = nd;
// fdo2 = FDORDER + 1;
fdo2 = 2*fdo;
......@@ -69,8 +68,6 @@ void exchange_s(int nd, float ** sxx, float ** syy,
bufferbot_to_top[i][n++] = syy[NY-l+1][i];
}
}
/* send and receive values for points at inner boundaries */
/* blocking communication */
/*
......@@ -117,16 +114,17 @@ void exchange_s(int nd, float ** sxx, float ** syy,
}
/* left - right */
if ((BOUNDARY) || (POS[1]!=0)) /* no boundary exchange at left edge of global grid */
for (j=1;j<=NY;j++){
/* storage of left edge of local volume into buffer */
n = 1;
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
bufferlef_to_rig[j][n++] = sxy[j][l];
}
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
bufferlef_to_rig[j][n++] = sxx[j][l];
}
}
......@@ -136,10 +134,10 @@ void exchange_s(int nd, float ** sxx, float ** syy,
for (j=1;j<=NY;j++){
/* storage of right edge of local volume into buffer */
n = 1;
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
bufferrig_to_lef[j][n++] = sxy[j][NX-l+1];
}
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
bufferrig_to_lef[j][n++] = sxx[j][NX-l+1];
}
}
......@@ -173,10 +171,10 @@ void exchange_s(int nd, float ** sxx, float ** syy,
if ((BOUNDARY) || (POS[1]!=NPROCX-1)) /* no boundary exchange at right edge of global grid */
for (j=1;j<=NY;j++){
n = 1;
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
sxy[j][NX+l] = bufferlef_to_rig[j][n++];
}
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
sxx[j][NX+l] = bufferlef_to_rig[j][n++];
}
}
......@@ -184,10 +182,10 @@ void exchange_s(int nd, float ** sxx, float ** syy,
if ((BOUNDARY) || (POS[1]!=0)) /* no boundary exchange at left edge of global grid */
for (j=1;j<=NY;j++){
n = 1;
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
sxy[j][1-l] = bufferrig_to_lef[j][n++];
}
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
sxx[j][1-l] = bufferrig_to_lef[j][n++];
}
}
......
......@@ -40,7 +40,6 @@ void exchange_v(int nd, float ** vx, float ** vy,
fdo = nd;
// fdo2= FDORDER + 1;
fdo2 = 2*fdo;
/* top - bottom */
......@@ -127,10 +126,10 @@ void exchange_v(int nd, float ** vx, float ** vy,
for (j=1;j<=NY;j++){
/* storage of left edge of local volume into buffer */
n = 1;
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
bufferlef_to_rig[j][n++] = vy[j][l];
}
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
bufferlef_to_rig[j][n++] = vx[j][l];
}
}
......@@ -141,10 +140,10 @@ void exchange_v(int nd, float ** vx, float ** vy,
for (j=1;j<=NY;j++){
/* storage of right edge of local volume into buffer */
n = 1;
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
bufferrig_to_lef[j][n++] = vy[j][NX-l+1];
}
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
bufferrig_to_lef[j][n++] = vx[j][NX-l+1];
}
}
......@@ -178,10 +177,10 @@ void exchange_v(int nd, float ** vx, float ** vy,
if ((BOUNDARY) || (POS[1]!=NPROCX-1)) /* no boundary exchange at right edge of global grid */
for (j=1;j<=NY;j++){
n = 1;
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
vy[j][NX+l] = bufferlef_to_rig[j][n++];
}
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
vx[j][NX+l] = bufferlef_to_rig[j][n++];
}
}
......@@ -190,10 +189,10 @@ void exchange_v(int nd, float ** vx, float ** vy,
if ((BOUNDARY) || (POS[1]!=0)) /* no boundary exchange at left edge of global grid */
for (j=1;j<=NY;j++){
n = 1;
for (l=1;l<fdo-1;l++) {
for (l=1;l<=fdo-1;l++) {
vy[j][1-l] = bufferrig_to_lef[j][n++];
}
for (l=1;l<fdo;l++) {
for (l=1;l<=fdo;l++) {
vx[j][1-l] = bufferrig_to_lef[j][n++];
}
}
......
......@@ -161,7 +161,7 @@ void snap(FILE *fp,int nt, int nsnap, float **vx, float **vy, float **sxx,
fpy2=fopen(snapfile_rot,"a");
}
nd = FDORDER/2+1;
nd = FDORDER/2;
curlfield = matrix(-nd+1,NY+nd,-nd+1,NX+nd);
......
......@@ -262,7 +262,7 @@ int main ( int argc, char **argv )
}
/*allocate memory for dynamic, static and buffer arrays */
nd = FDORDER/2 + 1;
nd = FDORDER/2;
fdo3 = 2 * nd;
fac1 = ( NX + fdo3 ) * ( NY + fdo3 );
......@@ -650,24 +650,21 @@ int main ( int argc, char **argv )
if(FDORDER_TIME==4){
if(L) {
zero_visco_4(-nd+1,NY+nd,-nd+1,NX+nd,vxx_1,vxx_2,vxx_3,vxx_4,vyy_1,vyy_2,vyy_3,vyy_4,vxy_1,vxy_2,vxy_3,vxy_4,vyx_1,vyx_2,vyx_3,vyx_4,svx_1,svx_2,svx_3,svx_4,svy_1,svy_2,svy_3,svy_4,pr_2,pr_3,pr_4,pp_2,pp_3,pp_4,pq_2,pq_3,pq_4);
zero_visco_4(-nd+1,NY+nd,-nd+1,NX+nd,vxx_1,vxx_2,vxx_3,vxx_4,vyy_1,vyy_2,vyy_3,vyy_4,vxy_1,vxy_2,vxy_3,vxy_4,
vyx_1,vyx_2,vyx_3,vyx_4,svx_1,svx_2,svx_3,svx_4,
svy_1,svy_2,svy_3,svy_4,pr_2,pr_3,pr_4,pp_2,pp_3,pp_4,pq_2,pq_3,pq_4);
} else {
zero_elastic_4(-nd+1,NY+nd,-nd+1,NX+nd,vxx_1,vxx_2,vxx_3,vxx_4,vyy_1,vyy_2,vyy_3,vyy_4,vxy_1,vxy_2,vxy_3,vxy_4,vyx_1,vyx_2,vyx_3,vyx_4,svx_1,svx_2,svx_3,svx_4,svy_1,svy_2,svy_3,svy_4);
zero_elastic_4(-nd+1,NY+nd,-nd+1,NX+nd,vxx_1,vxx_2,vxx_3,vxx_4,vyy_1,vyy_2,vyy_3,vyy_4,vxy_1,vxy_2,vxy_3,vxy_4,
vyx_1,vyx_2,vyx_3,vyx_4,svx_1,svx_2,svx_3,svx_4,svy_1,svy_2,svy_3,svy_4);
}
}
if ( ABS_TYPE != 1 ) {
if ( L )
zero_visc ( -FDORDER / 2, NX + FDORDER / 2 + 1, -FDORDER / 2,
NY + FDORDER / 2 + 1, pvx, pvy, psxx, psyy, psxy, pr, pp,
pq );
zero_visc ( -nd+1, NX+nd, -nd+1,NY+nd, pvx, pvy, psxx, psyy, psxy, pr, pp,pq );
else
zero_elastic ( -FDORDER / 2, NX + FDORDER / 2 + 1, 1 - FDORDER / 2,
NY + FDORDER / 2 + 1, pvx, pvy, psxx, psyy, psxy );
zero_elastic ( -nd+1, NX + nd, -nd+1, NY+nd, pvx, pvy, psxx, psyy, psxy );
}
/* above arguments for function call have been defined differently -> minus 1 in every dimension, not quite clear why
old version:
zero_elastic(1-FDORDER/2,NX+FDORDER/2,1-FDORDER/2,NY+FDORDER/2,pvx,pvy,psxx,psyy,psxy); */
/* Reseting lsmap to NDT for saving seismograms */
lsamp = NDT;
......
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