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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1410 - (hide annotations)
Thu Feb 7 04:24:00 2008 UTC (11 years, 9 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 gross 1361 /* $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 gross 1407 /* 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 gross 1361
45 gross 1407 void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) {
46     dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i;
47 gross 1361 index_t color, iptr_ij,j,iptr_ji;
48 gross 1407 Paso_SystemMatrixPattern *pattern;
49     register double d_ij, sum, rtmp1, rtmp2;
50 gross 1361
51     if (fc==NULL) return;
52 gross 1407 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 gross 1361
59 gross 1407 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 gross 1401 }
81     }
82 gross 1407 /* TODO process couple block */
83    
84     /* set main diagonal entry */
85     fc->main_diagonal_low_order_transport_matrix[i]=sum;
86     }
87 gross 1361 }
88     }
89 gross 1407 /*
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 gross 1370
107 gross 1407 n=Paso_SystemMatrix_getTotalNumRows(L);
108 gross 1370
109 gross 1407 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 gross 1370
152 gross 1407 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 gross 1410 u_max_i=MAX(u_max_i,u_j);
167 gross 1407 }
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 gross 1410 u_max_i=MAX(u_max_i,u_j);
174 gross 1407 }
175     QN[i]=u_min_i-u_i;
176     QP[i]=u_max_i-u_i;
177     }
178 gross 1370 }
179 gross 1361
180 gross 1407 /*
181     *
182 gross 1410 * 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 gross 1407 * 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 gross 1361
371    
372 gross 1407 }
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 gross 1361
394 gross 1407 }
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 gross 1361
409 gross 1407 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 gross 1361
430 gross 1370
431 gross 1407 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 gross 1371
439 gross 1407 n=Paso_SystemMatrix_getTotalNumRows(f);
440     pattern=f->pattern;
441 gross 1408 #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
442 gross 1407 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 gross 1410 if (f_ij <=0) {
451 gross 1407 PN_i+=f_ij;
452     } else {
453     PP_i+=f_ij;
454     }
455 gross 1401 }
456 gross 1407 }
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 gross 1371
479 gross 1407 }
480 gross 1371
481 gross 1407 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 gross 1375
490 gross 1407 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 gross 1401 }
516 gross 1410 f[i]+=f_i;
517 gross 1407 }
518     Paso_SystemMatrix_startCollect(flux_matrix,RP);
519     remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix);
520 gross 1408 #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i)
521 gross 1407 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 gross 1361 }

  ViewVC Help
Powered by ViewVC 1.1.26