solver
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
void iterate_hyd_p() | |
{ | |
int i,j,k, mythread; | |
double dummy; | |
normres_sparse = 0.0; | |
pipi_sparse = 0.0; | |
riri_sparse2 = 0.0; | |
# pragma omp parallel num_threads(NTt) default(none) private(i,j,k, mythread, dummy) shared(STA,res_sparse_s,COEFF,p_sparse_s, ap_sparse_s,h_sparse_s,RLL, pipi_sparse, normres_sparse, riri_sparse,riri_sparse2,noemer_sparse, nx, ny, nz, nv, PeriodicBoundaryX, PeriodicBoundaryY, PeriodicBoundaryZ) | |
{ | |
mythread = omp_get_thread_num();//0 | |
// loop 1 | |
#pragma omp for reduction(+:pipi_sparse) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
dummy = COEFF[i][j][k][6] * p_sparse_s[i][j][k]; | |
if (PeriodicBoundaryX && i == 1) dummy += COEFF[i][j][k][0] * p_sparse_s[nx ][j][k]; | |
else dummy += COEFF[i][j][k][0] * p_sparse_s[i-1][j][k]; | |
if (PeriodicBoundaryX && i == nx) dummy += COEFF[i][j][k][1] * p_sparse_s[1 ][j][k]; | |
else dummy += COEFF[i][j][k][1] * p_sparse_s[i+1][j][k]; | |
if (PeriodicBoundaryY && j == 1) dummy += COEFF[i][j][k][2] * p_sparse_s[i][ny ][k]; | |
else dummy += COEFF[i][j][k][2] * p_sparse_s[i][j-1][k]; | |
if (PeriodicBoundaryY && j == ny) dummy += COEFF[i][j][k][3] * p_sparse_s[i][ 1][k]; | |
else dummy += COEFF[i][j][k][3] * p_sparse_s[i][j+1][k]; | |
if (PeriodicBoundaryZ && k == 1) dummy += COEFF[i][j][k][4] * p_sparse_s[i][j][nz ]; | |
else dummy += COEFF[i][j][k][4] * p_sparse_s[i][j][k-1]; | |
if (PeriodicBoundaryZ && k == nz) dummy += COEFF[i][j][k][5] * p_sparse_s[i][j][ 1]; | |
else dummy += COEFF[i][j][k][5] * p_sparse_s[i][j][k+1]; | |
ap_sparse_s[i][j][k] = dummy; | |
pipi_sparse += p_sparse_s[i][j][k] * ap_sparse_s[i][j][k]; | |
} | |
// loop 2 | |
#pragma omp for reduction(+:normres_sparse) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
STA[i][j][k] += (riri_sparse / pipi_sparse) * p_sparse_s[i][j][k]; | |
res_sparse_s[i][j][k] -= (riri_sparse / pipi_sparse) * ap_sparse_s[i][j][k]; | |
normres_sparse += (res_sparse_s[i][j][k] * res_sparse_s[i][j][k])/ nv;// need to take square root separately | |
} | |
// loop 3 | |
// FORWARD | |
#pragma omp for schedule(static, nx/NTt) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
dummy = res_sparse_s[i][j][k]; | |
dummy -= COEFF[i][j][k][7] * RLL[i-1][j][k]; | |
if (PeriodicBoundaryX && i==nx)dummy -= COEFF[i][j][k][8] * RLL[1 ][j][k]; | |
dummy -= COEFF[i][j][k][2] * RLL[i][j-1][k]; | |
if (PeriodicBoundaryY && j==ny) dummy -= COEFF[i][j][k][3] * RLL[i][1 ][k]; | |
dummy -= COEFF[i][j][k][4] * RLL[i][j][k-1]; | |
if (PeriodicBoundaryZ && k==nz) dummy -= COEFF[i][j][k][5] * RLL[i][j][1 ]; | |
RLL[i][j][k] = dummy / h_sparse_s[i][j][k]; | |
} | |
// loop 4 | |
// BACKWARD | |
#pragma omp for schedule(static, nx/NTt) | |
for (i=nx; i>=1;i--) for (j=ny; j>=1;j--) for (k=nz; k>=1;k--) | |
{ | |
dummy = RLL[i][j][k]*h_sparse_s[i][j][k]; | |
if (PeriodicBoundaryX && i==1) dummy -= COEFF[i][j][k][7] * RLL[nx ][j][k]; | |
dummy -= COEFF[i][j][k][8] * RLL[i+1][j][k]; | |
if (PeriodicBoundaryY && j==1) dummy -= COEFF[i][j][k][2] * RLL[i][ny ][k]; | |
dummy -= COEFF[i][j][k][3] * RLL[i][j+1][k]; | |
if (PeriodicBoundaryZ && k==1) dummy -= COEFF[i][j][k][4] * RLL[i][j][nz ]; | |
dummy -= COEFF[i][j][k][5] * RLL[i][j][k+1]; | |
RLL[i][j][k] = dummy / h_sparse_s[i][j][k]; | |
} | |
// loop 5 | |
#pragma omp for reduction(+:riri_sparse2) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
riri_sparse2 += res_sparse_s[i][j][k] * RLL[i][j][k]; | |
} | |
// loop 6 | |
#pragma omp for | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
p_sparse_s[i][j][k] = RLL[i][j][k] + (riri_sparse2 / noemer_sparse) * p_sparse_s[i][j][k]; | |
} | |
} | |
noemer_sparse = riri_sparse2; | |
riri_sparse = riri_sparse2; | |
return; | |
} | |
void initiate_hyd_p() | |
{ | |
int i,j,k, mythread; | |
double dummy; | |
normres_sparse = 0.0; | |
riri_sparse = 0.0; | |
// loop 1 | |
for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
COEFF[1 ][j][k][7] = 0.0;// -ve X | |
COEFF[nx][j][k][8] = 0.0;// +ve X | |
} | |
# pragma omp parallel num_threads(NTt) default(none) private(i,j,k, mythread, dummy) shared(ap_sparse_s, STA, res_sparse_s, COEFF, RLL, p_sparse_s, h_sparse_s, normres_sparse, riri_sparse, nx, ny, nz, nv, PeriodicBoundaryX, PeriodicBoundaryY, PeriodicBoundaryZ) | |
{ | |
mythread = omp_get_thread_num();//0 | |
// loop 2 | |
if(mythread != NTt-1)// if NTt no. of processor division = NTt - 1 | |
{ | |
for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
COEFF[(mythread + 1)*(nx/NTt) + 1][j][k][7] = 0.0;// -ve X | |
COEFF[(mythread + 1)*(nx/NTt) ][j][k][8] = 0.0;// +ve X | |
} | |
} | |
// loop 3 | |
#pragma omp for | |
for (i = 0; i <= nx+1; i++) for (j = 0; j <= ny+1; j++) for (k = 0; k <= nz+1; k++) | |
{ | |
h_sparse_s[i][j][k] = 1.0; | |
res_sparse_s[i][j][k] = 0.0; | |
p_sparse_s[i][j][k] = 0.0; | |
ap_sparse_s[i][j][k] = 0.0; | |
} | |
// loop 4 | |
#pragma omp for schedule(static, nx/NTt) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
dummy = COEFF[i ][j][k][6]; | |
if (PeriodicBoundaryX && i==nx) dummy -= COEFF[nx][j][k][8] * COEFF[1 ][j][k][7] / h_sparse_s[1 ][j][k]; | |
dummy -= COEFF[i ][j][k][7] * COEFF[i-1][j][k][8] / h_sparse_s[i-1][j][k]; | |
if (PeriodicBoundaryY && j==ny) dummy -= COEFF[i][ny][k][3] * COEFF[i][1 ][k][2] / h_sparse_s[i][1 ][k]; | |
dummy -= COEFF[i][j ][k][2] * COEFF[i][j-1][k][3] / h_sparse_s[i][j-1][k]; | |
if (PeriodicBoundaryZ && k==nz) dummy -= COEFF[i][j][nz][5] * COEFF[i][j][1 ][4] / h_sparse_s[i][j][1 ]; | |
dummy -= COEFF[i][j][k ][4] * COEFF[i][j][k-1][5] / h_sparse_s[i][j][k-1]; | |
h_sparse_s[i][j][k] = dummy; | |
} | |
// loop 5 | |
#pragma omp for reduction(+:normres_sparse) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
dummy = COEFF[i][j][k][6] * STA[i][j][k]; | |
if (PeriodicBoundaryX && i == 1) dummy += COEFF[i][j][k][0] * STA[nx ][j][k]; | |
else dummy += COEFF[i][j][k][0] * STA[i-1][j][k]; | |
if (PeriodicBoundaryX && i == nx) dummy += COEFF[i][j][k][1] * STA[1 ][j][k]; | |
else dummy += COEFF[i][j][k][1] * STA[i+1][j][k]; | |
if (PeriodicBoundaryY && j == 1) dummy += COEFF[i][j][k][2] * STA[i][ny ][k]; | |
else dummy += COEFF[i][j][k][2] * STA[i][j-1][k]; | |
if (PeriodicBoundaryY && j == ny) dummy += COEFF[i][j][k][3] * STA[i][1 ][k]; | |
else dummy += COEFF[i][j][k][3] * STA[i][j+1][k]; | |
if (PeriodicBoundaryZ && k == 1) dummy += COEFF[i][j][k][4] * STA[i][j][nz ]; | |
else dummy += COEFF[i][j][k][4] * STA[i][j][k-1]; | |
if (PeriodicBoundaryZ && k == nz) dummy += COEFF[i][j][k][5] * STA[i][j][1 ]; | |
else dummy += COEFF[i][j][k][5] * STA[i][j][k+1]; | |
p_sparse_s[i][j][k] = dummy; | |
res_sparse_s[i][j][k] = RLL[i][j][k] + (-1.0) * p_sparse_s[i][j][k]; //coef_a_sparse = -1.0; | |
normres_sparse += (res_sparse_s[i][j][k] * res_sparse_s[i][j][k])/nv; // need to take square separately | |
} | |
// loop 6 | |
// FORWARD | |
#pragma omp for schedule(static, nx/NTt) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
dummy = res_sparse_s[i][j][k]; | |
dummy -= COEFF[i][j][k][7] * p_sparse_s[i-1][j][k]; | |
if (PeriodicBoundaryX && i==nx) dummy -= COEFF[i][j][k][8] * p_sparse_s[1 ][j][k]; | |
dummy -= COEFF[i][j][k][2] * p_sparse_s[i][j-1][k]; | |
if (PeriodicBoundaryY && j==ny) dummy -= COEFF[i][j][k][3] * p_sparse_s[i][1 ][k]; | |
dummy -= COEFF[i][j][k][4] * p_sparse_s[i][j][k-1]; | |
if (PeriodicBoundaryZ && k==nz) dummy -= COEFF[i][j][k][5] * p_sparse_s[i][j][1 ]; | |
p_sparse_s[i][j][k] = dummy / h_sparse_s[i][j][k]; | |
} | |
// loop 7 | |
// BACKWARD | |
#pragma omp for schedule(static, nx/NTt) | |
for (i=nx; i>=1;i--) for (j=ny; j>=1;j--) for (k=nz; k>=1;k--) | |
{ | |
dummy = p_sparse_s[i][j][k]*h_sparse_s[i][j][k]; | |
if (PeriodicBoundaryX && i==1) dummy -= COEFF[i][j][k][7] * p_sparse_s[nx ][j][k]; | |
dummy -= COEFF[i][j][k][8] * p_sparse_s[i+1][j][k]; | |
if (PeriodicBoundaryY && j==1) dummy -= COEFF[i][j][k][2] * p_sparse_s[i][ny ][k]; | |
dummy -= COEFF[i][j][k][3] * p_sparse_s[i][j+1][k]; | |
if (PeriodicBoundaryZ && k==1) dummy -= COEFF[i][j][k][4] * p_sparse_s[i][j][nz ]; | |
dummy -= COEFF[i][j][k][5] * p_sparse_s[i][j][k+1]; | |
p_sparse_s[i][j][k] = dummy / h_sparse_s[i][j][k]; | |
} | |
// loop 8 | |
# pragma omp for reduction(+:riri_sparse) | |
for (i=1; i<=nx; i++) for (j=1; j<=ny; j++) for (k=1; k<=nz; k++) | |
{ | |
riri_sparse += res_sparse_s[i][j][k]* p_sparse_s[i][j][k]; | |
} | |
} | |
noemer_sparse = riri_sparse; | |
return; | |
} | |
void biccg_periodic(double convg_criteria) | |
{ | |
int ite_start; | |
initiate_hyd_p(); | |
ite_start = ite_hydro; | |
while ((sqrt(normres_sparse) > convg_criteria*MatScale) && (ite_hydro -ite_start <= itm_icg)) | |
{ | |
ite_hydro++; | |
iterate_hyd_p(); | |
} | |
return; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment