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 |
|
33 |
#define FLUX_L(a,b) SUPERBEE(a,b) /* alter for other flux limiter */ |
34 |
|
35 |
#define FLUX_LIMITER(a) FLUX_L(a,1) |
36 |
|
37 |
/**************************************************************/ |
38 |
|
39 |
/* adds A plus stabelising diffusion into the matrix B */ |
40 |
|
41 |
/* d_ij=alpha*max(0,-a[i,j],-a[j,i]) */ |
42 |
/* b[i,j]+=alpha*(a[i,j]+d_ij) */ |
43 |
/* b[j,i]+=alpha*(a[j,i]+d_ij) */ |
44 |
/* b[i,i]-=alpha*d_ij */ |
45 |
/* b[j,j]-=alpha*d_ij */ |
46 |
|
47 |
void Paso_FCTransportProblem_addAdvectivePart(Paso_FCTransportProblem * fc, double alpha) { |
48 |
dim_t n,i; |
49 |
index_t color, iptr_ij,j,iptr_ji; |
50 |
register double d_ij; |
51 |
|
52 |
if (fc==NULL) return; |
53 |
n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix); |
54 |
|
55 |
#pragma omp parallel private(color) |
56 |
{ |
57 |
/* process main block */ |
58 |
for (color=0;color<fc->num_colors;++color) { |
59 |
#pragma omp for private(i,iptr_ij,j,iptr_ji,d_ij) schedule(static) |
60 |
for (i = 0; i < n; ++i) { |
61 |
if (fc->colorOf[i]==color) { |
62 |
fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]+=alpha*fc->flux_matrix->mainBlock->val[fc->main_iptr[i]]; |
63 |
|
64 |
for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) { |
65 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
66 |
if (j<i) { |
67 |
/* find entry a[j,i] */ |
68 |
for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[j]; iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]; ++iptr_ji) { |
69 |
if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) { |
70 |
d_ij=(-alpha)*MIN3(0.,fc->flux_matrix->mainBlock->val[iptr_ij],fc->flux_matrix->mainBlock->val[iptr_ji]); |
71 |
printf("%d %d -> %e\n",i,j,d_ij); |
72 |
fc->transport_matrix->mainBlock->val[iptr_ij]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ij]+d_ij; |
73 |
fc->transport_matrix->mainBlock->val[iptr_ji]+=alpha*fc->flux_matrix->mainBlock->val[iptr_ji]+d_ij; |
74 |
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]); |
75 |
fc->transport_matrix->mainBlock->val[fc->main_iptr[i]]-=d_ij; |
76 |
fc->transport_matrix->mainBlock->val[fc->main_iptr[j]]-=d_ij; |
77 |
break; |
78 |
} |
79 |
} |
80 |
} |
81 |
|
82 |
} |
83 |
/* TODO process couple block */ |
84 |
} |
85 |
} |
86 |
#pragma omp barrier |
87 |
} |
88 |
} |
89 |
|
90 |
} |
91 |
|
92 |
void Paso_FCTransportProblem_setFlux(Paso_FCTransportProblem * fc, double * u, double* fa) { |
93 |
|
94 |
double *remote_u=NULL; |
95 |
|
96 |
if (fc==NULL) return; |
97 |
Paso_SystemMatrix_startCollect(fc->transport_matrix,u); |
98 |
/* process main block */ |
99 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->mainBlock,u,0.,fa); |
100 |
/* finish exchange */ |
101 |
remote_u=Paso_SystemMatrix_finishCollect(fc->transport_matrix); |
102 |
/* process couple block */ |
103 |
Paso_SparseMatrix_MatrixVector_CSR_OFFSET0(1.,fc->transport_matrix->coupleBlock,remote_u,1.,fa); |
104 |
|
105 |
Paso_FCTransportProblem_setAntiDiffusiveFlux(fc,u,remote_u,fa); |
106 |
} |
107 |
/**************************************************************/ |
108 |
|
109 |
/* adds antidiffusion to fa |
110 |
|
111 |
P_p[i] + = sum_j min(0,a[i,j]) min(0,u[j]-u[i]) |
112 |
P_n[i] + = sum_j min(0,a[i,j]) max(0,u[j]-u[i]) |
113 |
Q_p[i] + = sum_j max(0,a[i,j]) max(0,u[j]-u[i]) |
114 |
Q_n[i] + = sum_j max(0,a[i,j]) min(0,u[j]-u[i]) |
115 |
d_ij=max(0,-a[i,j],-a[j,i]) |
116 |
l_ji=max(0,a[j,i],a[j,i]-a[i,j]) |
117 |
if a[j,i] >= a[i,j] and 0>a[i,j] : (i.e d_ij>0 and l_ji>=l_ij) |
118 |
r_ij = u[i]>u[j] ? Q_p[i]/P_p[i] : Q_n[i]/Q_n[i] |
119 |
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]) |
120 |
fa[i]+=f_ij |
121 |
fa[j]-=f_ij |
122 |
|
123 |
*/ |
124 |
|
125 |
void Paso_FCTransportProblem_setAntiDiffusiveFlux(Paso_FCTransportProblem * fc, double * u, double *u_remote, double* fa) { |
126 |
|
127 |
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; |
128 |
index_t color, iptr_ij,j,iptr_ji, i; |
129 |
dim_t n; |
130 |
|
131 |
|
132 |
if (fc==NULL) return; |
133 |
n=Paso_SystemMatrix_getTotalNumRows(fc->flux_matrix); |
134 |
/* exchange */ |
135 |
|
136 |
|
137 |
#pragma omp parallel private(color) |
138 |
{ |
139 |
for (color=0;color<fc->num_colors;++color) { |
140 |
#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) |
141 |
for (i = 0; i < n; ++i) { |
142 |
if (fc->colorOf[i]==color) { |
143 |
u_i=u[i]; |
144 |
/* gather the smoothness sensor */ |
145 |
P_p=0.; |
146 |
P_n=0.; |
147 |
Q_p=0.; |
148 |
Q_n=0.; |
149 |
#pragma ivdep |
150 |
for (iptr_ij=(fc->flux_matrix->mainBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->mainBlock->pattern->ptr[i+1]); ++iptr_ij) { |
151 |
a_ij=fc->flux_matrix->mainBlock->val[iptr_ij]; |
152 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
153 |
d=u[j]-u_i; |
154 |
if (a_ij<0.) { |
155 |
if (d<0.) { |
156 |
P_p+=a_ij*d; |
157 |
} else { |
158 |
P_n+=a_ij*d; |
159 |
} |
160 |
} else { |
161 |
if (d>0.) { |
162 |
Q_p+=a_ij*d; |
163 |
} else { |
164 |
Q_n+=a_ij*d; |
165 |
} |
166 |
} |
167 |
} |
168 |
#pragma ivdep |
169 |
for (iptr_ij=(fc->flux_matrix->coupleBlock->pattern->ptr[i]);iptr_ij<(fc->flux_matrix->coupleBlock->pattern->ptr[i+1]); ++iptr_ij) { |
170 |
a_ij=fc->flux_matrix->coupleBlock->val[iptr_ij]; |
171 |
j=fc->flux_matrix->coupleBlock->pattern->index[iptr_ij]; |
172 |
d=u_remote[j]-u_i; |
173 |
if (a_ij<0.) { |
174 |
if (d<0.) { |
175 |
P_p+=a_ij*d; |
176 |
} else { |
177 |
P_n+=a_ij*d; |
178 |
} |
179 |
} else { |
180 |
if (d>0.) { |
181 |
Q_p+=a_ij*d; |
182 |
} else { |
183 |
Q_n+=a_ij*d; |
184 |
} |
185 |
} |
186 |
} |
187 |
/* set the smoothness indicators */ |
188 |
r_p = (P_p > 0.) ? FLUX_LIMITER(Q_p/P_p) : 0.; |
189 |
r_n = (P_n < 0 ) ? FLUX_LIMITER(Q_n/P_n) : 0.; |
190 |
/* anti diffusive flux from main block */ |
191 |
for (iptr_ij=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ij<fc->flux_matrix->mainBlock->pattern->ptr[i+1]; ++iptr_ij) { |
192 |
a_ij=fc->flux_matrix->mainBlock->val[iptr_ij]; |
193 |
j=fc->flux_matrix->mainBlock->pattern->index[iptr_ij]; |
194 |
if (a_ij < 0 && i!=j) { |
195 |
/* find entry a[j,i] */ |
196 |
for (iptr_ji=fc->flux_matrix->mainBlock->pattern->ptr[i];iptr_ji<fc->flux_matrix->mainBlock->pattern->ptr[j+1]-1; ++iptr_ji) { |
197 |
if (fc->flux_matrix->mainBlock->pattern->index[iptr_ji]==i) { |
198 |
a_ji=fc->flux_matrix->mainBlock->val[iptr_ji]; |
199 |
if (a_ji > a_ij || (a_ji == a_ij && j<i) ) { |
200 |
u_j=u[j]; |
201 |
r_ij = u_i>u_j ? r_p : r_n; |
202 |
f_ij =MIN(r_ij*a_ij,a_ji-a_ij)*(u_i-u_j); |
203 |
fa[i]+=f_ij; |
204 |
fa[j]-=f_ij; |
205 |
break; |
206 |
} |
207 |
} |
208 |
} |
209 |
} |
210 |
} |
211 |
/* anti diffusive flux from couple block */ |
212 |
|
213 |
/* TODO */ |
214 |
} |
215 |
} |
216 |
#pragma omp barrier |
217 |
} |
218 |
} |
219 |
} |