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

Revision 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 9 months ago) by caltinay
File MIME type: text/plain
File size: 18923 byte(s)
Assorted spelling/comment fixes in paso.


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