Contents of /branches/stage3.1/paso/src/SolverFCT_FluxControl.c

Revision 2945 - (show annotations)
Wed Feb 24 00:17:46 2010 UTC (8 years, 11 months ago) by jfenwick
File MIME type: text/plain
File size: 19021 byte(s)
Bringing release stage up to trunk version 2944


 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 "SolverFCT.h" 28 #include "PasoUtil.h" 29 30 31 /**************************************************************/ 32 33 /* 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 46 void Paso_FCTransportProblem_setLowOrderOperator(Paso_FCTransportProblem * fc) { 47 dim_t n=Paso_SystemMatrix_getTotalNumRows(fc->transport_matrix),i; 48 index_t iptr_ij,j,iptr_ji; 49 Paso_SystemMatrixPattern *pattern; 50 register double d_ij, sum, rtmp1, rtmp2; 51 52 if (fc==NULL) return; 53 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, TRUE); 58 } 59 60 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 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 #pragma ivdep 72 for (iptr_ji=pattern->mainPattern->ptr[j]; iptr_jimainPattern->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 } 83 } 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 #pragma ivdep 88 for (iptr_ji=pattern->row_couplePattern->ptr[j]; iptr_jirow_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 /* set main diagonal entry */ 101 fc->main_diagonal_low_order_transport_matrix[i]=sum; 102 } 103 } 104 } 105 /* 106 * out_i=m_i u_i + a * \sum_{j <> i} l_{ij} (u_j-u_i) + b * q_i 107 * 108 */ 109 void Paso_SolverFCT_setMuPaLuPbQ(double* out, 110 const double* M, 111 const Paso_Coupler* u_coupler, 112 const double a, 113 const Paso_SystemMatrix *L, 114 const double b, 115 const double* Q) 116 { 117 dim_t n, i, j; 118 Paso_SystemMatrixPattern *pattern; 119 const double *u=Paso_Coupler_borrowLocalData(u_coupler); 120 const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler); 121 register double sum, u_i, l_ij; 122 register index_t iptr_ij; 123 n=Paso_SystemMatrix_getTotalNumRows(L); 124 125 #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 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 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 void Paso_SolverFCT_setQs(const Paso_Coupler* u_coupler,double* QN, double* QP, const Paso_SystemMatrix *L) 157 { 158 dim_t n, i, j; 159 Paso_SystemMatrixPattern *pattern; 160 const double *u=Paso_Coupler_borrowLocalData(u_coupler); 161 const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler); 162 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 u_max_i=MAX(u_max_i,u_j); 177 } 178 #pragma ivdep 179 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 u_j=remote_u[j]; 182 u_min_i=MIN(u_min_i,u_j); 183 u_max_i=MAX(u_max_i,u_j); 184 } 185 QN[i]=u_min_i-u_i; 186 QP[i]=u_max_i-u_i; 187 } 188 } 189 190 /* 191 * 192 * 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 void Paso_FCTransportProblem_setAntiDiffusionFlux(const double dt, const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix, const Paso_Coupler* u_coupler) 198 { 199 dim_t n, j, i; 200 index_t iptr_ij; 201 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 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 const Paso_SystemMatrixPattern *pattern=fc->iteration_matrix->pattern; 209 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 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 d_u=remote_u[j]-u_i; 233 d_u_last=remote_u_last[j]-u_last_i; 234 flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last- (m_ij+f2*d_ij)*d_u; 235 } 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 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 d_u=remote_u[j]-u_i; 259 d_u_last=remote_u_last[j]-u_last_i; 260 flux_matrix->col_coupleBlock->val[iptr_ij]=(m_ij+f1*d_ij)*d_u_last-m_ij*d_u; 261 } 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 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 d_u=remote_u[j]-u_i; 287 d_u_last=remote_u_last[j]-u_last_i; 288 flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*d_u_last- (m_ij+f2*d_ij)*d_u; 289 } 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 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 d_u=remote_u[j]-u_i; 313 d_u_last=remote_u_last[j]-u_last_i; 314 flux_matrix->col_coupleBlock->val[iptr_ij]=m_ij*(d_u_last-d_u); 315 } 316 } 317 } 318 } 319 } 320 /* 321 * apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0 322 * 323 */ 324 void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const Paso_Coupler* u_coupler) 325 { 326 dim_t n, i, j; 327 Paso_SystemMatrixPattern *pattern; 328 const double *u=Paso_Coupler_borrowLocalData(u_coupler); 329 const double *remote_u=Paso_Coupler_borrowRemoteData(u_coupler); 330 register double u_i, f_ij; 331 register index_t iptr_ij; 332 333 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 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 } 350 } 351 } 352 353 354 void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix, 355 const Paso_Coupler* QN_coupler,const Paso_Coupler* QP_coupler,double* RN,double* RP) 356 { 357 dim_t n, i, j; 358 Paso_SystemMatrixPattern *pattern; 359 const double *QN=Paso_Coupler_borrowLocalData(QN_coupler); 360 const double *QP=Paso_Coupler_borrowLocalData(QP_coupler); 361 register double f_ij, PP_i, PN_i; 362 register index_t iptr_ij; 363 364 n=Paso_SystemMatrix_getTotalNumRows(f); 365 pattern=f->pattern; 366 #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i) 367 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 if (f_ij <=0) { 376 PN_i+=f_ij; 377 } else { 378 PP_i+=f_ij; 379 } 380 } 381 } 382 #pragma ivdep 383 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 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 404 } 405 406 void Paso_FCTransportProblem_addCorrectedFluxes(double* f,const Paso_SystemMatrix *flux_matrix, const Paso_Coupler* RN_coupler, const Paso_Coupler* RP_coupler) 407 { 408 dim_t n, i, j; 409 Paso_SystemMatrixPattern *pattern; 410 register double RN_i, RP_i, f_i, f_ij; 411 double alpha_i; 412 register index_t iptr_ij; 413 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 n=Paso_SystemMatrix_getTotalNumRows(flux_matrix); 418 419 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 423 alpha_i=0; 424 RN_i=RN[i]; 425 RP_i=RP[i]; 426 f_i=0; 427 #pragma ivdep 428 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ijmainPattern->ptr[i+1]; ++iptr_ij) { 429 j=pattern->mainPattern->index[iptr_ij]; 430 f_ij=flux_matrix->mainBlock->val[iptr_ij]; 431 if (f_ij >=0) { 432 f_i+=f_ij*MIN(RP_i,RN[j]); 433 /* printf("alpha %d %d = %e %e\n",i,j,MIN(RP_i,RN[j]),f_ij); 434 alpha_i=MAX(alpha_i,MIN(RP_i,RN[j]) ); */ 435 } else { 436 f_i+=f_ij*MIN(RN_i,RP[j]); 437 /* 438 printf("alpha %d %d = %e %e\n",i,j,MIN(RN_i,RP[j]),f_ij); 439 alpha_i=MAX(alpha_i,MIN(RN_i,RP[j]) ); 440 */ 441 } 442 } 443 #pragma ivdep 444 for (iptr_ij=(pattern->col_couplePattern->ptr[i]);iptr_ijcol_couplePattern->ptr[i+1]; ++iptr_ij) { 445 j=pattern->col_couplePattern->index[iptr_ij]; 446 f_ij=flux_matrix->col_coupleBlock->val[iptr_ij]; 447 if (f_ij >=0) { 448 f_i+=f_ij*MIN(RP_i,remote_RN[j]); 449 }else { 450 f_i+=f_ij*MIN(RN_i,remote_RP[j]); 451 } 452 } 453 /* printf("alpha %d = %e -> %e (%e, %e)\n",i,alpha_i, f_i, f[i],f[i]+f_i); */ 454 f[i]+=f_i; 455 } 456 }