Changeset 6902 for branches/numpy/anuga/shallow_water/shallow_water_ext.c
 Timestamp:
 Apr 24, 2009, 5:22:14 PM (14 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/numpy/anuga/shallow_water/shallow_water_ext.c
r6780 r6902 28 28 29 29 // Computational function for rotation 30 // FIXME: Perhaps inline this and profile 30 31 int _rotate(double *q, double n1, double n2) { 31 32 /*Rotate the momentum component q (q[1], q[2]) … … 52 53 53 54 int find_qmin_and_qmax(double dq0, double dq1, double dq2, 54 55 double *qmin, double *qmax){ 55 56 // Considering the centroid of an FV triangle and the vertices of its 56 57 // auxiliary triangle, find … … 67 68 if (dq1>=dq2){ 68 69 if (dq1>=0.0) 69 70 *qmax=dq0+dq1; 70 71 else 71 72 72 *qmax=dq0; 73 73 74 *qmin=dq0+dq2; 74 75 if (*qmin>=0.0) *qmin = 0.0; … … 76 77 else{// dq1<dq2 77 78 if (dq2>0) 78 79 *qmax=dq0+dq2; 79 80 else 80 81 82 *qmin=dq0+dq1; 81 *qmax=dq0; 82 83 *qmin=dq0+dq1; 83 84 if (*qmin>=0.0) *qmin=0.0; 84 85 } … … 87 88 if (dq1<=dq2){ 88 89 if (dq1<0.0) 89 90 *qmin=dq0+dq1; 90 91 else 91 92 93 *qmax=dq0+dq2; 92 *qmin=dq0; 93 94 *qmax=dq0+dq2; 94 95 if (*qmax<=0.0) *qmax=0.0; 95 96 } 96 97 else{// dq1>dq2 97 98 if (dq2<0.0) 98 99 *qmin=dq0+dq2; 99 100 else 100 101 101 *qmin=dq0; 102 102 103 *qmax=dq0+dq1; 103 104 if (*qmax<=0.0) *qmax=0.0; … … 133 134 134 135 phi=min(r*beta_w,1.0); 135 for (i=0;i<3;i++) 136 dqv[i]=dqv[i]*phi; 136 //for (i=0;i<3;i++) 137 dqv[0]=dqv[0]*phi; 138 dqv[1]=dqv[1]*phi; 139 dqv[2]=dqv[2]*phi; 137 140 138 141 return 0; … … 141 144 142 145 double compute_froude_number(double uh, 143 144 145 146 146 double h, 147 double g, 148 double epsilon) { 149 147 150 // Compute Froude number; v/sqrt(gh) 148 151 // FIXME (Ole): Not currently in use … … 166 169 // This is used by flux functions 167 170 // Input parameters uh and h may be modified by this function. 171 172 // FIXME: Perhaps inline this and profile 168 173 double _compute_speed(double *uh, 169 170 171 174 double *h, 175 double epsilon, 176 double h0) { 172 177 173 178 double u; … … 188 193 } 189 194 195 196 // Optimised squareroot computation 197 // 198 //See 199 //http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf 200 //and http://mail.opensolaris.org/pipermail/toolsgcc/2005August/000047.html 201 float fast_squareroot_approximation(float number) { 202 float x; 203 const float f = 1.5; 204 205 // Allow i and y to occupy the same memory 206 union 207 { 208 long i; 209 float y; 210 } u; 211 212 // Good initial guess 213 x = number * 0.5; 214 u.y = number; 215 u.i = 0x5f3759df  ( u.i >> 1 ); 216 217 // Take a few iterations 218 u.y = u.y * ( f  ( x * u.y * u.y ) ); 219 u.y = u.y * ( f  ( x * u.y * u.y ) ); 220 221 return number * u.y; 222 } 223 224 225 226 // Optimised squareroot computation (double version, slower) 227 double Xfast_squareroot_approximation(double number) { 228 double x; 229 const double f = 1.5; 230 231 // Allow i and y to occupy the same memory 232 union 233 { 234 long long i; 235 double y; 236 } u; 237 238 // Good initial guess 239 x = number * 0.5; 240 u.y = number; 241 u.i = 0x5fe6ec85e7de30daLL  ( u.i >> 1 ); 242 243 // Take a few iterations 244 u.y = u.y * ( f  ( x * u.y * u.y ) ); 245 u.y = u.y * ( f  ( x * u.y * u.y ) ); 246 247 return number * u.y; 248 } 249 250 190 251 // Innermost flux function (using stage w=z+h) 191 252 int _flux_function_central(double *q_left, double *q_right, 192 double z_left, double z_right, 193 double n1, double n2, 194 double epsilon, double H0, double g, 195 double *edgeflux, double *max_speed) { 253 double z_left, double z_right, 254 double n1, double n2, 255 double epsilon, double H0, double g, 256 double *edgeflux, double *max_speed) 257 { 196 258 197 259 /*Compute fluxes between volumes for the shallow water wave equation … … 212 274 double v_left, v_right; 213 275 double s_min, s_max, soundspeed_left, soundspeed_right; 214 double denom, z;276 double denom, inverse_denominator, z; 215 277 216 278 // Workspace (allocate once, use many) … … 219 281 double h0 = H0*H0; // This ensures a good balance when h approaches H0. 220 282 // But evidence suggests that h0 can be as little as 221 283 // epsilon! 222 284 223 285 // Copy conserved quantities to protect from modification 224 for (i=0; i<3; i++) { 225 q_left_rotated[i] = q_left[i]; 226 q_right_rotated[i] = q_right[i]; 227 } 286 q_left_rotated[0] = q_left[0]; 287 q_right_rotated[0] = q_right[0]; 288 q_left_rotated[1] = q_left[1]; 289 q_right_rotated[1] = q_right[1]; 290 q_left_rotated[2] = q_left[2]; 291 q_right_rotated[2] = q_right[2]; 228 292 229 293 // Align x and ymomentum with xaxis … … 231 295 _rotate(q_right_rotated, n1, n2); 232 296 233 z = (z_left+z_right)/2; // Average elevation values.234 // Even though this will nominally allow for discontinuities235 236 297 z = 0.5*(z_left + z_right); // Average elevation values. 298 // Even though this will nominally allow for discontinuities 299 // in the elevation data, there is currently no numerical 300 // support for this so results may be strange near jumps in the bed. 237 301 238 302 // Compute speeds in xdirection 239 303 w_left = q_left_rotated[0]; 240 h_left = w_left z;304 h_left = w_left  z; 241 305 uh_left = q_left_rotated[1]; 242 306 u_left = _compute_speed(&uh_left, &h_left, epsilon, h0); 243 307 244 308 w_right = q_right_rotated[0]; 245 h_right = w_right z;309 h_right = w_right  z; 246 310 uh_right = q_right_rotated[1]; 247 311 u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); … … 257 321 // Maximal and minimal wave speeds 258 322 soundspeed_left = sqrt(g*h_left); 259 soundspeed_right = sqrt(g*h_right); 260 261 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right); 262 if (s_max < 0.0) s_max = 0.0; 263 264 s_min = min(u_leftsoundspeed_left, u_rightsoundspeed_right); 265 if (s_min > 0.0) s_min = 0.0; 266 323 soundspeed_right = sqrt(g*h_right); 324 325 // Code to use fast square root optimisation if desired. 326 //soundspeed_left = fast_squareroot_approximation(g*h_left); 327 //soundspeed_right = fast_squareroot_approximation(g*h_right); 328 329 s_max = max(u_left + soundspeed_left, u_right + soundspeed_right); 330 if (s_max < 0.0) 331 { 332 s_max = 0.0; 333 } 334 335 s_min = min(u_left  soundspeed_left, u_right  soundspeed_right); 336 if (s_min > 0.0) 337 { 338 s_min = 0.0; 339 } 267 340 268 341 // Flux formulas … … 275 348 flux_right[2] = u_right*vh_right; 276 349 277 278 350 // Flux computation 279 denom = s_maxs_min; 280 if (denom < epsilon) { // FIXME (Ole): Try using H0 here 281 for (i=0; i<3; i++) edgeflux[i] = 0.0; 351 denom = s_max  s_min; 352 if (denom < epsilon) 353 { // FIXME (Ole): Try using H0 here 354 memset(edgeflux, 0, 3*sizeof(double)); 282 355 *max_speed = 0.0; 283 } else { 284 for (i=0; i<3; i++) { 356 } 357 else 358 { 359 inverse_denominator = 1.0/denom; 360 for (i = 0; i < 3; i++) 361 { 285 362 edgeflux[i] = s_max*flux_left[i]  s_min*flux_right[i]; 286 edgeflux[i] += s_max*s_min*(q_right_rotated[i] q_left_rotated[i]);287 edgeflux[i] /= denom;363 edgeflux[i] += s_max*s_min*(q_right_rotated[i]  q_left_rotated[i]); 364 edgeflux[i] *= inverse_denominator; 288 365 } 289 366 … … 315 392 // Computational function for flux computation (using stage w=z+h) 316 393 int flux_function_kinetic(double *q_left, double *q_right, 317 318 319 320 394 double z_left, double z_right, 395 double n1, double n2, 396 double epsilon, double H0, double g, 397 double *edgeflux, double *max_speed) { 321 398 322 399 /*Compute fluxes between volumes for the shallow water wave equation … … 413 490 414 491 void _manning_friction(double g, double eps, int N, 415 416 417 492 double* w, double* z, 493 double* uh, double* vh, 494 double* eta, double* xmom, double* ymom) { 418 495 419 496 int k; … … 426 503 S = g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k])); 427 504 S /= pow(h, 7.0/3); //Expensive (on Ole's home computer) 428 //S /= exp( 7.0/3.0*log(h)); //seems to save about 15% over manning_friction505 //S /= exp((7.0/3.0)*log(h)); //seems to save about 15% over manning_friction 429 506 //S /= h*h*(1 + h/3.0  h*h/9.0); //FIXME: Could use a Taylor expansion 430 507 … … 441 518 /* 442 519 void _manning_friction_explicit(double g, double eps, int N, 443 444 445 520 double* w, double* z, 521 double* uh, double* vh, 522 double* eta, double* xmom, double* ymom) { 446 523 447 524 int k; … … 452 529 h = w[k]z[k]; 453 530 if (h >= eps) { 454 455 456 457 458 459 460 461 462 531 S = g * eta[k]*eta[k] * sqrt((uh[k]*uh[k] + vh[k]*vh[k])); 532 S /= pow(h, 7.0/3); //Expensive (on Ole's home computer) 533 //S /= exp(7.0/3.0*log(h)); //seems to save about 15% over manning_friction 534 //S /= h*h*(1 + h/3.0  h*h/9.0); //FIXME: Could use a Taylor expansion 535 536 537 //Update momentum 538 xmom[k] += S*uh[k]; 539 ymom[k] += S*vh[k]; 463 540 } 464 541 } … … 471 548 472 549 double velocity_balance(double uh_i, 473 474 475 476 477 550 double uh, 551 double h_i, 552 double h, 553 double alpha, 554 double epsilon) { 478 555 // Find alpha such that speed at the vertex is within one 479 // order of magnitude of the centroid speed 556 // order of magnitude of the centroid speed 480 557 481 558 // FIXME(Ole): Work in progress … … 486 563 487 564 printf("alpha = %f, uh_i=%f, uh=%f, h_i=%f, h=%f\n", 488 565 alpha, uh_i, uh, h_i, h); 489 566 490 567 … … 500 577 if (h < epsilon) { 501 578 b = 1.0e10; // Limit 502 } else { 579 } else { 503 580 b = m*fabs(h_i  h)/h; 504 581 } … … 507 584 508 585 if (a > b) { 509 estimate = (m1)/(ab); 586 estimate = (m1)/(ab); 510 587 511 588 printf("Alpha %f, estimate %f\n", 512 513 589 alpha, estimate); 590 514 591 if (alpha < estimate) { 515 592 printf("Adjusting alpha from %f to %f\n", 516 593 alpha, estimate); 517 594 alpha = estimate; 518 595 } … … 520 597 521 598 if (h < h_i) { 522 estimate = (m1)/(ab); 599 estimate = (m1)/(ab); 523 600 524 601 printf("Alpha %f, estimate %f\n", 525 526 602 alpha, estimate); 603 527 604 if (alpha < estimate) { 528 529 530 605 printf("Adjusting alpha from %f to %f\n", 606 alpha, estimate); 607 alpha = estimate; 531 608 } 532 609 } … … 540 617 541 618 int _balance_deep_and_shallow(int N, 542 543 544 545 546 547 548 549 550 551 552 553 554 619 double* wc, 620 double* zc, 621 double* wv, 622 double* zv, 623 double* hvbar, // Retire this 624 double* xmomc, 625 double* ymomc, 626 double* xmomv, 627 double* ymomv, 628 double H0, 629 int tight_slope_limiters, 630 int use_centroid_velocities, 631 double alpha_balance) { 555 632 556 633 int k, k3, i; … … 579 656 // FIXME: Try with this one precomputed 580 657 for (i=0; i<3; i++) { 581 658 dz = max(dz, fabs(zv[k3+i]zc[k])); 582 659 } 583 660 } … … 592 669 593 670 //if (hmin < 0.0 ) { 594 // 671 // printf("hmin = %f\n",hmin); 595 672 //} 596 673 … … 607 684 608 685 if (dz > 0.0) { 609 686 alpha = max( min( alpha_balance*hmin/dz, 1.0), 0.0 ); 610 687 } else { 611 688 alpha = 1.0; // Flat bed 612 689 } 613 690 //printf("Using old style limiter\n"); … … 622 699 623 700 if (hmin < H0) { 624 625 626 627 h_diff = hc_k  hv[i];628 629 630 631 632 633 634 635 636 alpha = min(alpha, (hc_k  H0)/h_diff);637 638 639 640 641 642 643 701 alpha = 1.0; 702 for (i=0; i<3; i++) { 703 704 h_diff = hc_k  hv[i]; 705 if (h_diff <= 0) { 706 // Deep water triangle is further away from bed than 707 // shallow water (hbar < h). Any alpha will do 708 709 } else { 710 // Denominator is positive which means that we need some of the 711 // hlimited stage. 712 713 alpha = min(alpha, (hc_k  H0)/h_diff); 714 } 715 } 716 717 // Ensure alpha in [0,1] 718 if (alpha>1.0) alpha=1.0; 719 if (alpha<0.0) alpha=0.0; 720 644 721 } else { 645 646 722 // Use wlimited stage exclusively in deeper water. 723 alpha = 1.0; 647 724 } 648 725 } 649 726 650 727 651 // 728 // Let 652 729 // 653 // 654 // 730 // wvi be the wlimited stage (wvi = zvi + hvi) 731 // wvi be the hlimited state (wvi = zvi + hvi) 655 732 // 656 733 // 657 // 734 // where i=0,1,2 denotes the vertex ids 658 735 // 659 736 // Weighted balance between wlimited and hlimited stage is 660 737 // 661 // 738 // wvi := (1alpha)*(zvi+hvi) + alpha*(zvi+hvi) 662 739 // 663 740 // It follows that the updated wvi is 664 741 // wvi := zvi + (1alpha)*hvi + alpha*hvi 665 742 // 666 // 743 // Momentum is balanced between constant and limited 667 744 668 745 … … 670 747 for (i=0; i<3; i++) { 671 748 672 wv[k3+i] = zv[k3+i] + (1alpha)*hc_k + alpha*hv[i]; 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 749 wv[k3+i] = zv[k3+i] + (1alpha)*hc_k + alpha*hv[i]; 750 751 // Update momentum at vertices 752 if (use_centroid_velocities == 1) { 753 // This is a simple, efficient and robust option 754 // It uses first order approximation of velocities, but retains 755 // the order used by stage. 756 757 // Speeds at centroids 758 if (hc_k > epsilon) { 759 uc = xmomc[k]/hc_k; 760 vc = ymomc[k]/hc_k; 761 } else { 762 uc = 0.0; 763 vc = 0.0; 764 } 765 766 // Vertex momenta guaranteed to be consistent with depth guaranteeing 767 // controlled speed 768 hv[i] = wv[k3+i]  zv[k3+i]; // Recompute (balanced) vertex depth 769 xmomv[k3+i] = uc*hv[i]; 770 ymomv[k3+i] = vc*hv[i]; 771 772 } else { 773 // Update momentum as a linear combination of 774 // xmomc and ymomc (shallow) and momentum 775 // from extrapolator xmomv and ymomv (deep). 776 // This assumes that values from xmomv and ymomv have 777 // been established e.g. by the gradient limiter. 778 779 // FIXME (Ole): I think this should be used with vertex momenta 780 // computed above using centroid_velocities instead of xmomc 781 // and ymomc as they'll be more representative first order 782 // values. 783 784 xmomv[k3+i] = (1alpha)*xmomc[k] + alpha*xmomv[k3+i]; 785 ymomv[k3+i] = (1alpha)*ymomc[k] + alpha*ymomv[k3+i]; 786 787 } 711 788 } 712 789 } … … 719 796 720 797 int _protect(int N, 721 722 723 724 725 726 727 798 double minimum_allowed_height, 799 double maximum_allowed_speed, 800 double epsilon, 801 double* wc, 802 double* zc, 803 double* xmomc, 804 double* ymomc) { 728 805 729 806 int k; … … 738 815 739 816 if (hc < minimum_allowed_height) { 740 741 742 743 744 817 818 // Set momentum to zero and ensure h is non negative 819 xmomc[k] = 0.0; 820 ymomc[k] = 0.0; 821 if (hc <= 0.0) wc[k] = zc[k]; 745 822 } 746 823 } … … 757 834 758 835 if (hc <= 0.0) { 759 760 761 836 wc[k] = zc[k]; 837 xmomc[k] = 0.0; 838 ymomc[k] = 0.0; 762 839 } else { 763 840 //Reduce excessive speeds derived from division by small hc 764 765 841 //FIXME (Ole): This may be unnecessary with new slope limiters 842 //in effect. 766 843 767 844 u = xmomc[k]/hc; 768 769 770 771 //u, reduced_speed);772 773 845 if (fabs(u) > maximum_allowed_speed) { 846 reduced_speed = maximum_allowed_speed * u/fabs(u); 847 //printf("Speed (u) has been reduced from %.3f to %.3f\n", 848 // u, reduced_speed); 849 xmomc[k] = reduced_speed * hc; 850 } 774 851 775 852 v = ymomc[k]/hc; 776 777 778 779 //v, reduced_speed);780 781 853 if (fabs(v) > maximum_allowed_speed) { 854 reduced_speed = maximum_allowed_speed * v/fabs(v); 855 //printf("Speed (v) has been reduced from %.3f to %.3f\n", 856 // v, reduced_speed); 857 ymomc[k] = reduced_speed * hc; 858 } 782 859 } 783 860 } … … 790 867 791 868 int _assign_wind_field_values(int N, 792 793 794 795 796 869 double* xmom_update, 870 double* ymom_update, 871 double* s_vec, 872 double* phi_vec, 873 double cw) { 797 874 798 875 // Assign windfield values to momentum updates … … 839 916 840 917 if (!PyArg_ParseTuple(args, "OOOddOddd", 841 842 918 &normal, &ql, &qr, &zl, &zr, &edgeflux, 919 &epsilon, &g, &H0)) { 843 920 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: flux_function_central could not parse input arguments"); 844 921 return NULL; … … 847 924 848 925 _flux_function_central((double*) ql > data, 849 850 851 zr,852 853 ((double*) normal > data)[1],854 855 856 926 (double*) qr > data, 927 zl, 928 zr, 929 ((double*) normal > data)[0], 930 ((double*) normal > data)[1], 931 epsilon, H0, g, 932 (double*) edgeflux > data, 933 &max_speed); 857 934 858 935 return Py_BuildValue("d", max_speed); … … 874 951 875 952 if (!PyArg_ParseTuple(args, "dOOOOO", 876 877 953 &g, &h, &v, &x, 954 &xmom, &ymom)) { 878 955 //&epsilon)) { 879 956 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: gravity could not parse input arguments"); … … 940 1017 941 1018 if (!PyArg_ParseTuple(args, "ddOOOOOOO", 942 943 1019 &g, &eps, &w, &z, &uh, &vh, &eta, 1020 &xmom, &ymom)) { 944 1021 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: manning_friction could not parse input arguments"); 945 1022 return NULL; … … 957 1034 N = w > dimensions[0]; 958 1035 _manning_friction(g, eps, N, 959 960 961 962 963 964 965 1036 (double*) w > data, 1037 (double*) z > data, 1038 (double*) uh > data, 1039 (double*) vh > data, 1040 (double*) eta > data, 1041 (double*) xmom > data, 1042 (double*) ymom > data); 966 1043 967 1044 return Py_BuildValue(""); … … 981 1058 982 1059 if (!PyArg_ParseTuple(args, "ddOOOOOOO", 983 984 1060 &g, &eps, &w, &z, &uh, &vh, &eta, 1061 &xmom, &ymom)) 985 1062 return NULL; 986 1063 … … 996 1073 N = w > dimensions[0]; 997 1074 _manning_friction_explicit(g, eps, N, 998 999 1000 1001 1002 1003 1004 1075 (double*) w > data, 1076 (double*) z > data, 1077 (double*) uh > data, 1078 (double*) vh > data, 1079 (double*) eta > data, 1080 (double*) xmom > data, 1081 (double*) ymom > data); 1005 1082 1006 1083 return Py_BuildValue(""); … … 1012 1089 // Computational routine 1013 1090 int _extrapolate_second_order_sw(int number_of_elements, 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1091 double epsilon, 1092 double minimum_allowed_height, 1093 double beta_w, 1094 double beta_w_dry, 1095 double beta_uh, 1096 double beta_uh_dry, 1097 double beta_vh, 1098 double beta_vh_dry, 1099 long* surrogate_neighbours, 1100 long* number_of_boundaries, 1101 double* centroid_coordinates, 1102 double* stage_centroid_values, 1103 double* xmom_centroid_values, 1104 double* ymom_centroid_values, 1105 double* elevation_centroid_values, 1106 double* vertex_coordinates, 1107 double* stage_vertex_values, 1108 double* xmom_vertex_values, 1109 double* ymom_vertex_values, 1110 double* elevation_vertex_values, 1111 int optimise_dry_cells) { 1112 1113 1037 1114 1038 1115 // Local variables 1039 1116 double a, b; // Gradient vector used to calculate vertex values from centroids 1040 int k, k0,k1,k2,k3,k6,coord_index,i;1041 double x, y,x0,y0,x1,y1,x2,y2,xv0,yv0,xv1,yv1,xv2,yv2; // Vertices of the auxiliary triangle1042 double dx1, dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2;1117 int k, k0, k1, k2, k3, k6, coord_index, i; 1118 double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle 1119 double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2; 1043 1120 double dqv[3], qmin, qmax, hmin, hmax; 1044 1121 double hc, h0, h1, h2, beta_tmp, hfactor; 1045 1122 1046 1047 for (k=0; k<number_of_elements; k++) { 1123 1124 for (k = 0; k < number_of_elements; k++) 1125 { 1048 1126 k3=k*3; 1049 1127 k6=k*6; 1050 1128 1051 1052 if (number_of_boundaries[k]==3){1129 if (number_of_boundaries[k]==3) 1130 { 1053 1131 // No neighbours, set gradient on the triangle to zero 1054 1132 … … 1065 1143 continue; 1066 1144 } 1067 else { 1145 else 1146 { 1068 1147 // Triangle k has one or more neighbours. 1069 1148 // Get centroid and vertex coordinates of the triangle 1070 1149 1071 1150 // Get the vertex coordinates 1072 xv0 = vertex_coordinates[k6]; yv0=vertex_coordinates[k6+1]; 1073 xv1 = vertex_coordinates[k6+2]; yv1=vertex_coordinates[k6+3]; 1074 xv2 = vertex_coordinates[k6+4]; yv2=vertex_coordinates[k6+5]; 1151 xv0 = vertex_coordinates[k6]; 1152 yv0 = vertex_coordinates[k6+1]; 1153 xv1 = vertex_coordinates[k6+2]; 1154 yv1 = vertex_coordinates[k6+3]; 1155 xv2 = vertex_coordinates[k6+4]; 1156 yv2 = vertex_coordinates[k6+5]; 1075 1157 1076 1158 // Get the centroid coordinates 1077 coord_index =2*k;1078 x =centroid_coordinates[coord_index];1079 y =centroid_coordinates[coord_index+1];1159 coord_index = 2*k; 1160 x = centroid_coordinates[coord_index]; 1161 y = centroid_coordinates[coord_index+1]; 1080 1162 1081 1163 // Store x and y differentials for the vertices of 1082 1164 // triangle k relative to the centroid 1083 dxv0=xv0x; dxv1=xv1x; dxv2=xv2x; 1084 dyv0=yv0y; dyv1=yv1y; dyv2=yv2y; 1165 dxv0 = xv0  x; 1166 dxv1 = xv1  x; 1167 dxv2 = xv2  x; 1168 dyv0 = yv0  y; 1169 dyv1 = yv1  y; 1170 dyv2 = yv2  y; 1085 1171 } 1086 1172 1087 1173 1088 1174 1089 if (number_of_boundaries[k]<=1) {1090 1175 if (number_of_boundaries[k]<=1) 1176 { 1091 1177 //============================================== 1092 1178 // Number of boundaries <= 1 … … 1099 1185 // from this centroid and its two neighbours 1100 1186 1101 k0 =surrogate_neighbours[k3];1102 k1 =surrogate_neighbours[k3+1];1103 k2 =surrogate_neighbours[k3+2];1187 k0 = surrogate_neighbours[k3]; 1188 k1 = surrogate_neighbours[k3 + 1]; 1189 k2 = surrogate_neighbours[k3 + 2]; 1104 1190 1105 1191 // Get the auxiliary triangle's vertex coordinates 1106 1192 // (really the centroids of neighbouring triangles) 1107 coord_index =2*k0;1108 x0 =centroid_coordinates[coord_index];1109 y0 =centroid_coordinates[coord_index+1];1110 1111 coord_index =2*k1;1112 x1 =centroid_coordinates[coord_index];1113 y1 =centroid_coordinates[coord_index+1];1114 1115 coord_index =2*k2;1116 x2 =centroid_coordinates[coord_index];1117 y2 =centroid_coordinates[coord_index+1];1193 coord_index = 2*k0; 1194 x0 = centroid_coordinates[coord_index]; 1195 y0 = centroid_coordinates[coord_index+1]; 1196 1197 coord_index = 2*k1; 1198 x1 = centroid_coordinates[coord_index]; 1199 y1 = centroid_coordinates[coord_index+1]; 1200 1201 coord_index = 2*k2; 1202 x2 = centroid_coordinates[coord_index]; 1203 y2 = centroid_coordinates[coord_index+1]; 1118 1204 1119 1205 // Store x and y differentials for the vertices 1120 1206 // of the auxiliary triangle 1121 dx1=x1x0; dx2=x2x0; 1122 dy1=y1y0; dy2=y2y0; 1207 dx1 = x1  x0; 1208 dx2 = x2  x0; 1209 dy1 = y1  y0; 1210 dy2 = y2  y0; 1123 1211 1124 1212 // Calculate 2*area of the auxiliary triangle … … 1128 1216 // If the mesh is 'weird' near the boundary, 1129 1217 // the triangle might be flat or clockwise: 1130 if (area2<=0) { 1131 PyErr_SetString(PyExc_RuntimeError, 1132 "shallow_water_ext.c: negative triangle area encountered"); 1133 return 1; 1218 if (area2 <= 0) 1219 { 1220 PyErr_SetString(PyExc_RuntimeError, 1221 "shallow_water_ext.c: negative triangle area encountered"); 1222 1223 return 1; 1134 1224 } 1135 1225 … … 1139 1229 h1 = stage_centroid_values[k1]  elevation_centroid_values[k1]; 1140 1230 h2 = stage_centroid_values[k2]  elevation_centroid_values[k2]; 1141 hmin = min(min(h0, min(h1,h2)),hc);1231 hmin = min(min(h0, min(h1, h2)), hc); 1142 1232 //hfactor = hc/(hc + 1.0); 1143 1233 1144 1234 hfactor = 0.0; 1145 if (hmin > 0.001 ) { 1146 hfactor = (hmin0.001)/(hmin+0.004); 1147 } 1148 1149 if (optimise_dry_cells) { 1150 // Check if linear reconstruction is necessary for triangle k 1151 // This check will exclude dry cells. 1152 1153 hmax = max(h0,max(h1,h2)); 1154 if (hmax < epsilon) { 1155 continue; 1156 } 1157 } 1158 1159 1235 if (hmin > 0.001) 1236 { 1237 hfactor = (hmin  0.001)/(hmin + 0.004); 1238 } 1239 1240 if (optimise_dry_cells) 1241 { 1242 // Check if linear reconstruction is necessary for triangle k 1243 // This check will exclude dry cells. 1244 1245 hmax = max(h0, max(h1, h2)); 1246 if (hmax < epsilon) 1247 { 1248 continue; 1249 } 1250 } 1251 1160 1252 // 1161 1253 // stage … … 1164 1256 // Calculate the difference between vertex 0 of the auxiliary 1165 1257 // triangle and the centroid of triangle k 1166 dq0 =stage_centroid_values[k0]stage_centroid_values[k];1258 dq0 = stage_centroid_values[k0]  stage_centroid_values[k]; 1167 1259 1168 1260 // Calculate differentials between the vertices 1169 1261 // of the auxiliary triangle (centroids of neighbouring triangles) 1170 dq1=stage_centroid_values[k1]stage_centroid_values[k0]; 1171 dq2=stage_centroid_values[k2]stage_centroid_values[k0]; 1172 1262 dq1 = stage_centroid_values[k1]  stage_centroid_values[k0]; 1263 dq2 = stage_centroid_values[k2]  stage_centroid_values[k0]; 1264 1265 inv_area2 = 1.0/area2; 1173 1266 // Calculate the gradient of stage on the auxiliary triangle 1174 1267 a = dy2*dq1  dy1*dq2; 1175 a /=area2;1268 a *= inv_area2; 1176 1269 b = dx1*dq2  dx2*dq1; 1177 b /=area2;1270 b *= inv_area2; 1178 1271 1179 1272 // Calculate provisional jumps in stage from the centroid 1180 1273 // of triangle k to its vertices, to be limited 1181 dqv[0] =a*dxv0+b*dyv0;1182 dqv[1] =a*dxv1+b*dyv1;1183 dqv[2] =a*dxv2+b*dyv2;1274 dqv[0] = a*dxv0 + b*dyv0; 1275 dqv[1] = a*dxv1 + b*dyv1; 1276 dqv[2] = a*dxv2 + b*dyv2; 1184 1277 1185 1278 // Now we want to find min and max of the centroid and the 1186 1279 // vertices of the auxiliary triangle and compute jumps 1187 1280 // from the centroid to the min and max 1188 find_qmin_and_qmax(dq0, dq1,dq2,&qmin,&qmax);1281 find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax); 1189 1282 1190 1283 // Playing with dry wet interface … … 1193 1286 //if (hmin>minimum_allowed_height) 1194 1287 beta_tmp = beta_w_dry + (beta_w  beta_w_dry) * hfactor; 1195 1288 1196 1289 //printf("min_alled_height = %f\n",minimum_allowed_height); 1197 1290 //printf("hmin = %f\n",hmin); … … 1199 1292 //printf("beta_tmp = %f\n",beta_tmp); 1200 1293 // Limit the gradient 1201 limit_gradient(dqv,qmin,qmax,beta_tmp); 1202 1203 for (i=0;i<3;i++) 1204 stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i]; 1294 limit_gradient(dqv, qmin, qmax, beta_tmp); 1295 1296 //for (i=0;i<3;i++) 1297 stage_vertex_values[k3+0] = stage_centroid_values[k] + dqv[0]; 1298 stage_vertex_values[k3+1] = stage_centroid_values[k] + dqv[1]; 1299 stage_vertex_values[k3+2] = stage_centroid_values[k] + dqv[2]; 1205 1300 1206 1301 … … 1211 1306 // Calculate the difference between vertex 0 of the auxiliary 1212 1307 // triangle and the centroid of triangle k 1213 dq0 =xmom_centroid_values[k0]xmom_centroid_values[k];1308 dq0 = xmom_centroid_values[k0]  xmom_centroid_values[k]; 1214 1309 1215 1310 // Calculate differentials between the vertices 1216 1311 // of the auxiliary triangle 1217 dq1 =xmom_centroid_values[k1]xmom_centroid_values[k0];1218 dq2 =xmom_centroid_values[k2]xmom_centroid_values[k0];1312 dq1 = xmom_centroid_values[k1]  xmom_centroid_values[k0]; 1313 dq2 = xmom_centroid_values[k2]  xmom_centroid_values[k0]; 1219 1314 1220 1315 // Calculate the gradient of xmom on the auxiliary triangle 1221 1316 a = dy2*dq1  dy1*dq2; 1222 a /=area2;1317 a *= inv_area2; 1223 1318 b = dx1*dq2  dx2*dq1; 1224 b /=area2;1319 b *= inv_area2; 1225 1320 1226 1321 // Calculate provisional jumps in stage from the centroid 1227 1322 // of triangle k to its vertices, to be limited 1228 dqv[0] =a*dxv0+b*dyv0;1229 dqv[1] =a*dxv1+b*dyv1;1230 dqv[2] =a*dxv2+b*dyv2;1323 dqv[0] = a*dxv0+b*dyv0; 1324 dqv[1] = a*dxv1+b*dyv1; 1325 dqv[2] = a*dxv2+b*dyv2; 1231 1326 1232 1327 // Now we want to find min and max of the centroid and the 1233 1328 // vertices of the auxiliary triangle and compute jumps 1234 1329 // from the centroid to the min and max 1235 find_qmin_and_qmax(dq0, dq1,dq2,&qmin,&qmax);1330 find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax); 1236 1331 //beta_tmp = beta_uh; 1237 1332 //if (hmin<minimum_allowed_height) … … 1240 1335 1241 1336 // Limit the gradient 1242 limit_gradient(dqv,qmin,qmax,beta_tmp); 1243 1244 for (i=0;i<3;i++) 1245 xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i]; 1246 1337 limit_gradient(dqv, qmin, qmax, beta_tmp); 1338 1339 for (i=0; i < 3; i++) 1340 { 1341 xmom_vertex_values[k3+i] = xmom_centroid_values[k] + dqv[i]; 1342 } 1247 1343 1248 1344 // … … 1252 1348 // Calculate the difference between vertex 0 of the auxiliary 1253 1349 // triangle and the centroid of triangle k 1254 dq0 =ymom_centroid_values[k0]ymom_centroid_values[k];1350 dq0 = ymom_centroid_values[k0]  ymom_centroid_values[k]; 1255 1351 1256 1352 // Calculate differentials between the vertices 1257 1353 // of the auxiliary triangle 1258 dq1 =ymom_centroid_values[k1]ymom_centroid_values[k0];1259 dq2 =ymom_centroid_values[k2]ymom_centroid_values[k0];1354 dq1 = ymom_centroid_values[k1]  ymom_centroid_values[k0]; 1355 dq2 = ymom_centroid_values[k2]  ymom_centroid_values[k0]; 1260 1356 1261 1357 // Calculate the gradient of xmom on the auxiliary triangle 1262 1358 a = dy2*dq1  dy1*dq2; 1263 a /=area2;1359 a *= inv_area2; 1264 1360 b = dx1*dq2  dx2*dq1; 1265 b /=area2;1361 b *= inv_area2; 1266 1362 1267 1363 // Calculate provisional jumps in stage from the centroid 1268 1364 // of triangle k to its vertices, to be limited 1269 dqv[0] =a*dxv0+b*dyv0;1270 dqv[1] =a*dxv1+b*dyv1;1271 dqv[2] =a*dxv2+b*dyv2;1365 dqv[0] = a*dxv0 + b*dyv0; 1366 dqv[1] = a*dxv1 + b*dyv1; 1367 dqv[2] = a*dxv2 + b*dyv2; 1272 1368 1273 1369 // Now we want to find min and max of the centroid and the 1274 1370 // vertices of the auxiliary triangle and compute jumps 1275 1371 // from the centroid to the min and max 1276 find_qmin_and_qmax(dq0, dq1,dq2,&qmin,&qmax);1372 find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax); 1277 1373 1278 1374 //beta_tmp = beta_vh; … … 1280 1376 //if (hmin<minimum_allowed_height) 1281 1377 //beta_tmp = beta_vh_dry; 1282 beta_tmp = beta_vh_dry + (beta_vh  beta_vh_dry) * hfactor; 1378 beta_tmp = beta_vh_dry + (beta_vh  beta_vh_dry) * hfactor; 1283 1379 1284 1380 // Limit the gradient 1285 limit_gradient(dqv, qmin,qmax,beta_tmp);1381 limit_gradient(dqv, qmin, qmax, beta_tmp); 1286 1382 1287 1383 for (i=0;i<3;i++) 1288 ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i]; 1289 1384 { 1385 ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i]; 1386 } 1290 1387 } // End number_of_boundaries <=1 1291 else{ 1388 else 1389 { 1292 1390 1293 1391 //============================================== … … 1298 1396 1299 1397 // Find the only internal neighbour (k1?) 1300 for (k2=k3;k2<k3+3;k2++){ 1301 // Find internal neighbour of triangle k 1302 // k2 indexes the edges of triangle k 1303 1304 if (surrogate_neighbours[k2]!=k) 1305 break; 1306 } 1307 1308 if ((k2==k3+3)) { 1309 // If we didn't find an internal neighbour 1310 PyErr_SetString(PyExc_RuntimeError, 1311 "shallow_water_ext.c: Internal neighbour not found"); 1312 return 1; 1313 } 1314 1315 k1=surrogate_neighbours[k2]; 1398 for (k2 = k3; k2 < k3 + 3; k2++) 1399 { 1400 // Find internal neighbour of triangle k 1401 // k2 indexes the edges of triangle k 1402 1403 if (surrogate_neighbours[k2] != k) 1404 { 1405 break; 1406 } 1407 } 1408 1409 if ((k2 == k3 + 3)) 1410 { 1411 // If we didn't find an internal neighbour 1412 PyErr_SetString(PyExc_RuntimeError, 1413 "shallow_water_ext.c: Internal neighbour not found"); 1414 return 1; 1415 } 1416 1417 k1 = surrogate_neighbours[k2]; 1316 1418 1317 1419 // The coordinates of the triangle are already (x,y). 1318 1420 // Get centroid of the neighbour (x1,y1) 1319 coord_index =2*k1;1320 x1 =centroid_coordinates[coord_index];1321 y1 =centroid_coordinates[coord_index+1];1421 coord_index = 2*k1; 1422 x1 = centroid_coordinates[coord_index]; 1423 y1 = centroid_coordinates[coord_index + 1]; 1322 1424 1323 1425 // Compute x and y distances between the centroid of 1324 1426 // triangle k and that of its neighbour 1325 dx1=x1x; dy1=y1y; 1427 dx1 = x1  x; 1428 dy1 = y1  y; 1326 1429 1327 1430 // Set area2 as the square of the distance 1328 area2 =dx1*dx1+dy1*dy1;1431 area2 = dx1*dx1 + dy1*dy1; 1329 1432 1330 1433 // Set dx2=(x1x0)/((x1x0)^2+(y1y0)^2) … … 1332 1435 // respectively correspond to the x and y gradients 1333 1436 // of the conserved quantities 1334 dx2 =1.0/area2;1335 dy2 =dx2*dy1;1336 dx2 *=dx1;1437 dx2 = 1.0/area2; 1438 dy2 = dx2*dy1; 1439 dx2 *= dx1; 1337 1440 1338 1441 … … 1342 1445 1343 1446 // Compute differentials 1344 dq1 =stage_centroid_values[k1]stage_centroid_values[k];1447 dq1 = stage_centroid_values[k1]  stage_centroid_values[k]; 1345 1448 1346 1449 // Calculate the gradient between the centroid of triangle k 1347 1450 // and that of its neighbour 1348 a =dq1*dx2;1349 b =dq1*dy2;1451 a = dq1*dx2; 1452 b = dq1*dy2; 1350 1453 1351 1454 // Calculate provisional vertex jumps, to be limited 1352 dqv[0] =a*dxv0+b*dyv0;1353 dqv[1] =a*dxv1+b*dyv1;1354 dqv[2] =a*dxv2+b*dyv2;1455 dqv[0] = a*dxv0 + b*dyv0; 1456 dqv[1] = a*dxv1 + b*dyv1; 1457 dqv[2] = a*dxv2 + b*dyv2; 1355 1458 1356 1459 // Now limit the jumps 1357 if (dq1>=0.0){ 1358 qmin=0.0; 1359 qmax=dq1; 1360 } 1361 else{ 1362 qmin=dq1; 1363 qmax=0.0; 1460 if (dq1>=0.0) 1461 { 1462 qmin=0.0; 1463 qmax=dq1; 1464 } 1465 else 1466 { 1467 qmin = dq1; 1468 qmax = 0.0; 1364 1469 } 1365 1470 1366 1471 // Limit the gradient 1367 limit_gradient(dqv,qmin,qmax,beta_w); 1368 1369 for (i=0;i<3;i++) 1370 stage_vertex_values[k3+i]=stage_centroid_values[k]+dqv[i]; 1371 1472 limit_gradient(dqv, qmin, qmax, beta_w); 1473 1474 //for (i=0; i < 3; i++) 1475 //{ 1476 stage_vertex_values[k3] = stage_centroid_values[k] + dqv[0]; 1477 stage_vertex_values[k3 + 1] = stage_centroid_values[k] + dqv[1]; 1478 stage_vertex_values[k3 + 2] = stage_centroid_values[k] + dqv[2]; 1479 //} 1372 1480 1373 1481 // … … 1376 1484 1377 1485 // Compute differentials 1378 dq1 =xmom_centroid_values[k1]xmom_centroid_values[k];1486 dq1 = xmom_centroid_values[k1]  xmom_centroid_values[k]; 1379 1487 1380 1488 // Calculate the gradient between the centroid of triangle k 1381 1489 // and that of its neighbour 1382 a =dq1*dx2;1383 b =dq1*dy2;1490 a = dq1*dx2; 1491 b = dq1*dy2; 1384 1492 1385 1493 // Calculate provisional vertex jumps, to be limited 1386 dqv[0] =a*dxv0+b*dyv0;1387 dqv[1] =a*dxv1+b*dyv1;1388 dqv[2] =a*dxv2+b*dyv2;1494 dqv[0] = a*dxv0+b*dyv0; 1495 dqv[1] = a*dxv1+b*dyv1; 1496 dqv[2] = a*dxv2+b*dyv2; 1389 1497 1390 1498 // Now limit the jumps 1391 if (dq1>=0.0){ 1392 qmin=0.0; 1393 qmax=dq1; 1394 } 1395 else{ 1396 qmin=dq1; 1397 qmax=0.0; 1499 if (dq1 >= 0.0) 1500 { 1501 qmin = 0.0; 1502 qmax = dq1; 1503 } 1504 else 1505 { 1506 qmin = dq1; 1507 qmax = 0.0; 1398 1508 } 1399 1509 1400 1510 // Limit the gradient 1401 limit_gradient(dqv,qmin,qmax,beta_w); 1402 1403 for (i=0;i<3;i++) 1404 xmom_vertex_values[k3+i]=xmom_centroid_values[k]+dqv[i]; 1405 1511 limit_gradient(dqv, qmin, qmax, beta_w); 1512 1513 //for (i=0;i<3;i++) 1514 //xmom_vertex_values[k3] = xmom_centroid_values[k] + dqv[0]; 1515 //xmom_vertex_values[k3 + 1] = xmom_centroid_values[k] + dqv[1]; 1516 //xmom_vertex_values[k3 + 2] = xmom_centroid_values[k] + dqv[2]; 1517 1518 for (i = 0; i < 3;i++) 1519 { 1520 xmom_vertex_values[k3 + i] = xmom_centroid_values[k] + dqv[i]; 1521 } 1406 1522 1407 1523 // … … 1410 1526 1411 1527 // Compute differentials 1412 dq1 =ymom_centroid_values[k1]ymom_centroid_values[k];1528 dq1 = ymom_centroid_values[k1]  ymom_centroid_values[k]; 1413 1529 1414 1530 // Calculate the gradient between the centroid of triangle k 1415 1531 // and that of its neighbour 1416 a =dq1*dx2;1417 b =dq1*dy2;1532 a = dq1*dx2; 1533 b = dq1*dy2; 1418 1534 1419 1535 // Calculate provisional vertex jumps, to be limited 1420 dqv[0] =a*dxv0+b*dyv0;1421 dqv[1] =a*dxv1+b*dyv1;1422 dqv[2] =a*dxv2+b*dyv2;1536 dqv[0] = a*dxv0 + b*dyv0; 1537 dqv[1] = a*dxv1 + b*dyv1; 1538 dqv[2] = a*dxv2 + b*dyv2; 1423 1539 1424 1540 // Now limit the jumps 1425 if (dq1>=0.0){ 1426 qmin=0.0; 1427 qmax=dq1; 1428 } 1429 else{ 1430 qmin=dq1; 1431 qmax=0.0; 1541 if (dq1>=0.0) 1542 { 1543 qmin = 0.0; 1544 qmax = dq1; 1545 } 1546 else 1547 { 1548 qmin = dq1; 1549 qmax = 0.0; 1432 1550 } 1433 1551 1434 1552 // Limit the gradient 1435 limit_gradient(dqv,qmin,qmax,beta_w); 1436 1437 for (i=0;i<3;i++) 1438 ymom_vertex_values[k3+i]=ymom_centroid_values[k]+dqv[i]; 1439 1553 limit_gradient(dqv, qmin, qmax, beta_w); 1554 1555 //for (i=0;i<3;i++) 1556 //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0]; 1557 //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1]; 1558 //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2]; 1559 1560 for (i=0;i<3;i++) 1561 { 1562 ymom_vertex_values[k3 + i] = ymom_centroid_values[k] + dqv[i]; 1563 } 1564 //ymom_vertex_values[k3] = ymom_centroid_values[k] + dqv[0]; 1565 //ymom_vertex_values[k3 + 1] = ymom_centroid_values[k] + dqv[1]; 1566 //ymom_vertex_values[k3 + 2] = ymom_centroid_values[k] + dqv[2]; 1440 1567 } // else [number_of_boundaries==2] 1441 1568 } // for k=0 to number_of_elements1 1442 1569 1443 1570 return 0; 1444 } 1571 } 1445 1572 1446 1573 … … 1477 1604 Post conditions: 1478 1605 The vertices of each triangle have values from a 1479 1480 1606 limited linear reconstruction 1607 based on centroid values 1481 1608 1482 1609 */ … … 1505 1632 // Convert Python arguments to C 1506 1633 if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOi", 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 &optimise_dry_cells)) { 1521 1634 &domain, 1635 &surrogate_neighbours, 1636 &number_of_boundaries, 1637 ¢roid_coordinates, 1638 &stage_centroid_values, 1639 &xmom_centroid_values, 1640 &ymom_centroid_values, 1641 &elevation_centroid_values, 1642 &vertex_coordinates, 1643 &stage_vertex_values, 1644 &xmom_vertex_values, 1645 &ymom_vertex_values, 1646 &elevation_vertex_values, 1647 &optimise_dry_cells)) { 1648 1522 1649 PyErr_SetString(PyExc_RuntimeError, 1523 1650 "Input arguments to extrapolate_second_order_sw failed"); 1524 1651 return NULL; 1525 1652 } … … 1557 1684 // Call underlying computational routine 1558 1685 e = _extrapolate_second_order_sw(number_of_elements, 1559 1560 1561 1562 1563 1564 1565 1566 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1686 epsilon, 1687 minimum_allowed_height, 1688 beta_w, 1689 beta_w_dry, 1690 beta_uh, 1691 beta_uh_dry, 1692 beta_vh, 1693 beta_vh_dry, 1694 (long*) surrogate_neighbours > data, 1695 (long*) number_of_boundaries > data, 1696 (double*) centroid_coordinates > data, 1697 (double*) stage_centroid_values > data, 1698 (double*) xmom_centroid_values > data, 1699 (double*) ymom_centroid_values > data, 1700 (double*) elevation_centroid_values > data, 1701 (double*) vertex_coordinates > data, 1702 (double*) stage_vertex_values > data, 1703 (double*) xmom_vertex_values > data, 1704 (double*) ymom_vertex_values > data, 1705 (double*) elevation_vertex_values > data, 1706 optimise_dry_cells); 1580 1707 if (e == 1) { 1581 1708 // Use error string set inside computational routine 1582 1709 return NULL; 1583 } 1710 } 1584 1711 1585 1712 … … 1609 1736 // Convert Python arguments to C 1610 1737 if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOi", argnames, 1611 1738 &Q, &Normal, &direction)) { 1612 1739 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: rotate could not parse input arguments"); 1613 1740 return NULL; … … 1654 1781 // Computational function for flux computation 1655 1782 double _compute_fluxes_central(int number_of_elements, 1656 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671 1672 1673 1674 1675 1676 1677 1678 1679 int optimise_dry_cells) { 1680 1783 double timestep, 1784 double epsilon, 1785 double H0, 1786 double g, 1787 long* neighbours, 1788 long* neighbour_edges, 1789 double* normals, 1790 double* edgelengths, 1791 double* radii, 1792 double* areas, 1793 long* tri_full_flag, 1794 double* stage_edge_values, 1795 double* xmom_edge_values, 1796 double* ymom_edge_values, 1797 double* bed_edge_values, 1798 double* stage_boundary_values, 1799 double* xmom_boundary_values, 1800 double* ymom_boundary_values, 1801 double* stage_explicit_update, 1802 double* xmom_explicit_update, 1803 double* ymom_explicit_update, 1804 long* already_computed_flux, 1805 double* max_speed_array, 1806 int optimise_dry_cells) 1807 { 1681 1808 // Local variables 1682 double max_speed, length, area, zl, zr;1809 double max_speed, length, inv_area, zl, zr; 1683 1810 int k, i, m, n; 1684 1811 int ki, nm=0, ki2; // Index shorthands … … 1687 1814 double ql[3], qr[3], edgeflux[3]; // Work array for summing up fluxes 1688 1815 1689 static long call=1; // Static local variable flagging already computed flux 1690 1691 1816 static long call = 1; // Static local variable flagging already computed flux 1817 1692 1818 // Start computation 1693 1819 call++; // Flag 'id' of flux calculation for this timestep 1694 1820 1695 1821 // Set explicit_update to zero for all conserved_quantities. 1696 1822 // This assumes compute_fluxes called before forcing terms 1697 for (k=0; k<number_of_elements; k++) { 1698 stage_explicit_update[k]=0.0; 1699 xmom_explicit_update[k]=0.0; 1700 ymom_explicit_update[k]=0.0; 1701 } 1702 1823 memset((char*) stage_explicit_update, 0, number_of_elements*sizeof(double)); 1824 memset((char*) xmom_explicit_update, 0, number_of_elements*sizeof(double)); 1825 memset((char*) ymom_explicit_update, 0, number_of_elements*sizeof(double)); 1826 1703 1827 // For all triangles 1704 for (k =0; k<number_of_elements; k++) {1705 1828 for (k = 0; k < number_of_elements; k++) 1829 { 1706 1830 // Loop through neighbours and compute edge flux for each 1707 for (i=0; i<3; i++) { 1708 ki = k*3+i; // Linear index (triangle k, edge i) 1831 for (i = 0; i < 3; i++) 1832 { 1833 ki = k*3 + i; // Linear index to edge i of triangle k 1709 1834 1710 1835 if (already_computed_flux[ki] == call) 1836 { 1711 1837 // We've already computed the flux across this edge 1712 continue; 1713 1714 1838 continue; 1839 } 1840 1841 // Get left hand side values from triangle k, edge i 1715 1842 ql[0] = stage_edge_values[ki]; 1716 1843 ql[1] = xmom_edge_values[ki]; … … 1718 1845 zl = bed_edge_values[ki]; 1719 1846 1720 // Quantities at neighbour on nearest face 1847 // Get right hand side values either from neighbouring triangle 1848 // or from boundary array (Quantities at neighbour on nearest face). 1721 1849 n = neighbours[ki]; 1722 if (n < 0) { 1850 if (n < 0) 1851 { 1723 1852 // Neighbour is a boundary condition 1724 m = n1; // Convert negative flag to boundary index 1725 1726 qr[0] = stage_boundary_values[m]; 1727 qr[1] = xmom_boundary_values[m]; 1728 qr[2] = ymom_boundary_values[m]; 1729 zr = zl; // Extend bed elevation to boundary 1730 } else { 1731 // Neighbour is a real element 1732 m = neighbour_edges[ki]; 1733 nm = n*3+m; // Linear index (triangle n, edge m) 1734 1735 qr[0] = stage_edge_values[nm]; 1736 qr[1] = xmom_edge_values[nm]; 1737 qr[2] = ymom_edge_values[nm]; 1738 zr = bed_edge_values[nm]; 1739 } 1740 1741 1742 if (optimise_dry_cells) { 1743 // Check if flux calculation is necessary across this edge 1744 // This check will exclude dry cells. 1745 // This will also optimise cases where zl != zr as 1746 // long as both are dry 1747 1748 if ( fabs(ql[0]  zl) < epsilon && 1749 fabs(qr[0]  zr) < epsilon ) { 1750 // Cell boundary is dry 1751 1752 already_computed_flux[ki] = call; // #k Done 1753 if (n>=0) 1754 already_computed_flux[nm] = call; // #n Done 1755 1756 max_speed = 0.0; 1757 continue; 1758 } 1759 } 1760 1853 m = n  1; // Convert negative flag to boundary index 1854 1855 qr[0] = stage_boundary_values[m]; 1856 qr[1] = xmom_boundary_values[m]; 1857 qr[2] = ymom_boundary_values[m]; 1858 zr = zl; // Extend bed elevation to boundary 1859 } 1860 else 1861 { 1862 // Neighbour is a real triangle 1863 m = neighbour_edges[ki]; 1864 nm = n*3 + m; // Linear index (triangle n, edge m) 1865 1866 qr[0] = stage_edge_values[nm]; 1867 qr[1] = xmom_edge_values[nm]; 1868 qr[2] = ymom_edge_values[nm]; 1869 zr = bed_edge_values[nm]; 1870 } 1871 1872 // Now we have values for this edge  both from left and right side. 1873 1874 if (optimise_dry_cells) 1875 { 1876 // Check if flux calculation is necessary across this edge 1877 // This check will exclude dry cells. 1878 // This will also optimise cases where zl != zr as 1879 // long as both are dry 1880 1881 if (fabs(ql[0]  zl) < epsilon && 1882 fabs(qr[0]  zr) < epsilon) 1883 { 1884 // Cell boundary is dry 1885 1886 already_computed_flux[ki] = call; // #k Done 1887 if (n >= 0) 1888 { 1889 already_computed_flux[nm] = call; // #n Done 1890 } 1891 1892 max_speed = 0.0; 1893 continue; 1894 } 1895 } 1761 1896 1762 1897 // Outward pointing normal vector (domain.normals[k, 2*i:2*i+2]) … … 1765 1900 // Edge flux computation (triangle k, edge i) 1766 1901 _flux_function_central(ql, qr, zl, zr, 1767 1768 1769 1902 normals[ki2], normals[ki2+1], 1903 epsilon, H0, g, 1904 edgeflux, &max_speed); 1770 1905 1771 1906 … … 1786 1921 1787 1922 // Update neighbour n with same flux but reversed sign 1788 if (n >=0) {1789 stage_explicit_update[n] += edgeflux[0]; 1790 xmom_explicit_update[n] += edgeflux[1];1791 ymom_explicit_update[n] += edgeflux[2];1792 1793 already_computed_flux[nm] = call; // #n Done 1794 }1795 1923 if (n >= 0) 1924 { 1925 stage_explicit_update[n] += edgeflux[0]; 1926 xmom_explicit_update[n] += edgeflux[1]; 1927 ymom_explicit_update[n] += edgeflux[2]; 1928 1929 already_computed_flux[nm] = call; // #n Done 1930 } 1796 1931 1797 1932 // Update timestep based on edge i and possibly neighbour n 1798 if (tri_full_flag[k] == 1) { 1799 if (max_speed > epsilon) { 1800 1801 // Apply CFL condition for triangles joining this edge (triangle k and triangle n) 1802 1803 // CFL for triangle k 1804 timestep = min(timestep, radii[k]/max_speed); 1805 1806 if (n>=0) 1807 // Apply CFL condition for neigbour n (which is on the ith edge of triangle k) 1808 timestep = min(timestep, radii[n]/max_speed); 1809 1810 // Ted Rigby's suggested less conservative version 1811 //if (n>=0) { 1812 // timestep = min(timestep, (radii[k]+radii[n])/max_speed); 1813 //} else { 1814 // timestep = min(timestep, radii[k]/max_speed); 1815 // } 1816 } 1933 if (tri_full_flag[k] == 1) 1934 { 1935 if (max_speed > epsilon) 1936 { 1937 // Apply CFL condition for triangles joining this edge (triangle k and triangle n) 1938 1939 // CFL for triangle k 1940 timestep = min(timestep, radii[k]/max_speed); 1941 1942 if (n >= 0) 1943 { 1944 // Apply CFL condition for neigbour n (which is on the ith edge of triangle k) 1945 timestep = min(timestep, radii[n]/max_speed); 1946 } 1947 1948 // Ted Rigby's suggested less conservative version 1949 //if (n>=0) { 1950 // timestep = min(timestep, (radii[k]+radii[n])/max_speed); 1951 //} else { 1952 // timestep = min(timestep, radii[k]/max_speed); 1953 // } 1954 } 1817 1955 } 1818 1956 … … 1822 1960 // Normalise triangle k by area and store for when all conserved 1823 1961 // quantities get updated 1824 area =areas[k];1825 stage_explicit_update[k] /=area;1826 xmom_explicit_update[k] /=area;1827 ymom_explicit_update[k] /=area;1962 inv_area = 1.0/areas[k]; 1963 stage_explicit_update[k] *= inv_area; 1964 xmom_explicit_update[k] *= inv_area; 1965 ymom_explicit_update[k] *= inv_area; 1828 1966 1829 1967 … … 1832 1970 1833 1971 } // End triangle k 1834 1835 1836 1972 1837 1973 return timestep; 1838 1974 } … … 1869 2005 1870 2006 PyObject 1871 1872 1873 1874 1875 2007 *domain, 2008 *stage, 2009 *xmom, 2010 *ymom, 2011 *bed; 1876 2012 1877 2013 PyArrayObject 1878 1879 1880 1881 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 2014 *neighbours, 2015 *neighbour_edges, 2016 *normals, 2017 *edgelengths, 2018 *radii, 2019 *areas, 2020 *tri_full_flag, 2021 *stage_edge_values, 2022 *xmom_edge_values, 2023 *ymom_edge_values, 2024 *bed_edge_values, 2025 *stage_boundary_values, 2026 *xmom_boundary_values, 2027 *ymom_boundary_values, 2028 *stage_explicit_update, 2029 *xmom_explicit_update, 2030 *ymom_explicit_update, 2031 *already_computed_flux, //Tracks whether the flux across an edge has already been computed 2032 *max_speed_array; //Keeps track of max speeds for each triangle 1897 2033 1898 2034 … … 1902 2038 // Convert Python arguments to C 1903 2039 if (!PyArg_ParseTuple(args, "dOOOO", ×tep, &domain, &stage, &xmom, &ymom, &bed )) { 1904 1905 2040 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed"); 2041 return NULL; 1906 2042 } 1907 2043 … … 1940 2076 // the explicit update arrays 1941 2077 timestep = _compute_fluxes_central(number_of_elements, 1942 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959 1960 1961 1962 1963 1964 1965 2078 timestep, 2079 epsilon, 2080 H0, 2081 g, 2082 (long*) neighbours > data, 2083 (long*) neighbour_edges > data, 2084 (double*) normals > data, 2085 (double*) edgelengths > data, 2086 (double*) radii > data, 2087 (double*) areas > data, 2088 (long*) tri_full_flag > data, 2089 (double*) stage_edge_values > data, 2090 (double*) xmom_edge_values > data, 2091 (double*) ymom_edge_values > data, 2092 (double*) bed_edge_values > data, 2093 (double*) stage_boundary_values > data, 2094 (double*) xmom_boundary_values > data, 2095 (double*) ymom_boundary_values > data, 2096 (double*) stage_explicit_update > data, 2097 (double*) xmom_explicit_update > data, 2098 (double*) ymom_explicit_update > data, 2099 (long*) already_computed_flux > data, 2100 (double*) max_speed_array > data, 2101 optimise_dry_cells); 1966 2102 1967 2103 Py_DECREF(neighbours); … … 2013 2149 domain.timestep = compute_fluxes(timestep, 2014 2150 domain.epsilon, 2015 2151 domain.H0, 2016 2152 domain.g, 2017 2153 domain.neighbours, … … 2033 2169 Ymom.explicit_update, 2034 2170 already_computed_flux, 2035 optimise_dry_cells)2171 optimise_dry_cells) 2036 2172 2037 2173 … … 2066 2202 // Convert Python arguments to C 2067 2203 if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOOOi", 2068 2069 2070 2071 2072 2073 2074 2075 2076 2077 2078 2079 2080 2081 2082 2083 2084 2085 2086 2087 2088 2089 2204 ×tep, 2205 &epsilon, 2206 &H0, 2207 &g, 2208 &neighbours, 2209 &neighbour_edges, 2210 &normals, 2211 &edgelengths, &radii, &areas, 2212 &tri_full_flag, 2213 &stage_edge_values, 2214 &xmom_edge_values, 2215 &ymom_edge_values, 2216 &bed_edge_values, 2217 &stage_boundary_values, 2218 &xmom_boundary_values, 2219 &ymom_boundary_values, 2220 &stage_explicit_update, 2221 &xmom_explicit_update, 2222 &ymom_explicit_update, 2223 &already_computed_flux, 2224 &max_speed_array, 2225 &optimise_dry_cells)) { 2090 2226 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed"); 2091 2227 return NULL; … … 2118 2254 // the explicit update arrays 2119 2255 timestep = _compute_fluxes_central(number_of_elements, 2120 2121 2122 2123 2124 2125 2126 2127 2128 2129 2130 2131 2132 2133 2134 2135 2136 2137 2138 2139 2140 2141 2142 2143 2256 timestep, 2257 epsilon, 2258 H0, 2259 g, 2260 (long*) neighbours > data, 2261 (long*) neighbour_edges > data, 2262 (double*) normals > data, 2263 (double*) edgelengths > data, 2264 (double*) radii > data, 2265 (double*) areas > data, 2266 (long*) tri_full_flag > data, 2267 (double*) stage_edge_values > data, 2268 (double*) xmom_edge_values > data, 2269 (double*) ymom_edge_values > data, 2270 (double*) bed_edge_values > data, 2271 (double*) stage_boundary_values > data, 2272 (double*) xmom_boundary_values > data, 2273 (double*) ymom_boundary_values > data, 2274 (double*) stage_explicit_update > data, 2275 (double*) xmom_explicit_update > data, 2276 (double*) ymom_explicit_update > data, 2277 (long*) already_computed_flux > data, 2278 (double*) max_speed_array > data, 2279 optimise_dry_cells); 2144 2280 2145 2281 // Return updated flux timestep … … 2228 2364 // Convert Python arguments to C 2229 2365 if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOOOOOOOO", 2230 2231 2232 2233 2234 2235 2236 2237 2238 2239 2240 2241 2242 2243 2244 2245 2246 2247 2248 2249 2366 ×tep, 2367 &epsilon, 2368 &H0, 2369 &g, 2370 &neighbours, 2371 &neighbour_edges, 2372 &normals, 2373 &edgelengths, &radii, &areas, 2374 &tri_full_flag, 2375 &stage_edge_values, 2376 &xmom_edge_values, 2377 &ymom_edge_values, 2378 &bed_edge_values, 2379 &stage_boundary_values, 2380 &xmom_boundary_values, 2381 &ymom_boundary_values, 2382 &stage_explicit_update, 2383 &xmom_explicit_update, 2384 &ymom_explicit_update, 2385 &already_computed_flux)) { 2250 2386 PyErr_SetString(PyExc_RuntimeError, "Input arguments failed"); 2251 2387 return NULL; … … 2265 2401 ki = k*3+i; 2266 2402 if (((long *) already_computed_flux>data)[ki]==call)//we've already computed the flux across this edge 2267 2403 continue; 2268 2404 ql[0] = ((double *) stage_edge_values > data)[ki]; 2269 2405 ql[1] = ((double *) xmom_edge_values > data)[ki]; … … 2274 2410 n = ((long *) neighbours > data)[ki]; 2275 2411 if (n < 0) { 2276 2277 2278 2279 2280 2412 m = n1; //Convert negative flag to index 2413 qr[0] = ((double *) stage_boundary_values > data)[m]; 2414 qr[1] = ((double *) xmom_boundary_values > data)[m]; 2415 qr[2] = ((double *) ymom_boundary_values > data)[m]; 2416 zr = zl; //Extend bed elevation to boundary 2281 2417 } else { 2282 2283 2284 2285 2286 2287 2418 m = ((long *) neighbour_edges > data)[ki]; 2419 nm = n*3+m; 2420 qr[0] = ((double *) stage_edge_values > data)[nm]; 2421 qr[1] = ((double *) xmom_edge_values > data)[nm]; 2422 qr[2] = ((double *) ymom_edge_values > data)[nm]; 2423 zr = ((double *) bed_edge_values > data)[nm]; 2288 2424 } 2289 2425 // Outward pointing normal vector … … 2294 2430 //Edge flux computation 2295 2431 flux_function_kinetic(ql, qr, zl, zr, 2296 2297 2298 2432 normal[0], normal[1], 2433 epsilon, H0, g, 2434 edgeflux, &max_speed); 2299 2435 //update triangle k 2300 2436 ((long *) already_computed_flux>data)[ki]=call; … … 2304 2440 //update the neighbour n 2305 2441 if (n>=0){ 2306 2307 2308 2309 2442 ((long *) already_computed_flux>data)[nm]=call; 2443 ((double *) stage_explicit_update > data)[n] += edgeflux[0]*((double *) edgelengths > data)[nm]; 2444 ((double *) xmom_explicit_update > data)[n] += edgeflux[1]*((double *) edgelengths > data)[nm]; 2445 ((double *) ymom_explicit_update > data)[n] += edgeflux[2]*((double *) edgelengths > data)[nm]; 2310 2446 } 2311 2447 ///for (j=0; j<3; j++) { 2312 2313 2314 2315 2316 2448 ///flux[j] = edgeflux[j]*((double *) edgelengths > data)[ki]; 2449 ///} 2450 //Update timestep 2451 //timestep = min(timestep, domain.radii[k]/max_speed) 2452 //FIXME: SR Add parameter for CFL condition 2317 2453 if ( ((long *) tri_full_flag > data)[k] == 1) { 2318 2319 2320 2321 2322 2323 2454 if (max_speed > epsilon) { 2455 timestep = min(timestep, ((double *) radii > data)[k]/max_speed); 2456 //maxspeed in flux_function is calculated as max(u+a,ua) 2457 if (n>=0) 2458 timestep = min(timestep, ((double *) radii > data)[n]/max_speed); 2459 } 2324 2460 } 2325 2461 } // end for i … … 2350 2486 // Convert Python arguments to C 2351 2487 if (!PyArg_ParseTuple(args, "dddOOOO", 2352 2353 2354 2355 2488 &minimum_allowed_height, 2489 &maximum_allowed_speed, 2490 &epsilon, 2491 &wc, &zc, &xmomc, &ymomc)) { 2356 2492 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: protect could not parse input arguments"); 2357 2493 return NULL; … … 2361 2497 2362 2498 _protect(N, 2363 2364 2365 2366 2367 2368 2369 2499 minimum_allowed_height, 2500 maximum_allowed_speed, 2501 epsilon, 2502 (double*) wc > data, 2503 (double*) zc > data, 2504 (double*) xmomc > data, 2505 (double*) ymomc > data); 2370 2506 2371 2507 return Py_BuildValue(""); … … 2409 2545 // Convert Python arguments to C 2410 2546 if (!PyArg_ParseTuple(args, "OOOOOOOOOO", 2411 2412 2413 2414 2547 &domain, 2548 &wc, &zc, 2549 &wv, &zv, &hvbar, 2550 &xmomc, &ymomc, &xmomv, &ymomv)) { 2415 2551 PyErr_SetString(PyExc_RuntimeError, 2416 2552 "shallow_water_ext.c: balance_deep_and_shallow could not parse input arguments"); 2417 2553 return NULL; 2418 2554 } 2419 2420 2555 2556 2421 2557 // FIXME (Ole): I tested this without GetAttrString and got time down 2422 2558 // marginally from 4.0s to 3.8s. Consider passing everything in … … 2463 2599 N = wc > dimensions[0]; 2464 2600 _balance_deep_and_shallow(N, 2465 2466 2467 2468 2469 2470 2471 2472 2473 2474 2475 2476 (int) use_centroid_velocities,2477 2601 (double*) wc > data, 2602 (double*) zc > data, 2603 (double*) wv > data, 2604 (double*) zv > data, 2605 (double*) hvbar > data, 2606 (double*) xmomc > data, 2607 (double*) ymomc > data, 2608 (double*) xmomv > data, 2609 (double*) ymomv > data, 2610 H0, 2611 (int) tight_slope_limiters, 2612 (int) use_centroid_velocities, 2613 alpha_balance); 2478 2614 2479 2615 … … 2503 2639 // Convert Python arguments to C 2504 2640 if (!PyArg_ParseTuple(args, "OOOOd", 2505 2506 2507 2508 2641 &xmom_update, 2642 &ymom_update, 2643 &s_vec, &phi_vec, 2644 &cw)) { 2509 2645 PyErr_SetString(PyExc_RuntimeError, "shallow_water_ext.c: assign_windfield_values could not parse input arguments"); 2510 2646 return NULL; 2511 2647 } 2512 2648 2513 2649 2514 2650 N = xmom_update > dimensions[0]; 2515 2651 2516 2652 _assign_wind_field_values(N, 2517 2518 2519 2520 2521 2653 (double*) xmom_update > data, 2654 (double*) ymom_update > data, 2655 (double*) s_vec > data, 2656 (double*) phi_vec > data, 2657 cw); 2522 2658 2523 2659 return Py_BuildValue("");
Note: See TracChangeset
for help on using the changeset viewer.