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

Revision 3005 - (show annotations)
Thu Apr 22 05:59:31 2010 UTC (9 years, 5 months ago) by gross
File MIME type: text/plain
File size: 18994 byte(s)
early call of setPreconditioner in FCT solver removed.

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