/[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 1811 - (hide annotations)
Thu Sep 25 23:11:13 2008 UTC (11 years, 4 months ago) by ksteube
File MIME type: text/plain
File size: 18679 byte(s)
Copyright updated in all files

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

  ViewVC Help
Powered by ViewVC 1.1.26