- Timestamp:
- Jun 14, 2012, 12:31:32 AM (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c
r8397 r8446 227 227 edgeflux[i] = s_max*flux_left[i] - s_min*flux_right[i]; 228 228 //if(i<2) { 229 229 edgeflux[i] += s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3); 230 230 //}else{ 231 231 // edgeflux[i] += 0.5*s_max*s_min*(q_right_rotated[i] - q_left_rotated[i] - t3); … … 307 307 int useint; 308 308 double stage_edge_lim, scale_factor_shallow, bedtop, bedbot, pressure_flux, hc, hc_n; 309 double max_bed_edgevalue[number_of_elements]; 310 //double depthtemp, max_bed_edge_value; 309 //double min_bed_edgevalue[number_of_elements]; // FIXME: Set memory explicitly 311 310 static long call = 1; // Static local variable flagging already computed flux 312 311 312 double *min_bed_edgevalue; 313 314 min_bed_edgevalue = malloc(number_of_elements*sizeof(double)); 313 315 // Start computation 314 316 call++; // Flag 'id' of flux calculation for this timestep … … 321 323 322 324 323 // Compute maximum bed edge value on each triangle 324 //for (k = 0; k < number_of_elements; k++){ 325 // max_bed_edgevalue[k] = max(bed_edge_values[3*k], 326 // max(bed_edge_values[3*k+1], bed_edge_values[3*k+2])); 327 // 328 //} 325 // Compute minimum bed edge value on each triangle 326 for (k = 0; k < number_of_elements; k++){ 327 min_bed_edgevalue[k] = min(bed_edge_values[3*k], 328 min(bed_edge_values[3*k+1], bed_edge_values[3*k+2])); 329 //min_bed_edgevalue[k] = bed_centroid_values[k]; 330 331 } 329 332 330 333 … … 402 405 // NOTE: This also means no bedslope term -- which is arguably 403 406 // appropriate, since the depth is zero at the cell centroid 404 if(( hc == 0.0)&((hc_n == 0.0)|((n<0)&(qr[0]<zl)))){ 407 if(( hc == 0.0)&(hc_n == 0.0)){ 408 already_computed_flux[ki] = call; // #k Done 409 already_computed_flux[nm] = call; // #n Done 410 max_speed = 0.0; 405 411 continue; 406 412 } … … 410 416 if(n>=0){ 411 417 if(hc==0.0){ 412 ql[0]=max(min(qr[0],stage_centroid_values[k]),zl); 413 //ql[0]=qr[0]; // max(min(qr[0],stage_centroid_values[k]),zl); 418 //ql[0]=max(min(qr[0],stage_centroid_values[k]),zl); 419 //ql[0]=max(min(qr[0],0.5*(stage_centroid_values[k]+stage_centroid_values[n])),zl); 420 ql[0]=qr[0]; // max(min(qr[0],stage_centroid_values[k]),zl); 414 421 } 415 422 if(hc_n==0.0){ 416 qr[0]=max(min(ql[0],stage_centroid_values[n]),zr); 417 //qr[0]=ql[0]; 423 //qr[0]=max(min(ql[0],stage_centroid_values[n]),zr); 424 //qr[0]=max(min(ql[0],0.5*(stage_centroid_values[n]+stage_centroid_values[k])),zr); 425 qr[0]=ql[0]; 418 426 } 419 427 }else{ … … 421 429 if((hc==0.0)){ 422 430 ql[0]=max(min(qr[0],stage_centroid_values[k]),zl); 431 //ql[0]=max(min(qr[0],ql[0]),zl); 423 432 } 424 433 } … … 432 441 epsilon, h0, limiting_threshold, g, 433 442 edgeflux, &max_speed, &pressure_flux); 434 //_flux_function_fraccaro(ql, qr, zl, zr, 435 // normals[ki2], normals[ki2 + 1], 436 // epsilon, h0, limiting_threshold, g, 437 // edgeflux, &max_speed); 438 443 444 // Prevent outflow from 'seriously' dry cells 445 if(stage_centroid_values[k]<min_bed_edgevalue[k]){ 446 if(edgeflux[0]>0.0){ 447 edgeflux[0]=0.0; 448 edgeflux[1]=0.0; 449 edgeflux[2]=0.0; 450 } 451 } 452 if(n>=0){ 453 if(stage_centroid_values[n]<min_bed_edgevalue[n]){ 454 if(edgeflux[0]<0.0){ 455 edgeflux[0]=0.0; 456 edgeflux[1]=0.0; 457 edgeflux[2]=0.0; 458 } 459 } 460 } 461 439 462 // Multiply edgeflux by edgelength 440 463 length = edgelengths[ki]; … … 443 466 edgeflux[2] *= length; 444 467 445 if(n<0){ 446 if(boundary_flux_type[m]>0){ 447 // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX 448 if(boundary_flux_type[m]==1){ 449 // Zero_mass_flux_transmissive_momentum_flux boundary, or 450 // zero_mass_flux_zero_momentum_boundary 451 // Just need to set the mass flux to zero, as the 452 // transmissive momentum is ensured by the values of 453 // qr[0], qr[1], qr[2] 454 edgeflux[0] = 0.0; 455 } 456 } 457 } 468 // GET RID OF THIS: EXPERIMENTAL 469 //if(n<0){ 470 // if(boundary_flux_type[m]>0){ 471 // // IMPOSE BOUNDARY CONDITIONS DIRECTLY ON THE FLUX 472 // if(boundary_flux_type[m]==1){ 473 // // Zero_mass_flux_transmissive_momentum_flux boundary, or 474 // // zero_mass_flux_zero_momentum_boundary 475 // // Just need to set the mass flux to zero, as the 476 // // transmissive momentum is ensured by the values of 477 // // qr[0], qr[1], qr[2] 478 // printf("Warning: WEIRD EDGEFLUX THING IS ACTING"); 479 // edgeflux[0] = 0.0; 480 // } 481 // } 482 //} 458 483 459 484 // Update triangle k with flux from edge i … … 462 487 ymom_explicit_update[k] -= edgeflux[2]; 463 488 464 // Calculate the bed slope term, using the 'divergence form' from 465 // Alessandro Valiani and Lorenzo Begnudelli, Divergence Form for Bed Slope Source Term in Shallow Water Equations 466 // Journal of Hydraulic Engineering, 2006, 132:652-665 467 if(hc>0.0e-03){ 489 // Compute bed slope term 490 if(hc>-9999.0){ 468 491 //Bedslope approx 1: 469 492 //bedslope_work = g*hc*(ql[0])*length-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; … … 486 509 }else{ 487 510 // Treat nearly dry cells 488 //bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; //511 bedslope_work =-0.5*g*max(ql[0]-zl,0.)*(ql[0]-zl)*length; // 489 512 // 490 bedslope_work = -pressure_flux*length;513 //bedslope_work = -pressure_flux*length; 491 514 xmom_explicit_update[k] -= normals[ki2]*bedslope_work; 492 515 ymom_explicit_update[k] -= normals[ki2+1]*bedslope_work; … … 506 529 ymom_explicit_update[n] += edgeflux[2]; 507 530 //Add bed slope term here 508 if(hc_n> 0.0e-03){531 if(hc_n>-9999.0){ 509 532 //if(stage_centroid_values[n] > max_bed_edgevalue[n]){ 510 533 // Bedslope approx 1: … … 528 551 }else{ 529 552 // Treat nearly dry cells 530 //bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length;531 bedslope_work = -pressure_flux*length;553 bedslope_work = -0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; //-pressure_flux*length; //-0.5*g*max(qr[0]-zr,0.)*(qr[0]-zr)*length; 554 //bedslope_work = -pressure_flux*length; 532 555 xmom_explicit_update[n] += normals[ki2]*bedslope_work; 533 556 ymom_explicit_update[n] += normals[ki2+1]*bedslope_work; … … 578 601 } // End triangle k 579 602 603 free(min_bed_edgevalue); 580 604 581 605 return timestep; … … 599 623 600 624 // This acts like minimum_allowed height, but scales with the vertical 601 // distance between the bed_centroid_value and the max bed_edge_value of every triangle. 602 double minimum_relative_height=0.05; 625 // distance between the bed_centroid_value and the max bed_edge_value of 626 // every triangle. 627 double minimum_relative_height=0.1; 603 628 604 629 … … 619 644 if (hc <= 0.0){ 620 645 // Definine the minimum bed edge value on triangle k. 621 bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k]))); 646 // Setting = minimum edge value can lead to mass conservation problems 647 //bmin = 0.5*min((zv[3*k]+zv[3*k+1]),min((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k]))); 648 //bmin =0.5*bmin + 0.5*min(zv[3*k],min(zv[3*k+1],zv[3*k+2])); 649 // Setting = minimum vertex value seems better, but might tend to be less smooth 650 bmin =min(zv[3*k],min(zv[3*k+1],zv[3*k+2])); 622 651 // Minimum allowed stage = bmin 652 // WARNING: ADDING MASS if wc[k]<bmin 653 if(wc[k]<bmin){ 654 printf("Adding mass to dry cell\n"); 655 } 623 656 wc[k] = max(wc[k], bmin) ; 657 658 // Set vertex values as well. Seems that this shouldn't be 659 // needed. However from memory this is important at the first 660 // time step, for 'dry' areas where the designated stage is 661 // less than the bed centroid value 624 662 wv[3*k] = max(wv[3*k], bmin); 625 663 wv[3*k+1] = max(wv[3*k+1], bmin); … … 683 721 // dq2=q(vertex2)-q(vertex0) 684 722 685 // This is a simple implementation -- the one below is the 686 // original version -- I wonder if it is any faster? 723 // This is a simple implementation 687 724 *qmax = max(max(dq0, max(dq0+dq1, dq0+dq2)), 0.0) ; 688 725 *qmin = min(min(dq0, min(dq0+dq1, dq0+dq2)), 0.0) ; 689 726 690 //if (dq0>=0.0){691 // if (dq1>=dq2){692 // if (dq1>=0.0)693 // *qmax=dq0+dq1;694 // else695 // *qmax=dq0;696 //697 // *qmin=dq0+dq2;698 // if (*qmin>=0.0) *qmin = 0.0;699 // }700 // else{// dq1<dq2701 // if (dq2>0)702 // *qmax=dq0+dq2;703 // else704 // *qmax=dq0;705 //706 // *qmin=dq0+dq1;707 // if (*qmin>=0.0) *qmin=0.0;708 // }709 //}710 //else{//dq0<0711 // if (dq1<=dq2){712 // if (dq1<0.0)713 // *qmin=dq0+dq1;714 // else715 // *qmin=dq0;716 //717 // *qmax=dq0+dq2;718 // if (*qmax<=0.0) *qmax=0.0;719 // }720 // else{// dq1>dq2721 // if (dq2<0.0)722 // *qmin=dq0+dq2;723 // else724 // *qmin=dq0;725 //726 // *qmax=dq0+dq1;727 // if (*qmax<=0.0) *qmax=0.0;728 // }729 //}730 727 return 0; 731 728 } … … 759 756 phi=min(r*beta_w,1.0); 760 757 //phi=1.; 761 //for (i=0;i<3;i++)762 758 dqv[0]=dqv[0]*phi; 763 759 dqv[1]=dqv[1]*phi; … … 766 762 return 0; 767 763 } 764 765 //int limit_gradient2(double *dqv, double qmin, double qmax, double beta_w, double r0scale){ 766 // // EXPERIMENTAL LIMITER, DOESN'T WORK 767 // // Given provisional jumps dqv from the FV triangle centroid to its 768 // // vertices and jumps qmin (qmax) between the centroid of the FV 769 // // triangle and the minimum (maximum) of the values at the centroid of 770 // // the FV triangle and the auxiliary triangle vertices, 771 // // calculate a multiplicative factor phi by which the provisional 772 // // vertex jumps are to be limited 773 // 774 // int i; 775 // double r=1000.0, r0=1.0, phi=1.0; 776 // static double TINY = 1.0e-100; // to avoid machine accuracy problems. 777 // // FIXME: Perhaps use the epsilon used elsewhere. 778 // 779 // for (i=0;i<3;i++){ 780 // if (dqv[i]<-TINY) 781 // r0=qmin/dqv[i]*r0scale*2.0; 782 // 783 // if (dqv[i]>TINY) 784 // r0=qmax/dqv[i]*r0scale*2.0; 785 // 786 // r=min(r0,r); 787 // } 788 // 789 // phi=min(r*beta_w,1.0); 790 // //phi=1.; 791 // //for (i=0;i<3;i++) 792 // dqv[0]=dqv[0]*phi; 793 // dqv[1]=dqv[1]*phi; 794 // dqv[2]=dqv[2]*phi; 795 // 796 // return 0; 797 //} 768 798 769 799 // Computational routine … … 789 819 double* ymom_edge_values, 790 820 double* elevation_edge_values, 821 double* stage_vertex_values, 822 double* xmom_vertex_values, 823 double* ymom_vertex_values, 824 double* elevation_vertex_values, 791 825 int optimise_dry_cells, 792 826 int extrapolate_velocity_second_order) { 793 827 794 // Note that although this routine refers to EDGE values, it can equally well795 // be applied to vertex values -- it is a modified version of such a routine.796 // Also note that momentum can be extrapolated using velocity, via797 // modifications at the start and end of the routine.798 799 828 // Local variables 800 829 double a, b; // Gradient vector used to calculate edge values from centroids … … 804 833 double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin; 805 834 double hc, h0, h1, h2, beta_tmp, hfactor; 806 double dk, dv0, dv1, dv2, de[3]; 807 //double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements] 808 //double stage_centroid_store[number_of_elements]; 835 double dk, dv0, dv1, dv2, de[3], demin, dcmax, r0scale; 809 836 810 837 double *xmom_centroid_store, *ymom_centroid_store, *stage_centroid_store; … … 828 855 ymom_centroid_values[k] = ymom_centroid_values[k]/dk; 829 856 830 // Experimental -- replace stage with a weighting of 'stage' and831 // 'depth', where the weights depend on the froude number **2832 // stage_centroid_store[k]=stage_centroid_values[k];833 // hfactor =(xmom_centroid_values[k]*xmom_centroid_values[k] +834 // ymom_centroid_values[k]*ymom_centroid_values[k])/20.0;835 // //(9.8*max(stage_centroid_values[k] - elevation_centroid_values[k],minimum_allowed_height));836 //stage_centroid_values[k] -= min(1.0, hfactor)*elevation_centroid_values[k];837 857 } 838 858 } … … 887 907 dyv1 = yv1 - y; 888 908 dyv2 = yv2 - y; 909 // Compute the minimum distance from the centroid to an edge 910 //demin=min(dxv0*dxv0 +dyv0*dyv0, min(dxv1*dxv1+dyv1*dyv1, dxv2*dxv2+dyv2*dyv2)); 911 //demin=sqrt(demin); 889 912 } 890 913 … … 911 934 // triangle area will be zero, and a first order extrapolation will be 912 935 // used 913 //if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){914 //k2 = k ;915 //}916 //if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){917 //k0 = k ;918 //}919 //if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){920 //k1 = k ;921 //}936 if(stage_centroid_values[k2]<=elevation_centroid_values[k2]){ 937 k2 = k ; 938 } 939 if(stage_centroid_values[k0]<=elevation_centroid_values[k0]){ 940 k0 = k ; 941 } 942 if(stage_centroid_values[k1]<=elevation_centroid_values[k1]){ 943 k1 = k ; 944 } 922 945 // Alternative approach (BRUTAL) -- if the max neighbour bed elevation is greater 923 946 // than the min neighbour stage, then we use first order extrapolation 924 bedmax = max(elevation_centroid_values[k],925 max(elevation_centroid_values[k0],926 max(elevation_centroid_values[k1], elevation_centroid_values[k2])));927 stagemin =min(stage_centroid_values[k],928 min(stage_centroid_values[k0],929 min(stage_centroid_values[k1], stage_centroid_values[k2])));930 931 if(stagemin < bedmax){932 // This will cause first order extrapolation933 k2 = k;934 k0 = k;935 k1 = k;936 }947 //bedmax = max(elevation_centroid_values[k], 948 // max(elevation_centroid_values[k0], 949 // max(elevation_centroid_values[k1], elevation_centroid_values[k2]))); 950 //stagemin = min(stage_centroid_values[k], 951 // min(stage_centroid_values[k0], 952 // min(stage_centroid_values[k1], stage_centroid_values[k2]))); 953 // 954 //if(stagemin < bedmax){ 955 // // This will cause first order extrapolation 956 // k2 = k; 957 // k0 = k; 958 // k1 = k; 959 //} 937 960 938 961 // Get the auxiliary triangle's vertex coordinates … … 949 972 x2 = centroid_coordinates[coord_index]; 950 973 y2 = centroid_coordinates[coord_index+1]; 951 974 975 // compute the maximum distance from the centroid to a neighbouring 976 // centroid 977 //dcmax=max( (x0-x)*(x0-x) + (y0-y)*(y0-y), 978 // max((x1-x)*(x1-x) + (y1-y)*(y1-y), 979 // (x2-x)*(x2-x) + (y2-y)*(y2-y))); 980 //dcmax=sqrt(dcmax); 981 //// Ratio of centroid to edge distance -- useful in attempting to adapt limiter 982 //if(dcmax>0.){ 983 // r0scale=demin/dcmax; 984 // //printf("%f \n", r0scale); 985 //}else{ 986 // r0scale=0.5; 987 //} 988 952 989 // Store x- and y- differentials for the vertices 953 990 // of the auxiliary triangle … … 970 1007 //return -1; 971 1008 972 //stage_edge_values[k3] = 0.5*(stage_centroid_values[k]+stage_centroid_values[k0]);973 //stage_edge_values[k3+1] = 0.5*(stage_centroid_values[k]+stage_centroid_values[k1]);974 //stage_edge_values[k3+2] = 0.5*(stage_centroid_values[k]+stage_centroid_values[k2]);975 1009 stage_edge_values[k3] = stage_centroid_values[k]; 976 1010 stage_edge_values[k3+1] = stage_centroid_values[k]; … … 1002 1036 { 1003 1037 hfactor = 1.0 ;//hmin/(hmin + 0.004); 1038 //hfactor=hmin/(hmin + 0.004); 1004 1039 } 1005 1040 … … 1048 1083 find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax); 1049 1084 1050 // Playing with dry wet interface1051 //hmin = qmin;1052 //beta_tmp = beta_w_dry;1053 //if (hmin>minimum_allowed_height)1054 1085 beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor; 1055 1086 1056 //printf("min_alled_height = %f\n",minimum_allowed_height);1057 //printf("hmin = %f\n",hmin);1058 //printf("beta_w = %f\n",beta_w);1059 //printf("beta_tmp = %f\n",beta_tmp);1060 1087 // Limit the gradient 1061 1088 limit_gradient(dqv, qmin, qmax, beta_tmp); 1089 //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale); 1062 1090 1063 1091 //for (i=0;i<3;i++) … … 1066 1094 stage_edge_values[k3+2] = stage_centroid_values[k] + dqv[2]; 1067 1095 1068 //if(k==300){1069 // printf("Stage centroid values: %f, v1, %f, v2, %f, v3, %f \n", stage_centroid_values[k],1070 // stage_edge_values[k3+0], stage_edge_values[k3+1], stage_edge_values[k3+2]);1071 //}1072 1073 1096 //----------------------------------- 1074 1097 // xmomentum … … 1102 1125 find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax); 1103 1126 1104 //beta_tmp = beta_uh;1105 //if (hmin<minimum_allowed_height)1106 //beta_tmp = beta_uh_dry;1107 1127 beta_tmp = beta_uh_dry + (beta_uh - beta_uh_dry) * hfactor; 1108 1109 //if(k==300){1110 // printf("dqv0: %f, dqv1, %f, dqv2, %f \n", dqv[0],dqv[1], dqv[2]) ;1111 1112 //}1113 1128 1114 1129 // Limit the gradient 1115 1130 limit_gradient(dqv, qmin, qmax, beta_tmp); 1116 1117 //if(k==300){ 1118 // printf("dqv0: %f, dqv1, %f, dqv2, %f \n", dqv[0],dqv[1], dqv[2]) ; 1119 1120 //} 1131 //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale); 1132 1121 1133 1122 1134 for (i=0; i < 3; i++) … … 1124 1136 xmom_edge_values[k3+i] = xmom_centroid_values[k] + dqv[i]; 1125 1137 } 1126 1127 //if(k==300){1128 // printf("xmom centroid values: %f, v1, %f, v2, %f, v3, %f \n", xmom_centroid_values[k],1129 // xmom_edge_values[k3+0], xmom_edge_values[k3+1], xmom_edge_values[k3+2]);1130 //}1131 1138 1132 1139 //----------------------------------- … … 1161 1168 find_qmin_and_qmax(dq0, dq1, dq2, &qmin, &qmax); 1162 1169 1163 //beta_tmp = beta_vh;1164 //1165 //if (hmin<minimum_allowed_height)1166 //beta_tmp = beta_vh_dry;1167 1170 beta_tmp = beta_vh_dry + (beta_vh - beta_vh_dry) * hfactor; 1168 1171 1169 1172 // Limit the gradient 1170 1173 limit_gradient(dqv, qmin, qmax, beta_tmp); 1174 //limit_gradient2(dqv, qmin, qmax, beta_tmp,r0scale); 1171 1175 1172 1176 for (i=0;i<3;i++) … … 1175 1179 } 1176 1180 1177 //if(k==300){1178 // printf("ymom centroid values: %f, v1, %f, v2, %f, v3, %f \n", ymom_centroid_values[k],1179 // ymom_edge_values[k3+0], ymom_edge_values[k3+1], ymom_edge_values[k3+2]);1180 //}1181 1181 } // End number_of_boundaries <=1 1182 1182 else … … 1346 1346 limit_gradient(dqv, qmin, qmax, beta_w); 1347 1347 1348 //for (i=0;i<3;i++) 1349 //ymom_edge_values[k3] = ymom_centroid_values[k] + dqv[0]; 1350 //ymom_edge_values[k3 + 1] = ymom_centroid_values[k] + dqv[1]; 1351 //ymom_edge_values[k3 + 2] = ymom_centroid_values[k] + dqv[2]; 1352 1353 for (i=0;i<3;i++) 1354 { 1355 ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i]; 1356 } 1357 //ymom_edge_values[k3] = ymom_centroid_values[k] + dqv[0]; 1358 //ymom_edge_values[k3 + 1] = ymom_centroid_values[k] + dqv[1]; 1359 //ymom_edge_values[k3 + 2] = ymom_centroid_values[k] + dqv[2]; 1348 for (i=0;i<3;i++) 1349 { 1350 ymom_edge_values[k3 + i] = ymom_centroid_values[k] + dqv[i]; 1351 } 1360 1352 } // else [number_of_boundaries==2] 1361 1353 } // for k=0 to number_of_elements-1 1362 1354 1363 if(extrapolate_velocity_second_order==1){ 1364 // Convert back from velocity to momentum 1365 1366 for (k=0; k<number_of_elements; k++){ 1367 k3=3*k; 1368 1369 // Convert energy back to stage 1370 //stage_centroid_values[k] = stage_centroid_store[k]; 1371 1372 //Convert velocity back to momenta 1355 // Compute vertex values of quantities 1356 for (k=0; k<number_of_elements; k++){ 1357 k3=3*k; 1358 1359 // Compute stage vertex values 1360 stage_vertex_values[k3] = stage_edge_values[k3+1] + stage_edge_values[k3+2] -stage_edge_values[k3] ; 1361 stage_vertex_values[k3+1] = stage_edge_values[k3] + stage_edge_values[k3+2]-stage_edge_values[k3+1]; 1362 stage_vertex_values[k3+2] = stage_edge_values[k3] + stage_edge_values[k3+1]-stage_edge_values[k3+2]; 1363 1364 // Compute xmom vertex values 1365 xmom_vertex_values[k3] = xmom_edge_values[k3+1] + xmom_edge_values[k3+2] -xmom_edge_values[k3] ; 1366 xmom_vertex_values[k3+1] = xmom_edge_values[k3] + xmom_edge_values[k3+2]-xmom_edge_values[k3+1]; 1367 xmom_vertex_values[k3+2] = xmom_edge_values[k3] + xmom_edge_values[k3+1]-xmom_edge_values[k3+2]; 1368 1369 // Compute ymom vertex values 1370 ymom_vertex_values[k3] = ymom_edge_values[k3+1] + ymom_edge_values[k3+2] -ymom_edge_values[k3] ; 1371 ymom_vertex_values[k3+1] = ymom_edge_values[k3] + ymom_edge_values[k3+2]-ymom_edge_values[k3+1]; 1372 ymom_vertex_values[k3+2] = ymom_edge_values[k3] + ymom_edge_values[k3+1]-ymom_edge_values[k3+2]; 1373 1374 // If needed, convert from velocity to momenta 1375 if(extrapolate_velocity_second_order==1){ 1376 //Convert velocity back to momenta at centroids 1373 1377 xmom_centroid_values[k] = xmom_centroid_store[k]; 1374 1378 ymom_centroid_values[k] = ymom_centroid_store[k]; 1375 1379 1380 // Re-compute momenta at edges 1376 1381 for (i=0; i<3; i++){ 1377 //stage_edge_values[k3+i] -=(xmom_edge_values[k3+i]*xmom_edge_values[k3+i]1378 // + ymom_edge_values[k3+i]*ymom_edge_values[k3+i])*0.0;1379 //stage_edge_values[k3+i] += elevation_edge_values[k3+i];1380 //de[i] = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height);1381 //for(ii=0; ii<1; ii++){1382 //hfactor =(xmom_edge_values[k3+i]*xmom_edge_values[k3+i] +1383 // ymom_edge_values[k3+i]*ymom_edge_values[k3+i])/20.0; // 20 ~= 2*g1384 //(9.8*de[i]);1385 //hfactor = 0.0; //max(hfactor/20.0 - 0.2, 0.0);1386 //stage_edge_values[k3+i] = stage_edge_values[k3+i];//+ min(hfactor, 1.0)*elevation_edge_values[k3+i];1387 //de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],1.0e-03);1388 //}1389 1390 1382 de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0); 1391 1383 xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i]; … … 1393 1385 } 1394 1386 1395 1396 } 1397 } 1387 // Re-compute momenta at vertices 1388 for (i=0; i<3; i++){ 1389 de[i] = max(stage_vertex_values[k3+i]-elevation_vertex_values[k3+i],0.0); 1390 xmom_vertex_values[k3+i]=xmom_vertex_values[k3+i]*de[i]; 1391 ymom_vertex_values[k3+i]=ymom_vertex_values[k3+i]*de[i]; 1392 } 1393 1394 } 1395 1396 1397 } 1398 1398 1399 1399 free(xmom_centroid_store); … … 1404 1404 } 1405 1405 1406 int extrapolate_from_edges_to_vertices(int number_of_elements,1407 double* stage_edge_values,1408 double* xmom_edge_values,1409 double* ymom_edge_values,1410 double* elevation_edge_values,1411 double* stage_vertex_values,1412 double* xmom_vertex_values,1413 double* ymom_vertex_values,1414 double* elevation_vertex_values,1415 int extrapolate_velocity_second_order) {1416 1417 int i, i3;1418 double v0, v1, v2;1419 1420 for(i=0; i<number_of_elements; i++){1421 i3 = 3*i;1422 stage_vertex_values[i3] = stage_edge_values[i3+1] + stage_edge_values[i3+2] -stage_edge_values[i3] ;1423 stage_vertex_values[i3+1] = stage_edge_values[i3] + stage_edge_values[i3+2]-stage_edge_values[i3+1];1424 stage_vertex_values[i3+2] = stage_edge_values[i3] + stage_edge_values[i3+1]-stage_edge_values[i3+2];1425 1426 if(extrapolate_velocity_second_order==1){1427 // Compute x-velocity, carefully to avoid division by zero1428 //v0 = xmom_edge_values[i3+0]/(stage_edge_values[i3+0]-elevation_edge_values[i3+0]);1429 //v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);1430 //v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);1431 if(stage_edge_values[i3]>elevation_edge_values[i3]){1432 v0 = xmom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);1433 }else{1434 v0 = 0.0;1435 }1436 if(stage_edge_values[i3+1]>elevation_edge_values[i3+1]){1437 v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);1438 }else{1439 v1 = 0.0;1440 }1441 1442 if(stage_edge_values[i3+2]>elevation_edge_values[i3+2]){1443 v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);1444 }else{1445 v2 = 0.0;1446 }1447 //v1 = xmom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);1448 //v2 = xmom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);1449 1450 xmom_vertex_values[i3] = (-v0 + v1 +v2)*(stage_vertex_values[i3]-elevation_vertex_values[i3]);1451 xmom_vertex_values[i3+1] = (-v1 + v0 +v2)*(stage_vertex_values[i3+1]-elevation_vertex_values[i3+1]);1452 xmom_vertex_values[i3+2] = (-v2 + v0 +v1)*(stage_vertex_values[i3+2]-elevation_vertex_values[i3+2]);1453 1454 // Compute y-velocity, carefully to avoid division by zero1455 //v0 = ymom_edge_values[i3+0]/(stage_edge_values[i3+0]-elevation_edge_values[i3+0]);1456 //v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);1457 //v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);1458 if(stage_edge_values[i3]>elevation_edge_values[i3]){1459 v0 = ymom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);1460 }else{1461 v0 = 0.0;1462 }1463 if(stage_edge_values[i3+1]>elevation_edge_values[i3+1]){1464 v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);1465 }else{1466 v1 = 0.0;1467 }1468 1469 if(stage_edge_values[i3+2]>elevation_edge_values[i3+2]){1470 v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);1471 }else{1472 v2 = 0.0;1473 }1474 //v0 = ymom_edge_values[i3]/(stage_edge_values[i3]-elevation_edge_values[i3]);1475 //v1 = ymom_edge_values[i3+1]/(stage_edge_values[i3+1]-elevation_edge_values[i3+1]);1476 //v2 = ymom_edge_values[i3+2]/(stage_edge_values[i3+2]-elevation_edge_values[i3+2]);1477 1478 ymom_vertex_values[i3] = (-v0 + v1 +v2)*(stage_vertex_values[i3]-elevation_vertex_values[i3]);1479 ymom_vertex_values[i3+1] = (-v1 + v0 +v2)*(stage_vertex_values[i3+1]-elevation_vertex_values[i3+1]);1480 ymom_vertex_values[i3+2] = (-v2 + v0 +v1)*(stage_vertex_values[i3+2]-elevation_vertex_values[i3+2]);1481 1482 }else{1483 1484 xmom_vertex_values[i3] = xmom_edge_values[i3+1] + xmom_edge_values[i3+2]-xmom_edge_values[i3] ;1485 xmom_vertex_values[i3+1] = xmom_edge_values[i3] + xmom_edge_values[i3+2]-xmom_edge_values[i3+1];1486 xmom_vertex_values[i3+2] = xmom_edge_values[i3] + xmom_edge_values[i3+1]-xmom_edge_values[i3+2];1487 1488 ymom_vertex_values[i3] = ymom_edge_values[i3+1] + ymom_edge_values[i3+2]-ymom_edge_values[i3];1489 ymom_vertex_values[i3+1] = ymom_edge_values[i3] + ymom_edge_values[i3+2] -ymom_edge_values[i3+1];1490 ymom_vertex_values[i3+2] = ymom_edge_values[i3] + ymom_edge_values[i3+1]-ymom_edge_values[i3+2] ;1491 1492 }1493 1494 }1495 return 0 ;1496 }1497 1406 //========================================================================= 1498 1407 // Python Glue 1499 1408 //========================================================================= 1500 1409 1501 1502 //PyObject *rotate(PyObject *self, PyObject *args, PyObject *kwargs) {1503 // //1504 // // r = rotate(q, normal, direction=1)1505 // //1506 // // Where q is assumed to be a Float numeric array of length 3 and1507 // // normal a Float numeric array of length 2.1508 //1509 // // FIXME(Ole): I don't think this is used anymore1510 //1511 // PyObject *Q, *Normal;1512 // PyArrayObject *q, *r, *normal;1513 //1514 // static char *argnames[] = {"q", "normal", "direction", NULL};1515 // int dimensions[1], i, direction=1;1516 // double n1, n2;1517 //1518 // // Convert Python arguments to C1519 // if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO|i", argnames,1520 // &Q, &Normal, &direction)) {1521 // report_python_error(AT, "could not parse input arguments");1522 // return NULL;1523 // }1524 //1525 // // Input checks (convert sequences into numeric arrays)1526 // q = (PyArrayObject *)1527 // PyArray_ContiguousFromObject(Q, PyArray_DOUBLE, 0, 0);1528 // normal = (PyArrayObject *)1529 // PyArray_ContiguousFromObject(Normal, PyArray_DOUBLE, 0, 0);1530 //1531 //1532 // if (normal -> dimensions[0] != 2) {1533 // report_python_error(AT, "Normal vector must have 2 components");1534 // return NULL;1535 // }1536 //1537 // // Allocate space for return vector r (don't DECREF)1538 // dimensions[0] = 3;1539 // r = (PyArrayObject *) anuga_FromDims(1, dimensions, PyArray_DOUBLE);1540 //1541 // // Copy1542 // for (i=0; i<3; i++) {1543 // ((double *) (r -> data))[i] = ((double *) (q -> data))[i];1544 // }1545 //1546 // // Get normal and direction1547 // n1 = ((double *) normal -> data)[0];1548 // n2 = ((double *) normal -> data)[1];1549 // if (direction == -1) n2 = -n2;1550 //1551 // // Rotate1552 // _rotate((double *) r -> data, n1, n2);1553 //1554 // // Release numeric arrays1555 // Py_DECREF(q);1556 // Py_DECREF(normal);1557 //1558 // // Return result using PyArray to avoid memory leak1559 // return PyArray_Return(r);1560 //}1561 1410 1562 1411 //======================================================================== … … 1981 1830 (double*) ymom_edge_values -> data, 1982 1831 (double*) elevation_edge_values -> data, 1832 (double*) stage_vertex_values -> data, 1833 (double*) xmom_vertex_values -> data, 1834 (double*) ymom_vertex_values -> data, 1835 (double*) elevation_vertex_values -> data, 1983 1836 optimise_dry_cells, 1984 1837 extrapolate_velocity_second_order); 1985 1986 //printf("In C before edges-to-vertices");1987 1988 //Extrapolate from edges to vertices1989 e2 = extrapolate_from_edges_to_vertices(number_of_elements,1990 (double *) stage_edge_values -> data,1991 (double *) xmom_edge_values -> data,1992 (double *) ymom_edge_values -> data,1993 (double *) elevation_edge_values -> data,1994 (double *) stage_vertex_values -> data,1995 (double *) xmom_vertex_values -> data,1996 (double *) ymom_vertex_values -> data,1997 (double *) elevation_vertex_values -> data,1998 extrapolate_velocity_second_order);1999 //printf("In C after edges-to-vertices");2000 1838 2001 1839 if (e == -1) {
Note: See TracChangeset
for help on using the changeset viewer.
