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

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_ijmainPattern->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_jimainPattern->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_ijcol_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_jirow_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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_ijmainPattern->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_ijcol_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 }