From 08b12eaf4c1addf64e2bb3001ceae9539c808d49 Mon Sep 17 00:00:00 2001
From: Jan Thorbecke <janth@xs4all.nl>
Date: Fri, 5 Oct 2018 15:54:15 +0200
Subject: [PATCH] fixed OpenMP bug in viscoacoustic4.c for Threads>1

---
 fdelmodc/demo/fdelmodc_visco.scr | 1 +
 fdelmodc/fdelmodc.c              | 4 +---
 fdelmodc/viscoacoustic4.c        | 8 +-------
 3 files changed, 3 insertions(+), 10 deletions(-)

diff --git a/fdelmodc/demo/fdelmodc_visco.scr b/fdelmodc/demo/fdelmodc_visco.scr
index f43ec8d..ecdcbb7 100755
--- a/fdelmodc/demo/fdelmodc_visco.scr
+++ b/fdelmodc/demo/fdelmodc_visco.scr
@@ -37,6 +37,7 @@ export filero=model_ro.su
 export fileqp=relax_cp.su
 export fileqs=relax_cs.su
 
+export OMP_NUM_THREADS=4
 
 ../fdelmodc \
 	file_cp=$filecp file_cs=$filecs file_den=$filero \
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 2510135..1ff284c 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -40,9 +40,7 @@ int acoustic2(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixs
 int acoustic4Block(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx,
 float *vz, float *p, float *rox, float *roz, float *l2m, int verbose);
 
-
-int viscoacoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx, float
-*vz, float *p, float *rox, float *roz, float *l2m, float *tss, float *tep, float *q, int verbose);
+int viscoacoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx, float *vz, float *p, float *rox, float *roz, float *l2m, float *tss, float *tep, float *q, int verbose);
 
 int elastic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int verbose);
 
diff --git a/fdelmodc/viscoacoustic4.c b/fdelmodc/viscoacoustic4.c
index 91601cd..63207c7 100644
--- a/fdelmodc/viscoacoustic4.c
+++ b/fdelmodc/viscoacoustic4.c
@@ -112,14 +112,13 @@ int viscoacoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, in
 	boundariesP(mod, bnd, vx, vz, p, NULL, NULL, rox, roz, l2m, NULL, NULL, itime, verbose);
 
 	/* calculate p/tzz for all grid points except on the virtual boundary */
+#pragma omp	for private (iz,ix, Tpp) nowait schedule(guided,1)
 	for (ix=mod.ioPx; ix<mod.iePx; ix++) {
-#pragma omp	for private (iz) nowait schedule(guided,1)
 #pragma ivdep
 		for (iz=mod.ioPz; iz<mod.iePz; iz++) {
 			dxvx[iz] = c1*(vx[(ix+1)*n1+iz] - vx[ix*n1+iz]) +
 					   c2*(vx[(ix+2)*n1+iz] - vx[(ix-1)*n1+iz]);
 		}
-#pragma omp	for private (iz) nowait schedule(guided,1)
 #pragma ivdep
 		for (iz=mod.ioPz; iz<mod.iePz; iz++) {
 			dzvz[iz] = c1*(vz[ix*n1+iz+1]   - vz[ix*n1+iz]) +
@@ -127,14 +126,12 @@ int viscoacoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, in
 		}
 
 		/* help variables to let the compiler vectorize the loops */
-#pragma omp	for private (iz) nowait schedule(guided,1)
 #pragma ivdep
 		for (iz=mod.ioPz; iz<mod.iePz; iz++) {
 			Tpp     = tep[ix*n1+iz]*tss[ix*n1+iz];
 			Tlm[iz] = (1.0-Tpp)*tss[ix*n1+iz]*l2m[ix*n1+iz]*0.5;
 			Tlp[iz] = l2m[ix*n1+iz]*Tpp;
 		}   
-#pragma omp	for private (iz)  schedule(guided,1)
 #pragma ivdep
 		for (iz=mod.ioPz; iz<mod.iePz; iz++) {
 			Tt1[iz] = 1.0/(ddt+0.5*tss[ix*n1+iz]);
@@ -142,18 +139,15 @@ int viscoacoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, in
 		}   
 
 		/* the update with the relaxation correction */
-#pragma omp	for private (iz)  schedule(guided,1)
 #pragma ivdep
 		for (iz=mod.ioPz; iz<mod.iePz; iz++) {
 			p[ix*n1+iz] -= Tlp[iz]*(dzvz[iz]+dxvx[iz]) + q[ix*n1+iz];
 		}
-#pragma omp	for private (iz)  schedule(guided,1)
 #pragma ivdep
 		for (iz=mod.ioPz; iz<mod.iePz; iz++) {
 			q[ix*n1+iz] = (Tt2[iz]*q[ix*n1+iz] + Tlm[iz]*(dxvx[iz]+dzvz[iz]))*Tt1[iz];
 			p[ix*n1+iz] -= q[ix*n1+iz];
 		}
-
 	}
 
 	/* Add stress source */
-- 
GitLab