/[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 1410 - (show annotations)
Thu Feb 7 04:24:00 2008 UTC (11 years, 8 months ago) by gross
File MIME type: text/plain
File size: 21643 byte(s)
a new version of the algebric upwinding a. flux limiter
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=MAX(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=MAX(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} = (m_{ij} - dt (1-theta) d_{ij}) (u_last[j]-u_last[i]) - (m_{ij} + dt theta 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_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const double* u, const double* u_last)
188 {
189 dim_t n, j, i;
190 double *remote_u=NULL, *remote_u_last=NULL;
191 index_t iptr_ij;
192 const double f1=- dt * (1.-fc->theta);
193 const double f2= dt * fc->theta;
194 register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
195 Paso_SystemMatrixPattern *pattern;
196 Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);
197 Paso_SystemMatrix_startCollect(fc->iteration_matrix,u_last);
198 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
199 pattern=fc->iteration_matrix->pattern;
200 remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
201 remote_u_last=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
202 if ( (ABS(f1) >0 ) ) {
203 if ( (ABS(f2) >0 ) ) {
204 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
205 for (i = 0; i < n; ++i) {
206 u_i=u[i];
207 u_last_i=u_last[i];
208 #pragma ivdep
209 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
210 j=pattern->mainPattern->index[iptr_ij];
211 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
212 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
213 fc->iteration_matrix->mainBlock->val[iptr_ij]);
214 d_u=u[j]-u_i;
215 d_u_last=u_last[j]-u_last_i;
216 flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
217 }
218 #pragma ivdep
219 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
220 j=pattern->couplePattern->index[iptr_ij];
221 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
222 d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
223 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
224 d_u=remote_u[j]-u_i;
225 d_u_last=remote_u_last[j]-u_last_i;
226 flux_matrix->coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
227 }
228 }
229 } else {
230 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
231 for (i = 0; i < n; ++i) {
232 u_i=u[i];
233 u_last_i=u_last[i];
234 #pragma ivdep
235 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
236 j=pattern->mainPattern->index[iptr_ij];
237 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
238 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
239 fc->iteration_matrix->mainBlock->val[iptr_ij]);
240 d_u=u[j]-u_i;
241 d_u_last=u_last[j]-u_last_i;
242 flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
243 }
244 #pragma ivdep
245 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
246 j=pattern->couplePattern->index[iptr_ij];
247 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
248 d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
249 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
250 d_u=remote_u[j]-u_i;
251 d_u_last=remote_u_last[j]-u_last_i;
252 flux_matrix->coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
253 }
254 }
255 }
256 } else {
257 if ( (ABS(f2) >0 ) ) {
258 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
259 for (i = 0; i < n; ++i) {
260 u_i=u[i];
261 u_last_i=u_last[i];
262 #pragma ivdep
263 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
264 j=pattern->mainPattern->index[iptr_ij];
265 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
266 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
267 fc->iteration_matrix->mainBlock->val[iptr_ij]);
268 d_u=u[j]-u_i;
269 d_u_last=u_last[j]-u_last_i;
270 flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
271 }
272 #pragma ivdep
273 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
274 j=pattern->couplePattern->index[iptr_ij];
275 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
276 d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
277 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
278 d_u=remote_u[j]-u_i;
279 d_u_last=remote_u_last[j]-u_last_i;
280 flux_matrix->coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
281 }
282 }
283 } else {
284 #pragma omp parallel for schedule(static) private(i, u_i, u_last_i, iptr_ij, j,m_ij,d_ij, d_u_last, d_u)
285 for (i = 0; i < n; ++i) {
286 u_i=u[i];
287 u_last_i=u_last[i];
288 #pragma ivdep
289 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
290 j=pattern->mainPattern->index[iptr_ij];
291 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
292 d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
293 fc->iteration_matrix->mainBlock->val[iptr_ij]);
294 d_u=u[j]-u_i;
295 d_u_last=u_last[j]-u_last_i;
296 flux_matrix->mainBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
297 }
298 #pragma ivdep
299 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
300 j=pattern->couplePattern->index[iptr_ij];
301 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
302 d_ij=-(fc->transport_matrix->coupleBlock->val[iptr_ij]+
303 fc->iteration_matrix->coupleBlock->val[iptr_ij]);
304 d_u=remote_u[j]-u_i;
305 d_u_last=remote_u_last[j]-u_last_i;
306 flux_matrix->coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
307 }
308 }
309 }
310 }
311 }
312 /*
313 *
314 * f_{ij} + = (a*m_{ij} + b* d_{ij}) (u[j]-u[i])
315 *
316 * m=fc->mass matrix
317 * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
318 */
319 void Paso_FCTransportProblem_updateAntiDiffusionFlux(const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix,const double a, const double b, const double* u)
320 {
321 dim_t n, j, i;
322 double *remote_u;
323 index_t iptr_ij;
324 const double b2=-b;
325 register double m_ij, ml_ij, k_ij, u_i;
326 Paso_SystemMatrixPattern *pattern;
327 Paso_SystemMatrix_startCollect(fc->iteration_matrix,u);
328 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
329 pattern=fc->iteration_matrix->pattern;
330 remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix);
331 if ( (ABS(a) >0 ) ) {
332 if ( (ABS(b) >0 ) ) {
333 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij)
334 for (i = 0; i < n; ++i) {
335 u_i=u[i];
336 #pragma ivdep
337 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
338 j=pattern->mainPattern->index[iptr_ij];
339 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
340 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
341 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
342 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
343 }
344 #pragma ivdep
345 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
346 j=pattern->couplePattern->index[iptr_ij];
347 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
348 k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
349 ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
350 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij));
351 }
352 }
353 } else {
354 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij)
355 for (i = 0; i < n; ++i) {
356 u_i=u[i];
357 #pragma ivdep
358 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
359 j=pattern->mainPattern->index[iptr_ij];
360 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
361 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*a*m_ij;
362 }
363 #pragma ivdep
364 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
365 j=pattern->couplePattern->index[iptr_ij];
366 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij];
367 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij;
368 }
369 }
370
371
372 }
373 } else {
374 if ( (ABS(b) >0 ) ) {
375 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,k_ij, ml_ij)
376 for (i = 0; i < n; ++i) {
377 u_i=u[i];
378 #pragma ivdep
379 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
380 j=pattern->mainPattern->index[iptr_ij];
381 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij];
382 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij];
383 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij);
384 }
385 #pragma ivdep
386 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
387 j=pattern->couplePattern->index[iptr_ij];
388 k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij];
389 ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij];
390 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*b2*(ml_ij+k_ij);
391 }
392 }
393
394 }
395 }
396 }
397 /*
398 * apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
399 *
400 */
401 void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u)
402 {
403 dim_t n, i, j;
404 Paso_SystemMatrixPattern *pattern;
405 double *remote_u;
406 register double u_i, f_ij;
407 register index_t iptr_ij;
408
409 n=Paso_SystemMatrix_getTotalNumRows(f);
410 pattern=f->pattern;
411 remote_u=Paso_SystemMatrix_finishCollect(f);
412 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
413 for (i = 0; i < n; ++i) {
414 u_i=u[i];
415 #pragma ivdep
416 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
417 j=pattern->mainPattern->index[iptr_ij];
418 f_ij=f->mainBlock->val[iptr_ij];
419 if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
420 }
421 #pragma ivdep
422 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
423 j=pattern->couplePattern->index[iptr_ij];
424 f_ij=f->coupleBlock->val[iptr_ij];
425 if (f_ij * (u_i-remote_u[j]) <= 0) f->coupleBlock->val[iptr_ij]=0;
426 }
427 }
428 }
429
430
431 void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
432 const double* QN,const double* QP,double* RN,double* RP)
433 {
434 dim_t n, i, j;
435 Paso_SystemMatrixPattern *pattern;
436 register double f_ij, PP_i, PN_i;
437 register index_t iptr_ij;
438
439 n=Paso_SystemMatrix_getTotalNumRows(f);
440 pattern=f->pattern;
441 #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
442 for (i = 0; i < n; ++i) {
443 PP_i=0.;
444 PN_i=0.;
445 #pragma ivdep
446 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
447 j=pattern->mainPattern->index[iptr_ij];
448 if (i != j ) {
449 f_ij=f->mainBlock->val[iptr_ij];
450 if (f_ij <=0) {
451 PN_i+=f_ij;
452 } else {
453 PP_i+=f_ij;
454 }
455 }
456 }
457 #pragma ivdep
458 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
459 j=pattern->couplePattern->index[iptr_ij];
460 f_ij=f->coupleBlock->val[iptr_ij];
461 if (f_ij <=0 ) {
462 PN_i+=f_ij;
463 } else {
464 PP_i+=f_ij;
465 }
466 }
467 if (PN_i<0) {
468 RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
469 } else {
470 RN[i]=1.;
471 }
472 if (PP_i>0) {
473 RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
474 } else {
475 RP[i]=1.;
476 }
477 }
478
479 }
480
481 void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP)
482 {
483 dim_t n, i, j;
484 Paso_SystemMatrixPattern *pattern;
485 double *remote_RN, *remote_RP;
486 register double RN_i, RP_i, f_i, f_ij;
487 register index_t iptr_ij;
488 n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
489
490 Paso_SystemMatrix_startCollect(flux_matrix,RN);
491 pattern=flux_matrix->pattern;
492 remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix);
493 #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
494 for (i = 0; i < n; ++i) {
495 RN_i=RN[i];
496 RP_i=RP[i];
497 f_i=0;
498 #pragma ivdep
499 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
500 j=pattern->mainPattern->index[iptr_ij];
501 f_ij=flux_matrix->mainBlock->val[iptr_ij];
502 if (f_ij >=0) {
503 f_i+=f_ij*MIN(RP_i,RN[j]);
504 } else {
505 f_i+=f_ij*MIN(RN_i,RP[j]);
506 }
507 }
508 #pragma ivdep
509 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
510 j=pattern->couplePattern->index[iptr_ij];
511 f_ij=flux_matrix->coupleBlock->val[iptr_ij];
512 if (f_ij >=0) {
513 f_i+=f_ij*MIN(RP_i,remote_RN[j]);
514 }
515 }
516 f[i]+=f_i;
517 }
518 Paso_SystemMatrix_startCollect(flux_matrix,RP);
519 remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);
520 #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i)
521 for (i = 0; i < n; ++i) {
522 RN_i=RN[i];
523 f_i=0;
524 #pragma ivdep
525 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ij<pattern->couplePattern->ptr[i+1]; ++iptr_ij) {
526 j=pattern->couplePattern->index[iptr_ij];
527 f_ij=flux_matrix->coupleBlock->val[iptr_ij];
528 if (f_ij < 0) {
529 f_i+=f_ij*MIN(RN_i,remote_RP[j]);
530 }
531 }
532 f[i]+=f_i;
533 }
534 }

  ViewVC Help
Powered by ViewVC 1.1.26