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

Revision 1408 - (show annotations)
Mon Feb 4 07:19:50 2008 UTC (13 years, 7 months ago) by gross
File MIME type: text/plain
File size: 14864 byte(s)
some problems with openmp fixed.

 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->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 131 j=pattern->couplePattern->index[iptr_ij]; 132 l_ij=L->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=MIN(u_max_i,u_j); 167 } 168 #pragma ivdep 169 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 170 j=pattern->couplePattern->index[iptr_ij]; 171 u_j=remote_u[j]; 172 u_min_i=MIN(u_min_i,u_j); 173 u_max_i=MIN(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} + = (a*m_{ij} + b* 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_updateAntiDiffusionFlux(const Paso_FCTransportProblem * fc, Paso_SystemMatrix *flux_matrix,const double a, const double b, const double* u) 188 { 189 dim_t n, j, i; 190 double *remote_u; 191 index_t iptr_ij; 192 const double b2=-b; 193 register double m_ij, ml_ij, k_ij, u_i; 194 Paso_SystemMatrixPattern *pattern; 195 Paso_SystemMatrix_startCollect(fc->iteration_matrix,u); 196 n=Paso_SystemMatrix_getTotalNumRows(fc->iteration_matrix); 197 pattern=fc->iteration_matrix->pattern; 198 remote_u=Paso_SystemMatrix_finishCollect(fc->iteration_matrix); 199 200 if ( (ABS(a) >0 ) ) { 201 if ( (ABS(b) >0 ) ) { 202 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij,k_ij,ml_ij) 203 for (i = 0; i < n; ++i) { 204 u_i=u[i]; 205 #pragma ivdep 206 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ijmainPattern->ptr[i+1]; ++iptr_ij) { 207 j=pattern->mainPattern->index[iptr_ij]; 208 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij]; 209 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij]; 210 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij]; 211 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij)); 212 } 213 #pragma ivdep 214 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 215 j=pattern->couplePattern->index[iptr_ij]; 216 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij]; 217 k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij]; 218 ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij]; 219 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*(a*m_ij+b2*(ml_ij+k_ij)); 220 } 221 } 222 } else { 223 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,m_ij) 224 for (i = 0; i < n; ++i) { 225 u_i=u[i]; 226 #pragma ivdep 227 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ijmainPattern->ptr[i+1]; ++iptr_ij) { 228 j=pattern->mainPattern->index[iptr_ij]; 229 m_ij=fc->mass_matrix->mainBlock->val[iptr_ij]; 230 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*a*m_ij; 231 } 232 #pragma ivdep 233 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 234 j=pattern->couplePattern->index[iptr_ij]; 235 m_ij=fc->mass_matrix->coupleBlock->val[iptr_ij]; 236 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*a*m_ij; 237 } 238 } 239 240 241 } 242 } else { 243 if ( (ABS(b) >0 ) ) { 244 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j,k_ij, ml_ij) 245 for (i = 0; i < n; ++i) { 246 u_i=u[i]; 247 #pragma ivdep 248 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ijmainPattern->ptr[i+1]; ++iptr_ij) { 249 j=pattern->mainPattern->index[iptr_ij]; 250 k_ij=fc->transport_matrix->mainBlock->val[iptr_ij]; 251 ml_ij=fc->iteration_matrix->mainBlock->val[iptr_ij]; 252 flux_matrix->mainBlock->val[iptr_ij]+=(u[j]-u_i)*b2*(ml_ij+k_ij); 253 } 254 #pragma ivdep 255 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 256 j=pattern->couplePattern->index[iptr_ij]; 257 k_ij=fc->transport_matrix->coupleBlock->val[iptr_ij]; 258 ml_ij=fc->iteration_matrix->coupleBlock->val[iptr_ij]; 259 flux_matrix->coupleBlock->val[iptr_ij]+=(remote_u[j]-u_i)*b2*(ml_ij+k_ij); 260 } 261 } 262 263 } 264 } 265 } 266 /* 267 * apply pre flux-correction: f_{ij}:=0 if f_{ij}*(u_[i]-u[j])<=0 268 * 269 */ 270 void Paso_FCTransportProblem_applyPreAntiDiffusionCorrection(Paso_SystemMatrix *f,const double* u) 271 { 272 dim_t n, i, j; 273 Paso_SystemMatrixPattern *pattern; 274 double *remote_u; 275 register double u_i, f_ij; 276 register index_t iptr_ij; 277 278 n=Paso_SystemMatrix_getTotalNumRows(f); 279 pattern=f->pattern; 280 remote_u=Paso_SystemMatrix_finishCollect(f); 281 #pragma omp parallel for schedule(static) private(i, u_i, iptr_ij, j, f_ij) 282 for (i = 0; i < n; ++i) { 283 u_i=u[i]; 284 #pragma ivdep 285 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ijmainPattern->ptr[i+1]; ++iptr_ij) { 286 j=pattern->mainPattern->index[iptr_ij]; 287 f_ij=f->mainBlock->val[iptr_ij]; 288 if (f_ij * (u_i-u[j]) <= 0) f->mainBlock->val[iptr_ij]=0; 289 } 290 #pragma ivdep 291 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 292 j=pattern->couplePattern->index[iptr_ij]; 293 f_ij=f->coupleBlock->val[iptr_ij]; 294 if (f_ij * (u_i-remote_u[j]) <= 0) f->coupleBlock->val[iptr_ij]=0; 295 } 296 } 297 } 298 299 300 void Paso_FCTransportProblem_setRs(const Paso_SystemMatrix *f,const double* lumped_mass_matrix, 301 const double* QN,const double* QP,double* RN,double* RP) 302 { 303 dim_t n, i, j; 304 Paso_SystemMatrixPattern *pattern; 305 register double f_ij, PP_i, PN_i; 306 register index_t iptr_ij; 307 308 n=Paso_SystemMatrix_getTotalNumRows(f); 309 pattern=f->pattern; 310 #pragma omp parallel for schedule(static) private(i, iptr_ij, j, f_ij, PP_i, PN_i) 311 for (i = 0; i < n; ++i) { 312 PP_i=0.; 313 PN_i=0.; 314 #pragma ivdep 315 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ijmainPattern->ptr[i+1]; ++iptr_ij) { 316 j=pattern->mainPattern->index[iptr_ij]; 317 if (i != j ) { 318 f_ij=f->mainBlock->val[iptr_ij]; 319 if (f_ij <=0 ) { 320 PN_i+=f_ij; 321 } else { 322 PP_i+=f_ij; 323 } 324 } 325 } 326 #pragma ivdep 327 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 328 j=pattern->couplePattern->index[iptr_ij]; 329 f_ij=f->coupleBlock->val[iptr_ij]; 330 if (f_ij <=0 ) { 331 PN_i+=f_ij; 332 } else { 333 PP_i+=f_ij; 334 } 335 } 336 if (PN_i<0) { 337 RN[i]=MIN(1,QN[i]*lumped_mass_matrix[i]/PN_i); 338 } else { 339 RN[i]=1.; 340 } 341 if (PP_i>0) { 342 RP[i]=MIN(1,QP[i]*lumped_mass_matrix[i]/PP_i); 343 } else { 344 RP[i]=1.; 345 } 346 } 347 348 } 349 350 void Paso_FCTransportProblem_addCorrectedFluxes(double* f,Paso_SystemMatrix *flux_matrix,const double* RN,const double* RP) 351 { 352 dim_t n, i, j; 353 Paso_SystemMatrixPattern *pattern; 354 double *remote_RN, *remote_RP; 355 register double RN_i, RP_i, f_i, f_ij; 356 register index_t iptr_ij; 357 n=Paso_SystemMatrix_getTotalNumRows(flux_matrix); 358 359 Paso_SystemMatrix_startCollect(flux_matrix,RN); 360 pattern=flux_matrix->pattern; 361 remote_RN=Paso_SystemMatrix_finishCollect(flux_matrix); 362 #pragma omp parallel for schedule(static) private(i, RN_i, RP_i, iptr_ij, j, f_ij, f_i) 363 for (i = 0; i < n; ++i) { 364 RN_i=RN[i]; 365 RP_i=RP[i]; 366 f_i=0; 367 #pragma ivdep 368 for (iptr_ij=(pattern->mainPattern->ptr[i]);iptr_ijmainPattern->ptr[i+1]; ++iptr_ij) { 369 j=pattern->mainPattern->index[iptr_ij]; 370 f_ij=flux_matrix->mainBlock->val[iptr_ij]; 371 if (f_ij >=0) { 372 f_i+=f_ij*MIN(RP_i,RN[j]); 373 } else { 374 f_i+=f_ij*MIN(RN_i,RP[j]); 375 } 376 } 377 #pragma ivdep 378 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 379 j=pattern->couplePattern->index[iptr_ij]; 380 f_ij=flux_matrix->coupleBlock->val[iptr_ij]; 381 if (f_ij >=0) { 382 f_i+=f_ij*MIN(RP_i,remote_RN[j]); 383 } 384 } 385 f[i]=f_i; 386 } 387 Paso_SystemMatrix_startCollect(flux_matrix,RP); 388 remote_RP=Paso_SystemMatrix_finishCollect(flux_matrix); 389 #pragma omp parallel for schedule(static) private(i, RN_i, iptr_ij, j, f_ij, f_i) 390 for (i = 0; i < n; ++i) { 391 RN_i=RN[i]; 392 f_i=0; 393 #pragma ivdep 394 for (iptr_ij=(pattern->couplePattern->ptr[i]);iptr_ijcouplePattern->ptr[i+1]; ++iptr_ij) { 395 j=pattern->couplePattern->index[iptr_ij]; 396 f_ij=flux_matrix->coupleBlock->val[iptr_ij]; 397 if (f_ij < 0) { 398 f_i+=f_ij*MIN(RN_i,remote_RP[j]); 399 } 400 } 401 f[i]+=f_i; 402 } 403 }