- Timestamp:
- May 12, 2013, 8:11:11 PM (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c
r8775 r8865 94 94 double g, 95 95 double *edgeflux, double *max_speed, 96 double *pressure_flux) 96 double *pressure_flux, double hc, 97 double hc_n, 98 double speed_max_last) 97 99 { 98 100 … … 141 143 h_left = max(w_left - z,0.); 142 144 uh_left = q_left_rotated[1]; 143 //if(h_left>1.e-06){144 //u_left = uh_left/max(h_left, 1.0e-06);145 //}else{146 //u_left=0.;147 //}148 u_left = _compute_speed(&uh_left, &h_left,149 epsilon, h0, limiting_threshold);145 if(h_left>h0){ 146 u_left = uh_left/max(h_left, 1.0e-06); 147 }else{ 148 u_left=0.; 149 } 150 //u_left = _compute_speed(&uh_left, &h_left, 151 // epsilon, h0, limiting_threshold); 150 152 151 153 w_right = q_right_rotated[0]; 152 154 h_right = max(w_right - z, 0.); 153 155 uh_right = q_right_rotated[1]; 154 u_right = _compute_speed(&uh_right, &h_right,155 epsilon, h0, limiting_threshold);156 //if(h_right>1.0e-06){157 //u_right = uh_right/max(h_right, 1.0e-06);158 //}else{159 //u_right=0.;160 //}156 //u_right = _compute_speed(&uh_right, &h_right, 157 // epsilon, h0, limiting_threshold); 158 if(h_right>h0){ 159 u_right = uh_right/max(h_right, 1.0e-06); 160 }else{ 161 u_right=0.; 162 } 161 163 162 164 // Momentum in y-direction … … 198 200 } 199 201 202 if( hc < h0){ 203 s_max = 0.0; 204 } 205 206 200 207 s_min = min(u_left - soundspeed_left, u_right - soundspeed_right); 201 208 if (s_min > 0.0) … … 204 211 } 205 212 213 if( hc_n < h0){ 214 s_min = 0.0; 215 } 216 206 217 // Flux formulas 207 218 flux_left[0] = u_left*h_left; 208 flux_left[1] = u_left*uh_left ; //+ 0.5*g*h_left*h_left; 209 flux_left[2] = u_left*vh_left; 219 //if(hc>h0){ 220 flux_left[1] = u_left*uh_left; //+ 0.5*g*h_left*h_left; 221 flux_left[2] = u_left*vh_left; 222 //}else{ 223 // flux_left[1] = 0.;//u_left*uh_left; //+ 0.5*g*h_left*h_left; 224 // flux_left[2] = 0.;//u_left*vh_left; 225 //} 210 226 211 227 flux_right[0] = u_right*h_right; 212 flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right; 213 flux_right[2] = u_right*vh_right; 228 //if(hc_n>h0){ 229 flux_right[1] = u_right*uh_right ; //+ 0.5*g*h_right*h_right; 230 flux_right[2] = u_right*vh_right; 231 //}else{ 232 // flux_right[1] = 0.; //u_right*uh_right ; //+ 0.5*g*h_right*h_right; 233 // flux_right[2] = 0.; //u_right*vh_right; 234 //} 214 235 215 236 // Flux computation … … 222 243 *max_speed = 0.0; 223 244 *pressure_flux = 0.0; 245 //*pressure_flux = 0.5*g*0.5*(h_left*h_left+h_right*h_right); 224 246 } 225 247 else 226 248 { 249 // Maximal wavespeed 250 *max_speed = max(s_max, -s_min); 251 227 252 inverse_denominator = 1.0/denom; 228 253 for (i = 0; i < 3; i++) 229 254 { 230 255 edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i]; 231 edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]); 256 // Standard smoothing term 257 //edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i]); 258 //GD HACK: Here, we supress the 'smoothing' part of the flux -- like scaling diffusion by local wave speed 259 edgeflux[i] += (s_max*s_min)*(q_right_rotated[i] - q_left_rotated[i])* 260 (min(*max_speed/(max(speed_max_last,1e-06)),1.0));//(min(hc/h_left,min(hc_n/h_right,1.0))); 232 261 edgeflux[i] *= inverse_denominator; 233 262 } … … 235 264 *pressure_flux = 0.5*g*( s_max*h_left*h_left -s_min*h_right*h_right)*inverse_denominator; 236 265 237 // Maximal wavespeed238 *max_speed = max(fabs(s_max), fabs(s_min));239 266 240 267 // Rotate back … … 274 301 double* max_speed_array, 275 302 int optimise_dry_cells, 303 int timestep_order, 276 304 double* stage_centroid_values, 277 305 double* xmom_centroid_values, … … 281 309 // Local variables 282 310 double max_speed, length, inv_area, zl, zr; 283 double h0 = H0* H0; // This ensures a good balance when h approaches H0.284 285 double limiting_threshold = H0; //10 * H0; // Avoid applying limiter below this311 double h0 = H0*2.0;//H0*H0;//H0*H0; // This ensures a good balance when h approaches H0. 312 313 double limiting_threshold = 10 * H0 ;//H0; //10 * H0; // Avoid applying limiter below this 286 314 // threshold for performance reasons. 287 315 // See ANUGA manual under flux limiting … … 297 325 double stage_edge_lim, outgoing_mass_edges, bedtop, bedbot, pressure_flux, hc, hc_n, tmp, tmp2; 298 326 static long call = 1; // Static local variable flagging already computed flux 327 double speed_max_last, vol; 299 328 300 329 double *max_bed_edgevalue, *min_bed_edgevalue, *edgeflux_store, *pressuregrad_store; … … 314 343 local_timestep=timestep; 315 344 345 //printf("timestep_order %lf \n", timestep_order*1.0); 316 346 // Compute minimum bed edge value on each triangle 347 speed_max_last=0.0; 317 348 for (k = 0; k < number_of_elements; k++){ 318 349 max_bed_edgevalue[k] = max(bed_edge_values[3*k], 319 350 max(bed_edge_values[3*k+1], bed_edge_values[3*k+2])); 351 speed_max_last=max(speed_max_last, max_speed_array[k]); 352 353 //if(stage_centroid_values[k]<bed_centroid_values[k]){ 354 // printf("Stage < bed"); 355 //} 320 356 } 321 357 358 //printf("SML: %.12lf, %.12lf, %.12lf \n", speed_max_last, fabs(-speed_max_last), max(speed_max_last*2.0, 0.5*speed_max_last)/2.0); 359 360 //printf("SPM: %lf \n", speed_max_last); 322 361 323 362 // For all triangles … … 375 414 //If one centroid is dry, then extrapolate its edge values from the neighbouring centroid, 376 415 // unless the local centroid value is smaller 377 if(n>=0){378 if((hc<=H0)&&(hc_n>H0)){379 ql[0]=max(min(qr[0],stage_centroid_values[k]),zl);380 }381 if((hc_n<=H0)&&(hc>H0)){382 qr[0]=max(min(ql[0],stage_centroid_values[n]),zr);383 }384 385 }else{386 // Treat the boundary case387 // ??388 }416 //if(n>=0){ 417 // if((hc<=H0)&&(hc_n>H0)){ 418 // ql[0]=max(min(qr[0],stage_centroid_values[k]),zl); 419 // } 420 // if((hc_n<=H0)&&(hc>H0)){ 421 // qr[0]=max(min(ql[0],stage_centroid_values[n]),zr); 422 // } 423 424 //}else{ 425 // // Treat the boundary case 426 // // ?? 427 //} 389 428 390 429 // Edge flux computation (triangle k, edge i) … … 392 431 normals[ki2], normals[ki2 + 1], 393 432 epsilon, h0, limiting_threshold, g, 394 edgeflux, &max_speed, &pressure_flux); 433 edgeflux, &max_speed, &pressure_flux, hc, hc_n, 434 speed_max_last); 395 435 396 436 … … 409 449 // bedslope_work contains all gravity related terms -- weighting of 410 450 // conservative and non-conservative versions. 411 bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux); 412 //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 413 pressuregrad_store[ki] = bedslope_work; 414 451 if(hc> -h0){ 452 bedslope_work = length*(g*(hc*ql[0]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux); 453 //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.5*max(ql[0]-zl,0.)*(ql[0]-zl)) + pressure_flux); 454 //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 455 pressuregrad_store[ki] = bedslope_work; 456 }else{ 457 bedslope_work = length*(g*(hc*ql[0]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux); 458 //bedslope_work = length*(g*(hc*stage_centroid_values[k]-0.0*max(ql[0]-zl,0.)*(ql[0]-zl)) + 0.0*pressure_flux); 459 //bedslope_work+= (1.0-tmp)*length*g*hc*ql[0]; // Non-conservative water surface slope 460 pressuregrad_store[ki] = bedslope_work; 461 } 415 462 416 463 already_computed_flux[ki] = call; // #k Done … … 424 471 edgeflux_store[nm3 + 2 ] = edgeflux[2]; 425 472 426 bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux); 427 //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0]; 428 pressuregrad_store[nm] = bedslope_work; 429 473 if(hc_n> -h0){ 474 bedslope_work = length*(g*(hc_n*qr[0]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux); 475 //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.5*max(qr[0]-zr,0.)*(qr[0]-zr)) + pressure_flux); 476 //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0]; 477 pressuregrad_store[nm] = bedslope_work; 478 }else{ 479 bedslope_work = length*(g*(hc_n*qr[0]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux); 480 //bedslope_work = length*(g*(hc_n*stage_centroid_values[n]-0.0*max(qr[0]-zr,0.)*(qr[0]-zr)) + 0.0*pressure_flux); 481 //bedslope_work+= (1.-tmp)*length*g*hc_n*qr[0]; 482 pressuregrad_store[nm] = bedslope_work; 483 } 430 484 431 485 already_computed_flux[nm] = call; // #n Done … … 439 493 if ((tri_full_flag[k] == 1) ) { 440 494 441 if(call% 2==1) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep495 if(call%timestep_order!=0) max_speed = max_speed_array[k]; // HACK to Ensure that local timestep is the same as the last timestep 442 496 443 497 if (max_speed > epsilon) { … … 465 519 // Limit edgefluxes, for mass conservation near wet/dry cells 466 520 for(k=0; k< number_of_elements; k++){ 521 //continue; 467 522 hc = max(stage_centroid_values[k] - bed_centroid_values[k],0.); 468 523 // Loop over every edge … … 474 529 if(edgeflux_store[3*(3*k+useint)]< 0.){ 475 530 //outgoing_mass_edges+=1.0; 476 outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)] *local_timestep);531 outgoing_mass_edges+=(edgeflux_store[3*(3*k+useint)]); 477 532 } 478 533 } 534 outgoing_mass_edges*=local_timestep; 479 535 } 480 536 … … 486 542 // Idea: The cell will not go dry if: 487 543 // total_outgoing_flux <= cell volume = Area_triangle*hc 488 if((edgeflux_store[ki3]< 0.0) && (outgoing_mass_edges<0.)){ 544 vol=areas[k]*hc; 545 if((edgeflux_store[ki3]< 0.0) && (-outgoing_mass_edges> vol)){ 489 546 490 547 // This bound could be improved (e.g. we could actually sum … … 493 550 // about subsequent changes to the + edgeflux caused by 494 551 // constraints associated with neighbouring triangles. 495 tmp = (areas[k]*hc)/(fabs(outgoing_mass_edges)+1.0e-10) ;552 tmp = vol/(-(outgoing_mass_edges)) ; 496 553 if(tmp< 1.0){ 497 554 edgeflux_store[ki3+0]*=tmp; … … 528 585 // GD HACK 529 586 // Option to limit advective fluxes 530 if(hc > H0){587 if(hc > h0){ 531 588 stage_explicit_update[k] += edgeflux_store[ki3+0]; 532 589 // Cut these terms for shallow flows, and protect stationary states from round-off error … … 546 603 // GD HACK 547 604 // Compute bed slope term 548 if(hc >H0){605 if(hc > h0){ 549 606 xmom_explicit_update[k] -= normals[ki2]*pressuregrad_store[ki]; 550 607 ymom_explicit_update[k] -= normals[ki2+1]*pressuregrad_store[ki]; … … 564 621 } // end cell k 565 622 566 if(call% 2==0) timestep=local_timestep;623 if(call%timestep_order==0) timestep=local_timestep; 567 624 568 625 // Look for 'dry' edges, and set momentum (& component of momentum_update) … … 668 725 669 726 //if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) { 670 if (hc < minimum_allowed_height ) {727 if (hc < minimum_allowed_height*2.0) { 671 728 // Set momentum to zero and ensure h is non negative 672 xmomc[k] = 0.0; 673 ymomc[k] = 0.0; 729 //xmomc[k] = 0.;//xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0; 730 //ymomc[k] = 0.;//ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0; 731 xmomc[k] = xmomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0; 732 ymomc[k] = ymomc[k]*hc*max(hc-minimum_allowed_height,0.)/(2.0*minimum_allowed_height*minimum_allowed_height);// 0.0; 674 733 675 734 if (hc <= 0.0){ … … 678 737 //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2])); 679 738 //bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])) -minimum_allowed_height; 680 bmin=zc[k] -minimum_allowed_height;739 bmin=zc[k];//-1.0e-03; 681 740 // Minimum allowed stage = bmin 682 741 … … 687 746 688 747 wc[k] = max(wc[k], bmin); 689 printf(" mass_add, %f, %f, %f,%f,%f,%d\n", xc[k], yc[k], wc[k],zc[k],wc[k]-zc[k], k) ;748 //printf(" mass_add, %f, %f, %f,%f,%f,%d\n", xc[k], yc[k], wc[k],zc[k],wc[k]-zc[k], k) ; 690 749 691 750 … … 801 860 double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle 802 861 double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2; 803 double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin;862 double dqv[3], qmin, qmax, hmin, hmax, bedmax,bedmin, stagemin; 804 863 double hc, h0, h1, h2, beta_tmp, hfactor, xtmp, ytmp, weight, tmp; 805 864 double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale, vect_norm, l1, l2, a_bs, b_bs, area_e, inv_area_e; … … 848 907 cell_wetness_scale[k] = 0.; 849 908 // Check if cell k is wet 850 if(stage_centroid_values[k] > elevation_centroid_values[k]){851 //if(stage_centroid_values[k] > elevation_centroid_values[k] +minimum_allowed_height+epsilon){909 //if(stage_centroid_values[k] > elevation_centroid_values[k]){ 910 if(stage_centroid_values[k] > elevation_centroid_values[k] + 2*minimum_allowed_height+epsilon){ 852 911 //if(stage_centroid_values[k] > max_elevation_edgevalue[k] + minimum_allowed_height+epsilon){ 853 912 cell_wetness_scale[k] = 1.; … … 935 994 936 995 // Compute bed slope as a reference 937 area_e = (yv2-yv0)*(xv1-xv0) - (yv1-yv0)*(xv2-xv0);938 inv_area_e=1.0/area_e;939 a_bs = (yv2 - yv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) -940 (yv1 - yv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);941 a_bs *= inv_area_e;942 943 b_bs = -(xv2 - xv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) +944 (xv1 - xv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]);945 b_bs *= inv_area_e;996 //area_e = (yv2-yv0)*(xv1-xv0) - (yv1-yv0)*(xv2-xv0); 997 //inv_area_e=1.0/area_e; 998 //a_bs = (yv2 - yv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) - 999 // (yv1 - yv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]); 1000 //a_bs *= inv_area_e; 1001 1002 //b_bs = -(xv2 - xv0)*(elevation_edge_values[k3+1]-elevation_edge_values[k3+0]) + 1003 // (xv1 - xv0)*(elevation_edge_values[k3+2]-elevation_edge_values[k3+0]); 1004 //b_bs *= inv_area_e; 946 1005 947 1006 //printf("slopes: %f, %f \n", a_bs, b_bs); … … 989 1048 // Take note if the max neighbour bed elevation is greater than the min 990 1049 // neighbour stage -- suggests a 'steep' bed relative to the flow 991 //bedmax = max(elevation_centroid_values[k], 992 // max(elevation_centroid_values[k0], 993 // max(elevation_centroid_values[k1], 994 // elevation_centroid_values[k2]))); 1050 bedmax = max(elevation_centroid_values[k], 1051 max(elevation_centroid_values[k0], 1052 max(elevation_centroid_values[k1], 1053 elevation_centroid_values[k2]))); 1054 bedmin = min(elevation_centroid_values[k], 1055 min(elevation_centroid_values[k0], 1056 min(elevation_centroid_values[k1], 1057 elevation_centroid_values[k2]))); 995 1058 //stagemin = min(max(stage_centroid_values[k], elevation_centroid_values[k]), 996 1059 // min(max(stage_centroid_values[k0], elevation_centroid_values[k0]), … … 1028 1091 1029 1092 if(0==1){ 1030 // Friction slope type extrapolation -- emulating steady flow1031 // d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth**(4/3)1032 hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0);1033 if(hc>1.0e-06){1034 tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] +1035 ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03;1036 tmp /=pow(hc,10.0/3.0);1037 }else{1038 tmp=0.0;1039 }1040 a = -xmom_centroid_values[k]*tmp;1041 b = -ymom_centroid_values[k]*tmp;1093 //// Friction slope type extrapolation -- emulating steady flow 1094 //// d(stage)/dx = - u*sqrt(u**2 + v**2)*n**2/depth**(4/3) 1095 //hc=max(stage_centroid_values[k] - elevation_centroid_values[k],0.0); 1096 //if(hc>1.0e-06){ 1097 // tmp = sqrt(xmom_centroid_values[k]*xmom_centroid_values[k] + 1098 // ymom_centroid_values[k]*ymom_centroid_values[k])*0.03*0.03; 1099 // tmp /=pow(hc,10.0/3.0); 1100 //}else{ 1101 // tmp=0.0; 1102 //} 1103 //a = -xmom_centroid_values[k]*tmp; 1104 //b = -ymom_centroid_values[k]*tmp; 1042 1105 1043 // Limit gradient to be between the bed slope, and zero1044 if(a*a_bs<0.) a=0.;1045 if(fabs(a)>fabs(a_bs)) a=a_bs;1046 if(b*b_bs<0.) b=0.;1047 if(fabs(b)>fabs(b_bs)) b=b_bs;1048 1049 1050 dqv[0] = a*dxv0 + b*dyv0;1051 dqv[1] = a*dxv1 + b*dyv1;1052 dqv[2] = a*dxv2 + b*dyv2;1053 1054 stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0];1055 stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1];1056 stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2];1106 //// Limit gradient to be between the bed slope, and zero 1107 //if(a*a_bs<0.) a=0.; 1108 //if(fabs(a)>fabs(a_bs)) a=a_bs; 1109 //if(b*b_bs<0.) b=0.; 1110 //if(fabs(b)>fabs(b_bs)) b=b_bs; 1111 1112 1113 //dqv[0] = a*dxv0 + b*dyv0; 1114 //dqv[1] = a*dxv1 + b*dyv1; 1115 //dqv[2] = a*dxv2 + b*dyv2; 1116 1117 //stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0]; 1118 //stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1]; 1119 //stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2]; 1057 1120 }else{ 1058 1121 // Constant stage extrapolation … … 1185 1248 ////if (hc>0.0) 1186 1249 //{ 1187 hfactor = 1.0 ;//hmin/(hmin + 0.004);1250 // hfactor = 1.0 ;//hmin/(hmin + 0.004); 1188 1251 //hfactor=hmin/(hmin + 0.004); 1189 1252 //} 1190 hfactor= min(2.0*max(hmin,0)/max(hc,1.0e-06), 1.0); 1253 hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height*10.)), 1.0); 1254 //hfactor= min(2.0*max(hmin,0.0)/max(hc,max(0.5*(bedmax-bedmin),minimum_allowed_height)), 1.0); 1191 1255 1192 1256 //----------------------------------- … … 1252 1316 //if( (area2>0) ){ 1253 1317 //if(count_wet_neighbours[k]>0){ 1254 1318 limit_gradient(dqv, qmin, qmax, beta_tmp); 1255 1319 //} 1256 1320 //} 1257 1258 1259 1321 stage_edge_values[k3+0] = stage_centroid_values[k] + dqv[0]; 1322 stage_edge_values[k3+1] = stage_centroid_values[k] + dqv[1]; 1323 stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2]; 1260 1324 // IDEA: extrapolate assuming slope is that of steady flow?? 1261 1325 // GD HACK - FIXME … … 1746 1810 1747 1811 double timestep, epsilon, H0, g; 1748 int optimise_dry_cells ;1812 int optimise_dry_cells, timestep_order; 1749 1813 1750 1814 // Convert Python arguments to C 1751 if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOi OOOOO",1815 if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOOOiiOOOOO", 1752 1816 ×tep, 1753 1817 &epsilon, … … 1774 1838 &max_speed_array, 1775 1839 &optimise_dry_cells, 1840 ×tep_order, 1776 1841 &stage_centroid_values, 1777 1842 &xmom_centroid_values, … … 1841 1906 (double*) max_speed_array -> data, 1842 1907 optimise_dry_cells, 1908 timestep_order, 1843 1909 (double*) stage_centroid_values -> data, 1844 1910 (double*) xmom_centroid_values -> data, … … 1871 1937 1872 1938 1873 h0 = H0*H0; // This ensures a good balance when h approaches H0. 1874 // But evidence suggests that h0 can be as little as 1875 // epsilon! 1876 1939 h0 = H0; 1877 1940 limiting_threshold = 10*H0; // Avoid applying limiter below this 1878 1941 // threshold for performance reasons. … … 1890 1953 (double*) edgeflux -> data, 1891 1954 &max_speed, 1892 &pressure_flux); 1955 &pressure_flux, 1956 ((double*) normal -> data)[0], 1957 ((double*) normal -> data)[1], 1958 ((double*) normal -> data)[1] 1959 ); 1893 1960 1894 1961 return Py_BuildValue("d", max_speed);
Note: See TracChangeset
for help on using the changeset viewer.
