-
Notifications
You must be signed in to change notification settings - Fork 2
/
sor.c
executable file
·102 lines (84 loc) · 2.42 KB
/
sor.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#include "sor.h"
#include <math.h>
#include "parallel.h"
void sor(
double omg,
double dx,
double dy,
int il,
int ir,
int jb,
int jt,
int rank_l,
int rank_r,
int rank_b,
int rank_t,
double *bufSend,
double *bufRecv,
int chunk,
double **P,
double **RS,
int myrank,
int imax,
int jmax,
double *res,
int omg_i,
int omg_j,
int iproc,
int jproc
) {
int i,j;
double rloc;
double coeff = omg/(2.0*(1.0/(dx*dx)+1.0/(dy*dy)));
double glob_res;
/* SOR iteration */
for(j = jb; j <= jt; j++) {
for(i = il; i<=ir; i++) {
P[i][j] = (1.0-omg)*P[i][j]
+ coeff*(( P[i+1][j]+P[i-1][j])/(dx*dx) + ( P[i][j+1]+P[i][j-1])/(dy*dy) - RS[i][j]);
}
}
/*Boundary values*/
if(omg_i==1){
for (j = jb; j <= jt+1; j++){
/*P values around left boundary*/
P[il-1][j] = P[il][j];
}
}
if(omg_i==iproc){
for (j = jb; j <= jt+1; j++){
/*P values around right boundary*/
P[ir+1][j] = P[ir][j];
}
}
if(omg_j==1){
for (i = il; i <= ir+1; i++){
/*P values around bottom boundary*/
P[i][jb-1] = P[i][jb];
}
}
if(omg_j==jproc){
for (i = il; i <= ir+1; i++){
/*P values around top boundary*/
P[i][jt+1] = P[i][jt];
}
}
/*Passing the pressure values*/
pressure_comm(P,il,ir,jb,jt ,rank_l,rank_r,rank_b,rank_t,bufSend,bufRecv,chunk);
/* compute the residual */
rloc = 0;
for(j = jb; j <= jt; j++) {
for(i = il; i <= ir; i++) {
rloc += ( (P[i+1][j]-2.0*P[i][j]+P[i-1][j])/(dx*dx) + ( P[i][j+1]-2.0*P[i][j]+P[i][j-1])/(dy*dy) - RS[i][j])*
( (P[i+1][j]-2.0*P[i][j]+P[i-1][j])/(dx*dx) + ( P[i][j+1]-2.0*P[i][j]+P[i][j-1])/(dy*dy) - RS[i][j]);
}
}
MPI_Reduce(&rloc, &glob_res, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
if ( myrank == 0 ) {
glob_res = glob_res/(jmax*imax) ;
glob_res = sqrt(glob_res) ;
/*printf("my rank = %i, global = %f\n", myrank, glob_res);*/
}
MPI_Bcast (&glob_res,1,MPI_DOUBLE,0,MPI_COMM_WORLD) ;
*res=glob_res;
}