29 |
#define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.) |
#define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.) |
30 |
#define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b))) |
#define MINMOD(a,b) (FLUX_S(a,b)*MIN(ABS(a),ABS(b))) |
31 |
#define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b)))) |
#define SUPERBEE(a,b) (FLUX_S(a,b)*MAX(MIN(2*ABS(a),ABS(b)),MIN(ABS(a),2*ABS(b)))) |
32 |
|
#define MC(a,b) (FLUX_S(a,b)*MIN3(ABS((a)+(b))/2,2*ABS(a),2*ABS(b))) |
33 |
|
|
34 |
#define FLUX_L(a,b) SUPERBEE(a,b) /* alter for other flux limiter */ |
#define FLUX_L(a,b) MC(a,b) /* alter for other flux limiter */ |
35 |
|
|
36 |
#define FLUX_LIMITER(a) FLUX_L(a,1) |
#define FLUX_LIMITER(a) FLUX_L(1,a) |
37 |
|
|
38 |
/**************************************************************/ |
/**************************************************************/ |
39 |
|
|
125 |
|
|
126 |
void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) { |
void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) { |
127 |
|
|
128 |
register double u_i,P_p,P_n,Q_p,Q_n,r_p,r_n, a_ij, d, u_j, r_ij, f_ij, a_ji; |
register double u_i,P_p,P_n,Q_p,Q_n,r_p,r_n, a_ij, d, u_j, r_ij, f_ij, a_ji, d_ij; |
129 |
index_t color, iptr_ij,j,iptr_ji, i; |
index_t color, iptr_ij,j,iptr_ji, i; |
130 |
dim_t n; |
dim_t n; |
131 |
|
|
138 |
#pragma omp parallel private(color) |
#pragma omp parallel private(color) |
139 |
{ |
{ |
140 |
for (color=0;color<fc->num_colors;++color) { |
for (color=0;color<fc->num_colors;++color) { |
141 |
#pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,r_p,r_n,iptr_ij,a_ij,d,j,iptr_ji, u_j, r_ij, f_ij, a_ji) |
#pragma omp for schedule(static) private(i, u_i,P_p,P_n,Q_p,Q_n,r_p,r_n,iptr_ij,a_ij,d,j,iptr_ji, u_j, r_ij, f_ij, a_ji, d_ij) |
142 |
for (i = 0; i < n; ++i) { |
for (i = 0; i < n; ++i) { |
143 |
if (fc->colorOf[i]==color) { |
if (fc->colorOf[i]==color) { |
144 |
u_i=u[i]; |
u_i=u[i]; |
147 |
P_n=0.; |
P_n=0.; |
148 |
Q_p=0.; |
Q_p=0.; |
149 |
Q_n=0.; |
Q_n=0.; |
150 |
#pragma ivdep |
r_p=0.; |
151 |
|
r_n=0.; |
152 |
|
/* #pragma ivdep */ |
153 |
for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) { |
for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) { |
154 |
a_ij=fc->flux_matrix->mainBlock->val[iptr_ij]; |
a_ij=fc->flux_matrix->mainBlock->val[iptr_ij]; |
155 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
156 |
d=u[j]-u_i; |
d=u[j]-u_i; |
157 |
|
printf("%d %d : %e %e :: %e \n",i,j,u_i,u[j],a_ij); |
158 |
if (a_ij<0.) { |
if (a_ij<0.) { |
159 |
if (d<0.) { |
if (d<0.) { |
160 |
P_p+=a_ij*d; |
P_p+=a_ij*d; |
174 |
a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij]; |
a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij]; |
175 |
j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij]; |
j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij]; |
176 |
d=u_remote[j]-u_i; |
d=u_remote[j]-u_i; |
177 |
|
|
178 |
if (a_ij<0.) { |
if (a_ij<0.) { |
179 |
if (d<0.) { |
if (d<0.) { |
180 |
P_p+=a_ij*d; |
P_p+=a_ij*d; |
191 |
} |
} |
192 |
/* set the smoothness indicators */ |
/* set the smoothness indicators */ |
193 |
r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.; |
r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.; |
194 |
r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.; |
r_n = (P_n < 0.) ? FLUX_LIMITER(Q_n/P_n) : 0.; |
195 |
|
printf("%d: %e %e %e : %e %e %e : %e\n",i,Q_p,P_p,r_p,Q_n,P_n,r_n,u_i); |
196 |
/* anti diffusive flux from main block */ |
/* anti diffusive flux from main block */ |
197 |
for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) { |
for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) { |
198 |
a_ij=fc->flux_matrix->mainBlock->val[iptr_ij]; |
a_ij=fc->flux_matrix->mainBlock->val[iptr_ij]; |
199 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
200 |
if (a_ij < 0 && i!=j) { |
if ( i!=j ) { |
201 |
/* find entry a[j,i] */ |
/* find entry a[j,i] */ |
202 |
for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) { |
for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) { |
203 |
if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) { |
if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) { |
204 |
a_ji=fc->flux_matrix->mainBlock->val[iptr_ji]; |
a_ji=fc->flux_matrix->mainBlock->val[iptr_ji]; |
205 |
if (a_ji > a_ij || (a_ji == a_ij && j<i) ) { |
d_ij=-MIN3(0,a_ji,a_ij); |
206 |
|
if ( (d_ij > 0.) && ( (a_ji > a_ij) || ( (a_ji == a_ij) && (j<i) ) ) ) { |
207 |
u_j=u[j]; |
u_j=u[j]; |
208 |
r_ij = u_i>u_j ? r_p : r_n; |
r_ij = u_i>u_j ? r_p : r_n; |
209 |
f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j); |
f_ij =MIN(r_ij*d_ij,a_ji+d_ij)*(u_i-u_j); |
210 |
|
|
211 |
fa[i]+=f_ij; |
fa[i]+=f_ij; |
212 |
fa[j]-=f_ij; |
fa[j]-=f_ij; |
213 |
break; |
printf("%d %d => %e %e : %e %e : %e %e : fa[%d]=%e fa[%d]=%e\n",i,j,d_ij,(u_i-u_j), a_ij, a_ji, r_ij,f_ij,i,fa[i],j,fa[j]); |
214 |
|
|
215 |
|
|
216 |
} |
} |
217 |
|
break; |
218 |
} |
} |
219 |
} |
} |
220 |
} |
} |