Changeset 6454
 Timestamp:
 Mar 4, 2009, 3:37:20 PM (14 years ago)
 Location:
 anuga_work/development/anuga_1d
 Files:

 2 added
 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/anuga_1d/channel_domain_ext.c
r6206 r6454 46 46 47 47 48 48 //WELL BALANCED VERSION 49 49 //Innermost flux function (using w=z+h) 50 int _flux_function_channel (double *q_leftm,double *q_leftp, double *q_rightm,50 int _flux_function_channel21(double *q_leftm,double *q_leftp, double *q_rightm, 51 51 double *q_rightp, double g, double epsilon, double h0, double *edgeflux, double *max_speed) { 52 52 … … 60 60 double zphalf,zmhalf,hleftstar,hrightstar; 61 61 double fluxtemp1,fluxtemp0,speedtemp; 62 double batemp ;62 double batemp,bphalf,bmhalf; 63 63 64 64 zmhalf = max(q_leftm[2],q_leftp[2]); … … 100 100 hrightstar = q_rightm[3]; 101 101 102 bphalf = 0.5*(b_rightm+b_rightp); 103 bmhalf = 0.5*(b_leftm+b_leftp); 104 //bphalf = min(b_rightm,b_rightp); 105 //bmhalf = min(b_leftm,b_leftp); 106 107 102 108 soundspeed_leftp = sqrt(g*h_leftp); 103 109 soundspeed_leftm = sqrt(g*h_leftm); … … 119 125 // Flux formulas for left hand side 120 126 flux_left[0] = d_leftm; 121 flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b _leftm;127 flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*bmhalf; 122 128 123 129 flux_right[0] = d_leftp; 124 flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b _leftp;130 flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*bmhalf; 125 131 126 132 … … 132 138 edgeflux[0] = s_maxl*flux_left[0]  s_minl*flux_right[0]; 133 139 134 batemp = ( h_leftp)*b_leftp(h_leftm)*b_leftm;135 140 batemp = (q_leftp[3]+q_leftp[2])*bmhalf(q_leftm[3]+q_leftm[2])*bmhalf; 141 edgeflux[0] += s_maxl*s_minl*batemp; 136 142 edgeflux[0] /= denom; 137 143 edgeflux[1] = s_maxl*flux_left[1]  s_minl*flux_right[1]; … … 147 153 // Flux formulas for right hand side 148 154 flux_left[0] = d_rightm; 149 flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b _rightm;155 flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*bphalf; 150 156 151 157 flux_right[0] = d_rightp; 152 flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b _rightp;158 flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*bphalf; 153 159 154 160 … … 161 167 edgeflux[0] = s_maxr*flux_left[0]  s_minr*flux_right[0]; 162 168 163 batemp = ( h_rightp)*b_rightp(h_rightm)*b_rightm;164 165 169 batemp = (q_rightp[3]+q_rightp[2])*bphalf(q_rightm[3]+q_rightm[2])*bphalf; 170 171 edgeflux[0] += s_maxr*s_minr*batemp; 166 172 edgeflux[0] /= denom; 167 173 edgeflux[1] = s_maxr*flux_left[1]  s_minr*flux_right[1]; … … 177 183 178 184 179 edgeflux[1]=0.5*g*h_rightm*h_rightm*b_rightm0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp0.5*g*h_leftp*h_leftp*b_leftp; 180 181 182 //edgeflux[1]=0.5*g*h_rightm*h_rightm0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar0.5*g*h_leftp*h_leftp; 185 edgeflux[1]=0.5*g*h_rightm*h_rightm*bphalf0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp0.5*g*h_leftp*h_leftp*bmhalf; 186 187 // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightmhleftstar*hleftstar*g*0.5*b_leftp); 188 189 edgeflux[1]=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp2*z_rightm)+hleftstar*(hrightstar+z_leftpz_rightm))b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstarz_leftp+z_rightm)+hleftstar*(hrightstar2*z_leftp+2*z_rightm))); 190 191 192 //edgeflux[1]=0.5*g*h_rightm*h_rightm0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar0.5*g*h_leftp*h_leftp; 183 193 184 194 //edgeflux[1]=0.5*g*b_rightm*h_rightm*h_rightm0.5*g*b_leftp*h_leftp*h_leftp; … … 196 206 return 0; 197 207 } 198 199 200 208 // GOOD BUT NOT WELL BALANCED VERSION 209 int _flux_function_channel(double *q_leftm,double *q_leftp, double *q_rightm, 210 double *q_rightp, double g, double epsilon, double h0, double *edgeflux, double *max_speed){ 211 212 int i; 213 double flux_left[2], flux_right[2]; 214 double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm; 215 double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp; 216 double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm; 217 double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp; 218 double s_maxl, s_minl,s_maxr,s_minr, denom; 219 double zphalf,zmhalf,hleftstar,hrightstar; 220 double fluxtemp1,fluxtemp0,speedtemp; 221 double batemp,bphalf,bmhalf; 222 223 zmhalf = max(q_leftm[2],q_leftp[2]); 224 zphalf = max(q_rightm[2],q_rightp[2]); 225 226 a_leftm = q_leftm[0]; 227 d_leftm = q_leftm[1]; 228 z_leftm = q_leftm[2]; 229 h_leftm = q_leftm[3]; 230 u_leftm = q_leftm[4]; 231 b_leftm = q_leftm[5]; 232 w_leftm = h_leftm+z_leftm; 233 234 a_leftp = q_leftp[0]; 235 d_leftp = q_leftp[1]; 236 z_leftp = q_leftp[2]; 237 h_leftp = q_leftp[3]; 238 u_leftp = q_leftp[4]; 239 b_leftp = q_leftp[5]; 240 w_leftp = h_leftp+z_leftp; 241 242 a_rightm = q_rightm[0]; 243 d_rightm = q_rightm[1]; 244 z_rightm = q_rightm[2]; 245 h_rightm = q_rightm[3]; 246 u_rightm = q_rightm[4]; 247 b_rightm = q_rightm[5]; 248 w_rightm = h_rightm+z_rightm; 249 250 a_rightp = q_rightp[0]; 251 d_rightp = q_rightp[1]; 252 z_rightp = q_rightp[2]; 253 h_rightp = q_rightp[3]; 254 u_rightp = q_rightp[4]; 255 b_rightp = q_rightp[5]; 256 w_rightp = h_rightp+z_rightp; 257 258 hleftstar = q_leftp[3]; 259 hrightstar = q_rightm[3]; 260 261 bphalf = min(b_rightm,b_rightp); 262 bmhalf = min(b_leftm,b_leftp); 263 264 soundspeed_leftp = sqrt(g*h_leftp); 265 soundspeed_leftm = sqrt(g*h_leftm); 266 soundspeed_rightp = sqrt(g*h_rightp); 267 soundspeed_rightm = sqrt(g*h_rightm); 268 269 s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp); 270 if (s_maxl < 0.0) s_maxl = 0.0; 271 272 s_minl = min(u_leftmsoundspeed_leftm, u_leftpsoundspeed_leftp); 273 if (s_minl > 0.0) s_minl = 0.0; 274 275 s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp); 276 if (s_maxr < 0.0) s_maxr = 0.0; 277 278 s_minr = min(u_rightmsoundspeed_rightm, u_rightpsoundspeed_rightp); 279 if (s_minr > 0.0) s_minr = 0.0; 280 281 // Flux formulas for left hand side 282 flux_left[0] = d_leftm; 283 flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm; 284 285 flux_right[0] = d_leftp; 286 flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp; 287 288 289 // Flux computation for left hand side 290 denom = s_maxls_minl; 291 if (denom < epsilon) { 292 for (i=0; i<2; i++) edgeflux[i] = 0.0; 293 } else { 294 edgeflux[0] = s_maxl*flux_left[0]  s_minl*flux_right[0]; 295 296 batemp = (q_leftp[3]+q_leftp[2])*b_leftp(q_leftm[3]+q_leftm[2])*b_leftm; 297 edgeflux[0] += s_maxl*s_minl*batemp; 298 edgeflux[0] /= denom; 299 edgeflux[1] = s_maxl*flux_left[1]  s_minl*flux_right[1]; 300 edgeflux[1] += s_maxl*s_minl*(d_leftpd_leftm); 301 edgeflux[1] /= denom; 302 303 304 } 305 fluxtemp0 = edgeflux[0]; 306 fluxtemp1 = edgeflux[1]; 307 308 309 // Flux formulas for right hand side 310 flux_left[0] = d_rightm; 311 flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b_rightm; 312 313 flux_right[0] = d_rightp; 314 flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b_rightp; 315 316 317 318 // Flux computation for right hand side 319 denom = s_maxrs_minr; 320 if (denom < epsilon) { 321 for (i=0; i<2; i++) edgeflux[i] = 0.0; 322 } else { 323 edgeflux[0] = s_maxr*flux_left[0]  s_minr*flux_right[0]; 324 325 batemp = (q_rightp[3]+q_rightp[2])*b_rightp(q_rightm[3]+q_rightm[2])*b_rightm; 326 327 edgeflux[0] += s_maxr*s_minr*batemp; 328 edgeflux[0] /= denom; 329 edgeflux[1] = s_maxr*flux_left[1]  s_minr*flux_right[1]; 330 edgeflux[1] += s_maxr*s_minr*(d_rightpd_rightm); 331 edgeflux[1] /= denom; 332 333 334 } 335 336 337 edgeflux[0]=edgeflux[0]fluxtemp0; 338 edgeflux[1]=edgeflux[1]fluxtemp1; 339 340 edgeflux[1]=0.5*0.5*g*(h_rightm+h_leftp)*(b_rightm+b_leftp)*(z_rightmz_leftp)+0.5*(h_rightm+h_leftp)*(h_rightm+h_leftp)*0.5*0.5*(b_rightmb_leftp)*g; 341 342 //edgeflux[1]=0.5*g*h_rightm*h_rightm*bphalf0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp0.5*g*h_leftp*h_leftp*bmhalf; 343 344 // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightmhleftstar*hleftstar*g*0.5*b_leftp); 345 346 //edgeflux[1]=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp2*z_rightm)+hleftstar*(hrightstar+z_leftpz_rightm))b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstarz_leftp+z_rightm)+hleftstar*(hrightstar2*z_leftp+2*z_rightm))); 347 348 349 //edgeflux[1]=0.5*g*h_rightm*h_rightm0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar0.5*g*h_leftp*h_leftp; 350 351 //edgeflux[1]=0.5*g*b_rightm*h_rightm*h_rightm0.5*g*b_leftp*h_leftp*h_leftp; 352 // Maximal wavespeed 353 if ( (s_maxls_minl)<epsilon && (s_maxrs_minr)<epsilon ){ 354 *max_speed = 0.0; 355 }else{ 356 speedtemp = max(fabs(s_maxl),fabs(s_minl)); 357 speedtemp = max(speedtemp,fabs(s_maxr)); 358 speedtemp = max(speedtemp,fabs(s_minr)); 359 *max_speed = speedtemp; 360 } 361 362 //printf("%f\n",h_right); 363 return 0; 364 } 365 // NAIEVE VERSION 366 int _flux_function_channel2(double *q_leftm,double *q_leftp, double *q_rightm, 367 double *q_rightp, double g, double epsilon, double h0, double *edgeflux, double *max_speed){ 368 369 int i; 370 double flux_left[2], flux_right[2]; 371 double a_leftm,w_leftm, h_leftm, d_leftm, z_leftm, u_leftm, b_leftm, soundspeed_leftm; 372 double a_leftp,w_leftp, h_leftp, d_leftp, z_leftp, u_leftp, b_leftp, soundspeed_leftp; 373 double a_rightm,w_rightm, h_rightm, d_rightm, z_rightm, u_rightm, b_rightm, soundspeed_rightm; 374 double a_rightp,w_rightp, h_rightp, d_rightp, z_rightp, u_rightp, b_rightp, soundspeed_rightp; 375 double s_maxl, s_minl,s_maxr,s_minr, denom; 376 double zphalf,zmhalf,hleftstar,hrightstar; 377 double fluxtemp1,fluxtemp0,speedtemp; 378 double batemp,bphalf,bmhalf; 379 380 zmhalf = max(q_leftm[2],q_leftp[2]); 381 zphalf = max(q_rightm[2],q_rightp[2]); 382 383 a_leftm = q_leftm[0]; 384 d_leftm = q_leftm[1]; 385 z_leftm = q_leftm[2]; 386 h_leftm = q_leftm[3]; 387 u_leftm = q_leftm[4]; 388 b_leftm = q_leftm[5]; 389 w_leftm = h_leftm+z_leftm; 390 391 a_leftp = q_leftp[0]; 392 d_leftp = q_leftp[1]; 393 z_leftp = q_leftp[2]; 394 h_leftp = q_leftp[3]; 395 u_leftp = q_leftp[4]; 396 b_leftp = q_leftp[5]; 397 w_leftp = h_leftp+z_leftp; 398 399 a_rightm = q_rightm[0]; 400 d_rightm = q_rightm[1]; 401 z_rightm = q_rightm[2]; 402 h_rightm = q_rightm[3]; 403 u_rightm = q_rightm[4]; 404 b_rightm = q_rightm[5]; 405 w_rightm = h_rightm+z_rightm; 406 407 a_rightp = q_rightp[0]; 408 d_rightp = q_rightp[1]; 409 z_rightp = q_rightp[2]; 410 h_rightp = q_rightp[3]; 411 u_rightp = q_rightp[4]; 412 b_rightp = q_rightp[5]; 413 w_rightp = h_rightp+z_rightp; 414 415 hleftstar = q_leftp[3]; 416 hrightstar = q_rightm[3]; 417 418 bphalf = min(b_rightm,b_rightp); 419 bmhalf = min(b_leftm,b_leftp); 420 421 soundspeed_leftp = sqrt(g*h_leftp); 422 soundspeed_leftm = sqrt(g*h_leftm); 423 soundspeed_rightp = sqrt(g*h_rightp); 424 soundspeed_rightm = sqrt(g*h_rightm); 425 426 s_maxl = max(u_leftm+soundspeed_leftm, u_leftp+soundspeed_leftp); 427 if (s_maxl < 0.0) s_maxl = 0.0; 428 429 s_minl = min(u_leftmsoundspeed_leftm, u_leftpsoundspeed_leftp); 430 if (s_minl > 0.0) s_minl = 0.0; 431 432 s_maxr = max(u_rightm+soundspeed_rightm, u_rightp+soundspeed_rightp); 433 if (s_maxr < 0.0) s_maxr = 0.0; 434 435 s_minr = min(u_rightmsoundspeed_rightm, u_rightpsoundspeed_rightp); 436 if (s_minr > 0.0) s_minr = 0.0; 437 438 // Flux formulas for left hand side 439 flux_left[0] = d_leftm; 440 flux_left[1] = u_leftm*d_leftm + 0.5*g*h_leftm*h_leftm*b_leftm; 441 442 flux_right[0] = d_leftp; 443 flux_right[1] = u_leftp*d_leftp + 0.5*g*h_leftp*h_leftp*b_leftp; 444 445 446 // Flux computation for left hand side 447 denom = s_maxls_minl; 448 if (denom < epsilon) { 449 for (i=0; i<2; i++) edgeflux[i] = 0.0; 450 } else { 451 edgeflux[0] = s_maxl*flux_left[0]  s_minl*flux_right[0]; 452 453 batemp = (q_leftp[3]+q_leftp[2])*b_leftp(q_leftm[3]+q_leftm[2])*b_leftm; 454 455 edgeflux[0] = 0.5*(flux_left[0]+flux_right[0]); 456 edgeflux[1] = 0.5*(flux_left[1]+flux_right[1]); 457 458 459 460 } 461 fluxtemp0 = edgeflux[0]; 462 fluxtemp1 = edgeflux[1]; 463 464 465 // Flux formulas for right hand side 466 flux_left[0] = d_rightm; 467 flux_left[1] = u_rightm*d_rightm + 0.5*g*h_rightm*h_rightm*b_rightm; 468 469 flux_right[0] = d_rightp; 470 flux_right[1] = u_rightp*d_rightp + 0.5*g*h_rightp*h_rightp*b_rightp; 471 472 473 474 // Flux computation for right hand side 475 denom = s_maxrs_minr; 476 if (denom < epsilon) { 477 for (i=0; i<2; i++) edgeflux[i] = 0.0; 478 } else { 479 edgeflux[0] = s_maxr*flux_left[0]  s_minr*flux_right[0]; 480 481 batemp = (q_rightp[3]+q_rightp[2])*b_rightp(q_rightm[3]+q_rightm[2])*b_rightm; 482 483 484 edgeflux[0] = 0.5*(flux_right[0]+flux_left[0]); 485 486 487 edgeflux[1] = 0.5*(flux_left[1]+flux_right[1]); 488 489 490 } 491 492 493 edgeflux[0]=edgeflux[0]fluxtemp0; 494 edgeflux[1]=edgeflux[1]fluxtemp1; 495 496 edgeflux[1]=0.5*0.5*g*(h_rightm+h_leftp)*(b_rightm+b_leftp)*(z_rightmz_leftp)+0.5*(h_rightm+h_leftp)*(h_rightm+h_leftp)*0.5*0.5*(b_rightmb_leftp)*g; 497 498 //edgeflux[1]=0.5*g*h_rightm*h_rightm*bphalf0.5*g*hrightstar*hrightstar*b_rightm+0.5*g*hleftstar*hleftstar*b_leftp0.5*g*h_leftp*h_leftp*bmhalf; 499 500 // printf("edgflux:%f expected:%f \n",edgeflux[1],hrightstar*hrightstar*g*0.5*b_rightmhleftstar*hleftstar*g*0.5*b_leftp); 501 502 //edgeflux[1]=g*(1.0/6.0)*(b_rightm*(hleftstar*hleftstar+hrightstar*(hrightstar+2*z_leftp2*z_rightm)+hleftstar*(hrightstar+z_leftpz_rightm))b_leftp*(hleftstar*hleftstar+hrightstar*(hrightstarz_leftp+z_rightm)+hleftstar*(hrightstar2*z_leftp+2*z_rightm))); 503 504 505 //edgeflux[1]=0.5*g*h_rightm*h_rightm0.5*g*hrightstar*hrightstar+0.5*g*hleftstar*hleftstar0.5*g*h_leftp*h_leftp; 506 507 //edgeflux[1]=0.5*g*b_rightm*h_rightm*h_rightm0.5*g*b_leftp*h_leftp*h_leftp; 508 // Maximal wavespeed 509 if ( (s_maxls_minl)<epsilon && (s_maxrs_minr)<epsilon ){ 510 *max_speed = 0.0; 511 }else{ 512 speedtemp = max(fabs(s_maxl),fabs(s_minl)); 513 speedtemp = max(speedtemp,fabs(s_maxr)); 514 speedtemp = max(speedtemp,fabs(s_minr)); 515 *max_speed = speedtemp; 516 } 517 518 //printf("%f\n",h_right); 519 return 0; 520 } 201 521 202 522
Note: See TracChangeset
for help on using the changeset viewer.