1 |
/* $Id:$ */ |
2 |
|
3 |
/******************************************************* |
4 |
* |
5 |
* Copyright 2007 by University of Queensland |
6 |
* |
7 |
* http://esscc.uq.edu.au |
8 |
* Primary Business: Queensland, Australia |
9 |
* Licensed under the Open Software License version 3.0 |
10 |
* http://www.opensource.org/licenses/osl-3.0.php |
11 |
* |
12 |
*******************************************************/ |
13 |
|
14 |
/**************************************************************/ |
15 |
|
16 |
/* Paso: FluxControl */ |
17 |
|
18 |
/**************************************************************/ |
19 |
|
20 |
/* Author: l.gross@uq.edu.au */ |
21 |
|
22 |
/**************************************************************/ |
23 |
|
24 |
|
25 |
#include "Paso.h" |
26 |
#include "SolverFCT.h" |
27 |
#include "PasoUtil.h" |
28 |
|
29 |
#define FLUX_S(a,b) ((SIGN(a)+SIGN(b))/2.) |
30 |
#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)))) |
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) MC(a,b) /* alter for other flux limiter */ |
35 |
|
36 |
#define FLUX_LIMITER(a) FLUX_L(1,a) |
37 |
|
38 |
/**************************************************************/ |
39 |
|
40 |
/* adds A plus stabelising diffusion into the matrix B */ |
41 |
|
42 |
/* d_ij=alpha*max(0,-a[i,j],-a[j,i]) */ |
43 |
/* b[i,j]+=alpha*(a[i,j]+d_ij) */ |
44 |
/* b[j,i]+=alpha*(a[j,i]+d_ij) */ |
45 |
/* b[i,i]-=alpha*d_ij */ |
46 |
/* b[j,j]-=alpha*d_ij */ |
47 |
|
48 |
void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) { |
49 |
dim_t n,i; |
50 |
index_t color, iptr_ij,j,iptr_ji; |
51 |
register double d_ij; |
52 |
|
53 |
if (fc==NULL) return; |
54 |
n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix); |
55 |
|
56 |
#pragma omp parallel private(color) |
57 |
{ |
58 |
/* process main block */ |
59 |
for (color=0;color<fc->num_colors;++color) { |
60 |
#pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij) schedule(static) |
61 |
for (i = 0; i < n; ++i) { |
62 |
if (fc->colorOf[i]==color) { |
63 |
fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]]; |
64 |
|
65 |
for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) { |
66 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
67 |
if (j<i) { |
68 |
/* find entry a[j,i] */ |
69 |
for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) { |
70 |
if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) { |
71 |
d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]); |
72 |
printf("%d %d -> %e\n",i,j,d_ij); |
73 |
fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij; |
74 |
fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij; |
75 |
printf("%d %d -> %e -> %e %e \n",i,j,d_ij,fc->transport_matrix->mainBlock->val[iptr_ij],fc->transport_matrix->mainBlock->val[iptr_ji]); |
76 |
fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij; |
77 |
fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij; |
78 |
break; |
79 |
} |
80 |
} |
81 |
} |
82 |
|
83 |
} |
84 |
/* TODO process couple block */ |
85 |
} |
86 |
} |
87 |
#pragma omp barrier |
88 |
} |
89 |
} |
90 |
|
91 |
} |
92 |
|
93 |
void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) { |
94 |
|
95 |
double *remote_u=NULL; |
96 |
|
97 |
if (fc==NULL) return; |
98 |
Paso_SystemMatrix_startCollect(fc->transport_matrix,u); |
99 |
/* process main block */ |
100 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa); |
101 |
/* finish exchange */ |
102 |
remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix); |
103 |
/* process couple block */ |
104 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa); |
105 |
|
106 |
Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa); |
107 |
} |
108 |
/**************************************************************/ |
109 |
|
110 |
/* adds antidiffusion to fa |
111 |
|
112 |
P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i]) |
113 |
P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i]) |
114 |
Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i]) |
115 |
Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i]) |
116 |
d_ij=max(0,-a[i,j],-a[j,i]) |
117 |
l_ji=max(0,a[j,i],a[j,i]-a[i,j]) |
118 |
if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij) |
119 |
r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i] |
120 |
f_ij =min(FLUX_LIMITER(r_ij)*d_ij,l_ji) (u[i]-u[j])=min(FLUX_LIMITER(r_ij)*a[i,j],a[j,i]-a[i,j]) (u[i]-u[j]) |
121 |
fa[i]+=f_ij |
122 |
fa[j]-=f_ij |
123 |
|
124 |
*/ |
125 |
|
126 |
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, d_ij; |
129 |
index_t color, iptr_ij,j,iptr_ji, i; |
130 |
dim_t n; |
131 |
|
132 |
|
133 |
if (fc==NULL) return; |
134 |
n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix); |
135 |
/* exchange */ |
136 |
|
137 |
|
138 |
#pragma omp parallel private(color) |
139 |
{ |
140 |
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, d_ij) |
142 |
for (i = 0; i < n; ++i) { |
143 |
if (fc->colorOf[i]==color) { |
144 |
u_i=u[i]; |
145 |
/* gather the smoothness sensor */ |
146 |
P_p=0.; |
147 |
P_n=0.; |
148 |
Q_p=0.; |
149 |
Q_n=0.; |
150 |
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) { |
154 |
a_ij=fc->flux_matrix->mainBlock->val[iptr_ij]; |
155 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
156 |
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.) { |
159 |
if (d<0.) { |
160 |
P_p+=a_ij*d; |
161 |
} else { |
162 |
P_n+=a_ij*d; |
163 |
} |
164 |
} else { |
165 |
if (d>0.) { |
166 |
Q_p+=a_ij*d; |
167 |
} else { |
168 |
Q_n+=a_ij*d; |
169 |
} |
170 |
} |
171 |
} |
172 |
#pragma ivdep |
173 |
for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) { |
174 |
a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij]; |
175 |
j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij]; |
176 |
d=u_remote[j]-u_i; |
177 |
|
178 |
if (a_ij<0.) { |
179 |
if (d<0.) { |
180 |
P_p+=a_ij*d; |
181 |
} else { |
182 |
P_n+=a_ij*d; |
183 |
} |
184 |
} else { |
185 |
if (d>0.) { |
186 |
Q_p+=a_ij*d; |
187 |
} else { |
188 |
Q_n+=a_ij*d; |
189 |
} |
190 |
} |
191 |
} |
192 |
/* set the smoothness indicators */ |
193 |
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.; |
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 */ |
197 |
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]; |
199 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
200 |
if ( i!=j ) { |
201 |
/* find entry a[j,i] */ |
202 |
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) { |
204 |
a_ji=fc->flux_matrix->mainBlock->val[iptr_ji]; |
205 |
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]; |
208 |
r_ij = u_i>u_j ? r_p : r_n; |
209 |
f_ij =MIN(r_ij*d_ij,a_ji+d_ij)*(u_i-u_j); |
210 |
|
211 |
fa[i]+=f_ij; |
212 |
fa[j]-=f_ij; |
213 |
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 |
} |
221 |
} |
222 |
/* anti diffusive flux from couple block */ |
223 |
|
224 |
/* TODO */ |
225 |
} |
226 |
} |
227 |
#pragma omp barrier |
228 |
} |
229 |
} |
230 |
} |