Skip to content

Instantly share code, notes, and snippets.

@data-panda
Created July 27, 2018 13:16
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save data-panda/079cfb076092a5289945c9b3b0881fa9 to your computer and use it in GitHub Desktop.
solver
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