/[escript]/trunk/paso/src/SolverFCT_FluxControl.c
ViewVC logotype

Contents of /trunk/paso/src/SolverFCT_FluxControl.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1408 - (show annotations)
Mon Feb 4 07:19:50 2008 UTC (11 years, 9 months ago) by gross
File MIME type: text/plain
File size: 14864 byte(s)
some problems with openmp fixed.
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
30 /**************************************************************/
31
32 /* create the low order transport matrix and stores it negative values
33 * into the iteration_matrix accept the main diagonal which is stored seperately
34 * if fc->iteration_matrix==NULL, fc->iteration_matrix is allocated
35 *
36 * a=transport_matrix
37 * b= low_order_transport_matrix = - iteration_matrix
38 * c=main diagonal low_order_transport_matrix
39 * initialize c[i] mit a[i,i]
40 * d_ij=max(0,-a[i,j],-a[j,i])
41 * b[i,j]=-(a[i,j]+d_ij)
42 * c[i]-=d_ij
43 */
44
45 void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) {
46 dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;
47 index_t color, iptr_ij,j,iptr_ji;
48 Paso_SystemMatrixPattern *pattern;
49 register double d_ij, sum, rtmp1, rtmp2;
50
51 if (fc==NULL) return;
52 if (fc->iteration_matrix==NULL) {
53 fc->iteration_matrix=Paso_SystemMatrix_alloc(fc->transport_matrix->type,
54 fc->transport_matrix->pattern,
55 fc->transport_matrix->row_block_size,
56 fc->transport_matrix->col_block_size);
57 }
58
59 if (Paso_noError()) {
60 pattern=fc->iteration_matrix->pattern;
61 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
62 #pragma omp parallel for private(i,sum,iptr_ij,j,iptr_ji,rtmp1, rtmp2,d_ij) schedule(static)
63 for (i = 0; i < n; ++i) {
64 sum=fc->transport_matrix->mainBlock->val[fc->main_iptr[i]];
65 /* look at a[i,j] */
66 for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
67 j=pattern->mainPattern->index[iptr_ij];
68 if (j!=i) {
69 /* find entry a[j,i] */
70 for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
71 if (pattern->mainPattern->index[iptr_ji]==i) {
72 rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
73 rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];
74 d_ij=-MIN3(0.,rtmp1,rtmp2);
75 fc->iteration_matrix->mainBlock->val[iptr_ij]=-(rtmp1+d_ij);
76 sum-=d_ij;
77 break;
78 }
79 }
80 }
81 }
82 /* TODO process couple block */
83
84 /* set main diagonal entry */
85 fc->main_diagonal_low_order_transport_matrix[i]=sum;
86 }
87 }
88 }
89 /*
90 * out_i=m_i u_i + alpha \sum_{j <> i} l_{ij} (u_j-u_i) + beta q_i
91 *
92 */
93 void Paso_SolverFCT_setMuPaLuPbQ(double* out,
94 const double* M,
95 const double* u,
96 const double a,
97 Paso_SystemMatrix *L,
98 const double b,
99 const double* Q)
100 {
101 dim_t n, i, j;
102 Paso_SystemMatrixPattern *pattern;
103 double *remote_u;
104 register double sum, u_i, l_ij;
105 register index_t iptr_ij;
106
107 n=Paso_SystemMatrix_getTotalNumRows(L);
108
109 if (ABS(a)>0) {
110 Paso_SystemMatrix_startCollect(L,u);
111 }
112 #pragma omp parallel for private(i) schedule(static)
113 for (i = 0; i < n; ++i) {
114 out[i]=M[i]*u[i]+b*Q[i];
115 }
116 if (ABS(a)>0) {
117 pattern=L->pattern;
118 remote_u=Paso_SystemMatrix_finishCollect(L);
119 #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
120 for (i = 0; i < n; ++i) {
121 sum=0;
122 u_i=u[i];
123 #pragma ivdep
124 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
125 j=pattern->mainPattern->index[iptr_ij];
126 l_ij=L->mainBlock->val[iptr_ij];
127 sum+=l_ij*(u[j]-u_i);
128 }
129 #pragma ivdep
130 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
131 j=pattern->couplePattern->index[iptr_ij];
132 l_ij=L->coupleBlock->val[iptr_ij];
133 sum+=l_ij*(remote_u[j]-u_i);
134 }
135 out[i]+=a*sum;
136 }
137 }
138 }
139 /*
140 * calculate QP[i] max_{i in L->pattern[i]} (u[j]-u[i] )
141 * QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
142 *
143 */
144 void Paso_SolverFCT_setQs(const double* u,double* QN, double* QP, Paso_SystemMatrix *L)
145 {
146 dim_t n, i, j;
147 Paso_SystemMatrixPattern *pattern;
148 double *remote_u;
149 register double u_i, u_min_i, u_max_i, u_j;
150 register index_t iptr_ij;
151
152 Paso_SystemMatrix_startCollect(L,u);
153 n=Paso_SystemMatrix_getTotalNumRows(L);
154 pattern=L->pattern;
155 remote_u=Paso_SystemMatrix_finishCollect(L);
156 #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)
157 for (i = 0; i < n; ++i) {
158 u_i=u[i];
159 u_min_i=u_i;
160 u_max_i=u_i;
161 #pragma ivdep
162 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
163 j=pattern->mainPattern->index[iptr_ij];
164 u_j=u[j];
165 u_min_i=MIN(u_min_i,u_j);
166 u_max_i=MIN(u_max_i,u_j);
167 }
168 #pragma ivdep
169 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
170 j=pattern->couplePattern->index[iptr_ij];
171 u_j=remote_u[j];
172 u_min_i=MIN(u_min_i,u_j);
173 u_max_i=MIN(u_max_i,u_j);
174 }
175 QN[i]=u_min_i-u_i;
176 QP[i]=u_max_i-u_i;
177 }
178 }
179
180 /*
181 *
182 * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])
183 *
184 * m=fc->mass matrix
185 * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
186 */
187 void Paso_FCTransportProblem_updateAntiDiffusionFlux(const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix,const double a, const double b, const double* u)
188 {
189 dim_t n, j, i;
190 double *remote_u;
191 index_t iptr_ij;
192 const double b2=-b;
193 register double m_ij, ml_ij, k_ij, u_i;
194 Paso_SystemMatrixPattern *pattern;
195 Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);
196 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
197 pattern=fc->iteration_matrix->pattern;
198 remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
199
200 if ( (ABS(a) >0 ) ) {
201 if ( (ABS(b) >0 ) ) {
202 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)
203 for (i = 0; i < n; ++i) {
204 u_i=u[i];
205 #pragma ivdep
206 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
207 j=pattern->mainPattern->index[iptr_ij];
208 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
209 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
210 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
211 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
212 }
213 #pragma ivdep
214 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
215 j=pattern->couplePattern->index[iptr_ij];
216 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
217 k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
218 ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
219 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
220 }
221 }
222 } else {
223 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij)
224 for (i = 0; i < n; ++i) {
225 u_i=u[i];
226 #pragma ivdep
227 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
228 j=pattern->mainPattern->index[iptr_ij];
229 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
230 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*a*m_ij;
231 }
232 #pragma ivdep
233 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
234 j=pattern->couplePattern->index[iptr_ij];
235 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
236 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij;
237 }
238 }
239
240
241 }
242 } else {
243 if ( (ABS(b) >0 ) ) {
244 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,k_ij, ml_ij)
245 for (i = 0; i < n; ++i) {
246 u_i=u[i];
247 #pragma ivdep
248 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
249 j=pattern->mainPattern->index[iptr_ij];
250 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
251 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
252 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij);
253 }
254 #pragma ivdep
255 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
256 j=pattern->couplePattern->index[iptr_ij];
257 k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
258 ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
259 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*b2*(ml_ij+k_ij);
260 }
261 }
262
263 }
264 }
265 }
266 /*
267 * apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
268 *
269 */
270 void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)
271 {
272 dim_t n, i, j;
273 Paso_SystemMatrixPattern *pattern;
274 double *remote_u;
275 register double u_i, f_ij;
276 register index_t iptr_ij;
277
278 n=Paso_SystemMatrix_getTotalNumRows(f);
279 pattern=f->pattern;
280 remote_u=Paso_SystemMatrix_finishCollect(f);
281 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
282 for (i = 0; i < n; ++i) {
283 u_i=u[i];
284 #pragma ivdep
285 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
286 j=pattern->mainPattern->index[iptr_ij];
287 f_ij=f->mainBlock->val[iptr_ij];
288 if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
289 }
290 #pragma ivdep
291 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
292 j=pattern->couplePattern->index[iptr_ij];
293 f_ij=f->coupleBlock->val[iptr_ij];
294 if (f_ij * (u_i-remote_u[j]) <= 0) f->coupleBlock->val[iptr_ij]=0;
295 }
296 }
297 }
298
299
300 void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
301 const double* QN,const double* QP,double* RN,double* RP)
302 {
303 dim_t n, i, j;
304 Paso_SystemMatrixPattern *pattern;
305 register double f_ij, PP_i, PN_i;
306 register index_t iptr_ij;
307
308 n=Paso_SystemMatrix_getTotalNumRows(f);
309 pattern=f->pattern;
310 #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
311 for (i = 0; i < n; ++i) {
312 PP_i=0.;
313 PN_i=0.;
314 #pragma ivdep
315 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
316 j=pattern->mainPattern->index[iptr_ij];
317 if (i != j ) {
318 f_ij=f->mainBlock->val[iptr_ij];
319 if (f_ij <=0 ) {
320 PN_i+=f_ij;
321 } else {
322 PP_i+=f_ij;
323 }
324 }
325 }
326 #pragma ivdep
327 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
328 j=pattern->couplePattern->index[iptr_ij];
329 f_ij=f->coupleBlock->val[iptr_ij];
330 if (f_ij <=0 ) {
331 PN_i+=f_ij;
332 } else {
333 PP_i+=f_ij;
334 }
335 }
336 if (PN_i<0) {
337 RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
338 } else {
339 RN[i]=1.;
340 }
341 if (PP_i>0) {
342 RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
343 } else {
344 RP[i]=1.;
345 }
346 }
347
348 }
349
350 void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP)
351 {
352 dim_t n, i, j;
353 Paso_SystemMatrixPattern *pattern;
354 double *remote_RN, *remote_RP;
355 register double RN_i, RP_i, f_i, f_ij;
356 register index_t iptr_ij;
357 n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
358
359 Paso_SystemMatrix_startCollect(flux_matrix,RN);
360 pattern=flux_matrix->pattern;
361 remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);
362 #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
363 for (i = 0; i < n; ++i) {
364 RN_i=RN[i];
365 RP_i=RP[i];
366 f_i=0;
367 #pragma ivdep
368 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
369 j=pattern->mainPattern->index[iptr_ij];
370 f_ij=flux_matrix->mainBlock->val[iptr_ij];
371 if (f_ij >=0) {
372 f_i+=f_ij*MIN(RP_i,RN[j]);
373 } else {
374 f_i+=f_ij*MIN(RN_i,RP[j]);
375 }
376 }
377 #pragma ivdep
378 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
379 j=pattern->couplePattern->index[iptr_ij];
380 f_ij=flux_matrix->coupleBlock->val[iptr_ij];
381 if (f_ij >=0) {
382 f_i+=f_ij*MIN(RP_i,remote_RN[j]);
383 }
384 }
385 f[i]=f_i;
386 }
387 Paso_SystemMatrix_startCollect(flux_matrix,RP);
388 remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);
389 #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i)
390 for (i = 0; i < n; ++i) {
391 RN_i=RN[i];
392 f_i=0;
393 #pragma ivdep
394 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
395 j=pattern->couplePattern->index[iptr_ij];
396 f_ij=flux_matrix->coupleBlock->val[iptr_ij];
397 if (f_ij < 0) {
398 f_i+=f_ij*MIN(RN_i,remote_RP[j]);
399 }
400 }
401 f[i]+=f_i;
402 }
403 }

  ViewVC Help
Powered by ViewVC 1.1.26