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

Revision 1552 - (show annotations)
Thu May 8 08:52:41 2008 UTC (11 years, 9 months ago) by gross
File MIME type: text/plain
File size: 21923 byte(s)
some changes to make the implementatiopn of a upwind MPI version easier

 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_ijmainPattern->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_jimainPattern->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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 131 j=pattern->col_couplePattern->index[iptr_ij]; 132 l_ij=L->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 170 j=pattern->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 220 j=pattern->col_couplePattern->index[iptr_ij]; 221 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij]; 222 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+ 223 fc->iteration_matrix->col_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->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 246 j=pattern->col_couplePattern->index[iptr_ij]; 247 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij]; 248 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+ 249 fc->iteration_matrix->col_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->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 274 j=pattern->col_couplePattern->index[iptr_ij]; 275 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij]; 276 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+ 277 fc->iteration_matrix->col_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->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 300 j=pattern->col_couplePattern->index[iptr_ij]; 301 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij]; 302 d_ij=-(fc->transport_matrix->col_coupleBlock->val[iptr_ij]+ 303 fc->iteration_matrix->col_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->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 346 j=pattern->col_couplePattern->index[iptr_ij]; 347 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij]; 348 k_ij=fc->transport_matrix->col_coupleBlock->val[iptr_ij]; 349 ml_ij=fc->iteration_matrix->col_coupleBlock->val[iptr_ij]; 350 flux_matrix->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 365 j=pattern->col_couplePattern->index[iptr_ij]; 366 m_ij=fc->mass_matrix->col_coupleBlock->val[iptr_ij]; 367 flux_matrix->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 387 j=pattern->col_couplePattern->index[iptr_ij]; 388 k_ij=fc->transport_matrix->col_coupleBlock->val[iptr_ij]; 389 ml_ij=fc->iteration_matrix->col_coupleBlock->val[iptr_ij]; 390 flux_matrix->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 423 j=pattern->col_couplePattern->index[iptr_ij]; 424 f_ij=f->col_coupleBlock->val[iptr_ij]; 425 if (f_ij * (u_i-remote_u[j]) <= 0) f->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 459 j=pattern->col_couplePattern->index[iptr_ij]; 460 f_ij=f->col_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_ijmainPattern->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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 510 j=pattern->col_couplePattern->index[iptr_ij]; 511 f_ij=flux_matrix->col_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->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 526 j=pattern->col_couplePattern->index[iptr_ij]; 527 f_ij=flux_matrix->col_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 }