/[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 1562 - (hide annotations)
Wed May 21 13:04:40 2008 UTC (11 years, 9 months ago) by gross
File MIME type: text/plain
File size: 18673 byte(s)
The algebraic upwinding with MPI. The case of boundary constraint needs still some attention. 


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 gross 1562 #pragma ivdep
67 gross 1407 for (iptr_ij=pattern->mainPattern->ptr[i];iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
68     j=pattern->mainPattern->index[iptr_ij];
69     if (j!=i) {
70     /* find entry a[j,i] */
71     for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_ji<pattern->mainPattern->ptr[j+1]; ++iptr_ji) {
72     if (pattern->mainPattern->index[iptr_ji]==i) {
73     rtmp1=fc->transport_matrix->mainBlock->val[iptr_ij];
74     rtmp2=fc->transport_matrix->mainBlock->val[iptr_ji];
75     d_ij=-MIN3(0.,rtmp1,rtmp2);
76     fc->iteration_matrix->mainBlock->val[iptr_ij]=-(rtmp1+d_ij);
77     sum-=d_ij;
78     break;
79     }
80     }
81 gross 1401 }
82     }
83 gross 1562 #pragma ivdep
84     for (iptr_ij=pattern->col_couplePattern->ptr[i];iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
85     j=pattern->col_couplePattern->index[iptr_ij];
86     /* find entry a[j,i] */
87     for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_ji<pattern->row_couplePattern->ptr[j+1]; ++iptr_ji) {
88     if (pattern->row_couplePattern->index[iptr_ji]==i) {
89     rtmp1=fc->transport_matrix->col_coupleBlock->val[iptr_ij];
90     rtmp2=fc->transport_matrix->row_coupleBlock->val[iptr_ji];
91     d_ij=-MIN3(0.,rtmp1,rtmp2);
92     fc->iteration_matrix->col_coupleBlock->val[iptr_ij]=-(rtmp1+d_ij);
93     fc->iteration_matrix->row_coupleBlock->val[iptr_ji]=-(rtmp2+d_ij);
94     sum-=d_ij;
95     break;
96     }
97     }
98     }
99 gross 1407 /* set main diagonal entry */
100     fc->main_diagonal_low_order_transport_matrix[i]=sum;
101     }
102 gross 1361 }
103     }
104 gross 1407 /*
105     * out_i=m_i u_i + alpha \sum_{j <> i} l_{ij} (u_j-u_i) + beta q_i
106     *
107     */
108     void Paso_SolverFCT_setMuPaLuPbQ(double* out,
109     const double* M,
110 gross 1562 const Paso_Coupler* u_coupler,
111 gross 1407 const double a,
112 gross 1562 const Paso_SystemMatrix *L,
113 gross 1407 const double b,
114     const double* Q)
115     {
116     dim_t n, i, j;
117     Paso_SystemMatrixPattern *pattern;
118 gross 1562 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
119     const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
120 gross 1407 register double sum, u_i, l_ij;
121     register index_t iptr_ij;
122 gross 1370
123 gross 1407 n=Paso_SystemMatrix_getTotalNumRows(L);
124 gross 1370
125 gross 1407 #pragma omp parallel for private(i) schedule(static)
126     for (i = 0; i < n; ++i) {
127     out[i]=M[i]*u[i]+b*Q[i];
128     }
129     if (ABS(a)>0) {
130     pattern=L->pattern;
131     #pragma omp parallel for schedule(static) private(i, sum, u_i, iptr_ij, j, l_ij)
132     for (i = 0; i < n; ++i) {
133     sum=0;
134     u_i=u[i];
135     #pragma ivdep
136     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
137     j=pattern->mainPattern->index[iptr_ij];
138     l_ij=L->mainBlock->val[iptr_ij];
139     sum+=l_ij*(u[j]-u_i);
140     }
141     #pragma ivdep
142 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
143     j=pattern->col_couplePattern->index[iptr_ij];
144     l_ij=L->col_coupleBlock->val[iptr_ij];
145 gross 1407 sum+=l_ij*(remote_u[j]-u_i);
146     }
147     out[i]+=a*sum;
148     }
149     }
150     }
151     /*
152     * calculate QP[i] max_{i in L->pattern[i]} (u[j]-u[i] )
153     * QN[i] min_{i in L->pattern[i]} (u[j]-u[i] )
154     *
155     */
156 gross 1562 void Paso_SolverFCT_setQs(const Paso_Coupler* u_coupler,double* QN, double* QP, const Paso_SystemMatrix *L)
157 gross 1407 {
158     dim_t n, i, j;
159     Paso_SystemMatrixPattern *pattern;
160 gross 1562 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
161     const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
162 gross 1407 register double u_i, u_min_i, u_max_i, u_j;
163     register index_t iptr_ij;
164     n=Paso_SystemMatrix_getTotalNumRows(L);
165     pattern=L->pattern;
166     #pragma omp parallel for schedule(static) private(i, u_i, u_min_i, u_max_i, j, u_j, iptr_ij)
167     for (i = 0; i < n; ++i) {
168     u_i=u[i];
169     u_min_i=u_i;
170     u_max_i=u_i;
171     #pragma ivdep
172     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
173     j=pattern->mainPattern->index[iptr_ij];
174     u_j=u[j];
175     u_min_i=MIN(u_min_i,u_j);
176 gross 1410 u_max_i=MAX(u_max_i,u_j);
177 gross 1407 }
178     #pragma ivdep
179 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
180     j=pattern->col_couplePattern->index[iptr_ij];
181 gross 1407 u_j=remote_u[j];
182     u_min_i=MIN(u_min_i,u_j);
183 gross 1410 u_max_i=MAX(u_max_i,u_j);
184 gross 1407 }
185     QN[i]=u_min_i-u_i;
186     QP[i]=u_max_i-u_i;
187     }
188 gross 1370 }
189 gross 1361
190 gross 1407 /*
191     *
192 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])
193     *
194     * m=fc->mass matrix
195     * d=artifical diffusion matrix = L - K = - fc->iteration matrix - fc->transport matrix (away from main diagonal)
196     */
197 gross 1562 void Paso_FCTransportProblem_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const Paso_Coupler* u_coupler)
198 gross 1410 {
199     dim_t n, j, i;
200     index_t iptr_ij;
201 gross 1562 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
202     const double *u_last= Paso_Coupler_borrowLocalData(fc->u_coupler);
203     const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
204     const double *remote_u_last=Paso_Coupler_borrowRemoteData(fc->u_coupler);
205 gross 1410 const double f1=- dt * (1.-fc->theta);
206     const double f2= dt * fc->theta;
207     register double m_ij, d_ij, u_i, u_last_i, d_u_last, d_u;
208 gross 1562 const Paso_SystemMatrixPattern *pattern=fc->iteration_matrix->pattern;
209 gross 1410 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix);
210     if ( (ABS(f1) >0 ) ) {
211     if ( (ABS(f2) >0 ) ) {
212     #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)
213     for (i = 0; i < n; ++i) {
214     u_i=u[i];
215     u_last_i=u_last[i];
216     #pragma ivdep
217     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
218     j=pattern->mainPattern->index[iptr_ij];
219     m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
220     d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
221     fc->iteration_matrix->mainBlock->val[iptr_ij]);
222     d_u=u[j]-u_i;
223     d_u_last=u_last[j]-u_last_i;
224     flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
225     }
226     #pragma ivdep
227 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
228     j=pattern->col_couplePattern->index[iptr_ij];
229     m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
230     d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
231     fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
232 gross 1410 d_u=remote_u[j]-u_i;
233     d_u_last=remote_u_last[j]-u_last_i;
234 gross 1552 flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u;
235 gross 1410 }
236     }
237     } else {
238     #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)
239     for (i = 0; i < n; ++i) {
240     u_i=u[i];
241     u_last_i=u_last[i];
242     #pragma ivdep
243     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
244     j=pattern->mainPattern->index[iptr_ij];
245     m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
246     d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
247     fc->iteration_matrix->mainBlock->val[iptr_ij]);
248     d_u=u[j]-u_i;
249     d_u_last=u_last[j]-u_last_i;
250     flux_matrix->mainBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
251     }
252     #pragma ivdep
253 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
254     j=pattern->col_couplePattern->index[iptr_ij];
255     m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
256     d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
257     fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
258 gross 1410 d_u=remote_u[j]-u_i;
259     d_u_last=remote_u_last[j]-u_last_i;
260 gross 1552 flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u;
261 gross 1410 }
262     }
263     }
264     } else {
265     if ( (ABS(f2) >0 ) ) {
266     #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)
267     for (i = 0; i < n; ++i) {
268     u_i=u[i];
269     u_last_i=u_last[i];
270     #pragma ivdep
271     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
272     j=pattern->mainPattern->index[iptr_ij];
273     m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
274     d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
275     fc->iteration_matrix->mainBlock->val[iptr_ij]);
276     d_u=u[j]-u_i;
277     d_u_last=u_last[j]-u_last_i;
278     flux_matrix->mainBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
279     }
280     #pragma ivdep
281 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
282     j=pattern->col_couplePattern->index[iptr_ij];
283     m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
284     d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
285     fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
286 gross 1410 d_u=remote_u[j]-u_i;
287     d_u_last=remote_u_last[j]-u_last_i;
288 gross 1552 flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u;
289 gross 1410 }
290     }
291     } else {
292     #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)
293     for (i = 0; i < n; ++i) {
294     u_i=u[i];
295     u_last_i=u_last[i];
296     #pragma ivdep
297     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
298     j=pattern->mainPattern->index[iptr_ij];
299     m_ij=fc->mass_matrix->mainBlock->val[iptr_ij];
300     d_ij=-(fc->transport_matrix->mainBlock->val[iptr_ij]+
301     fc->iteration_matrix->mainBlock->val[iptr_ij]);
302     d_u=u[j]-u_i;
303     d_u_last=u_last[j]-u_last_i;
304     flux_matrix->mainBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
305     }
306     #pragma ivdep
307 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
308     j=pattern->col_couplePattern->index[iptr_ij];
309     m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij];
310     d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+
311     fc->iteration_matrix->col_coupleBlock->val[iptr_ij]);
312 gross 1410 d_u=remote_u[j]-u_i;
313     d_u_last=remote_u_last[j]-u_last_i;
314 gross 1552 flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u);
315 gross 1410 }
316     }
317     }
318     }
319     }
320     /*
321 gross 1407 * apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0
322     *
323     */
324 gross 1562 void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler)
325 gross 1407 {
326     dim_t n, i, j;
327     Paso_SystemMatrixPattern *pattern;
328 gross 1562 const double *u=Paso_Coupler_borrowLocalData(u_coupler);
329     const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler);
330 gross 1407 register double u_i, f_ij;
331     register index_t iptr_ij;
332 gross 1361
333 gross 1407 n=Paso_SystemMatrix_getTotalNumRows(f);
334     pattern=f->pattern;
335     #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij)
336     for (i = 0; i < n; ++i) {
337     u_i=u[i];
338     #pragma ivdep
339     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
340     j=pattern->mainPattern->index[iptr_ij];
341     f_ij=f->mainBlock->val[iptr_ij];
342     if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0;
343     }
344     #pragma ivdep
345 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
346     j=pattern->col_couplePattern->index[iptr_ij];
347     f_ij=f->col_coupleBlock->val[iptr_ij];
348     if (f_ij * (u_i-remote_u[j]) <= 0) f->col_coupleBlock->val[iptr_ij]=0;
349 gross 1407 }
350     }
351     }
352 gross 1361
353 gross 1370
354 gross 1407 void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix,
355 gross 1562 const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP)
356 gross 1407 {
357     dim_t n, i, j;
358     Paso_SystemMatrixPattern *pattern;
359 gross 1562 const double *QN=Paso_Coupler_borrowLocalData(QN_coupler);
360     const double *QP=Paso_Coupler_borrowLocalData(QP_coupler);
361 gross 1407 register double f_ij, PP_i, PN_i;
362     register index_t iptr_ij;
363 gross 1371
364 gross 1407 n=Paso_SystemMatrix_getTotalNumRows(f);
365     pattern=f->pattern;
366 gross 1408 #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i)
367 gross 1407 for (i = 0; i < n; ++i) {
368     PP_i=0.;
369     PN_i=0.;
370     #pragma ivdep
371     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
372     j=pattern->mainPattern->index[iptr_ij];
373     if (i != j ) {
374     f_ij=f->mainBlock->val[iptr_ij];
375 gross 1410 if (f_ij <=0) {
376 gross 1407 PN_i+=f_ij;
377     } else {
378     PP_i+=f_ij;
379     }
380 gross 1401 }
381 gross 1407 }
382     #pragma ivdep
383 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
384     j=pattern->col_couplePattern->index[iptr_ij];
385     f_ij=f->col_coupleBlock->val[iptr_ij];
386 gross 1407 if (f_ij <=0 ) {
387     PN_i+=f_ij;
388     } else {
389     PP_i+=f_ij;
390     }
391     }
392     if (PN_i<0) {
393     RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i);
394     } else {
395     RN[i]=1.;
396     }
397     if (PP_i>0) {
398     RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i);
399     } else {
400     RP[i]=1.;
401     }
402     }
403 gross 1371
404 gross 1407 }
405 gross 1371
406 gross 1562 void Paso_FCTransportProblem_addCorrectedFluxes(double* f,const Paso_SystemMatrix *flux_matrix, const Paso_Coupler* RN_coupler, const Paso_Coupler* RP_coupler)
407 gross 1407 {
408     dim_t n, i, j;
409     Paso_SystemMatrixPattern *pattern;
410     register double RN_i, RP_i, f_i, f_ij;
411     register index_t iptr_ij;
412 gross 1562 const double *RN=Paso_Coupler_borrowLocalData(RN_coupler);
413     const double *remote_RN=Paso_Coupler_borrowRemoteData(RN_coupler);
414     const double *RP=Paso_Coupler_borrowLocalData(RP_coupler);
415     const double *remote_RP=Paso_Coupler_borrowRemoteData(RP_coupler);
416 gross 1407 n=Paso_SystemMatrix_getTotalNumRows(flux_matrix);
417 gross 1375
418 gross 1407 pattern=flux_matrix->pattern;
419     #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i)
420     for (i = 0; i < n; ++i) {
421     RN_i=RN[i];
422     RP_i=RP[i];
423     f_i=0;
424     #pragma ivdep
425     for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ij<pattern->mainPattern->ptr[i+1]; ++iptr_ij) {
426     j=pattern->mainPattern->index[iptr_ij];
427     f_ij=flux_matrix->mainBlock->val[iptr_ij];
428     if (f_ij >=0) {
429     f_i+=f_ij*MIN(RP_i,RN[j]);
430     } else {
431     f_i+=f_ij*MIN(RN_i,RP[j]);
432     }
433     }
434     #pragma ivdep
435 gross 1552 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ij<pattern->col_couplePattern->ptr[i+1]; ++iptr_ij) {
436     j=pattern->col_couplePattern->index[iptr_ij];
437     f_ij=flux_matrix->col_coupleBlock->val[iptr_ij];
438 gross 1407 if (f_ij >=0) {
439     f_i+=f_ij*MIN(RP_i,remote_RN[j]);
440 gross 1562 }else {
441     f_i+=f_ij*MIN(RN_i,remote_RP[j]);
442 gross 1407 }
443 gross 1401 }
444 gross 1410 f[i]+=f_i;
445 gross 1407 }
446 gross 1361 }

  ViewVC Help
Powered by ViewVC 1.1.26