Changeset 8359
- Timestamp:
- Mar 13, 2012, 3:46:12 PM (12 years ago)
- Location:
- trunk/anuga_work/development/gareth
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_work/development/gareth/experimental/balanced_dev/swb2_domain_ext.c
r8358 r8359 291 291 // Local variables 292 292 double max_speed, length, inv_area, zl, zr; 293 double h0 = H0*H0; // H0*H0; //This ensures a good balance when h approaches H0.293 double h0 = H0*H0; // This ensures a good balance when h approaches H0. 294 294 295 295 double limiting_threshold = 10 * H0; // Avoid applying limiter below this … … 580 580 581 581 int k; 582 double hc, wmin;582 double hc, bmin, bmax; 583 583 double u, v, reduced_speed; 584 585 // This acts like minimum_allowed height, but scales with the vertical 586 // distance between the bed_centroid_value and the max bed_edge_value of every triangle. 587 double minimum_relative_height=0.05; 584 588 585 589 … … 588 592 for (k=0; k<N; k++) { 589 593 hc = wc[k] - zc[k]; 590 591 if (hc < minimum_allowed_height) { 594 // Definine the maximum bed edge value on triangle k. 595 bmax = 0.5*max((zv[3*k]+zv[3*k+1]),max((zv[3*k+1]+zv[3*k+2]),(zv[3*k+2]+zv[3*k]))); 596 597 if (hc < max(minimum_relative_height*(bmax-zc[k]), minimum_allowed_height)) { 592 598 593 599 // Set momentum to zero and ensure h is non negative 600 // NOTE: THIS IS IMPORTANT -- WE ARE SETTING MOMENTUM TO ZERO 594 601 xmomc[k] = 0.0; 595 602 ymomc[k] = 0.0; 603 596 604 if (hc <= 0.0){ 597 // Definine the minimum possible stage value as the minimum bed edge value on triangle k. 598 wmin = 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]))); 599 wc[k] = max(wc[k], wmin) ; //max(wc[k], wmin); 600 wv[3*k] = max(wv[3*k], wmin); 601 wv[3*k+1] = max(wv[3*k+1], wmin); 602 wv[3*k+2] = max(wv[3*k+2], wmin); 603 //if(wc[k]< wmin){ 604 // //printf("Limiting wc[ %d ] = %f to %f, not %f \n", k, wc[k], max(wc[k], wmin), min(wc[k], wmin)); 605 // wc[k] = wmin ; //max(wc[k], wmin); 606 // wv[3*k] = max(wv[3*k], wmin); 607 // wv[3*k+1] = max(wv[3*k+1], wmin); 608 // wv[3*k+2] = max(wv[3*k+2], wmin); 609 //} 605 // Definine the minimum bed edge value on triangle k. 606 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]))); 607 // Minimum allowed stage = bmin 608 wc[k] = max(wc[k], bmin) ; 609 wv[3*k] = max(wv[3*k], bmin); 610 wv[3*k+1] = max(wv[3*k+1], bmin); 611 wv[3*k+2] = max(wv[3*k+2], bmin); 610 612 } 611 613 } … … 780 782 // Local variables 781 783 double a, b; // Gradient vector used to calculate edge values from centroids 782 int k, k0, k1, k2, k3, k6, coord_index, i ;784 int k, k0, k1, k2, k3, k6, coord_index, i, ii; 783 785 double x, y, x0, y0, x1, y1, x2, y2, xv0, yv0, xv1, yv1, xv2, yv2; // Vertices of the auxiliary triangle 784 786 double dx1, dx2, dy1, dy2, dxv0, dxv1, dxv2, dyv0, dyv1, dyv2, dq0, dq1, dq2, area2, inv_area2; 785 787 double dqv[3], qmin, qmax, hmin, hmax, bedmax, stagemin; 786 788 double hc, h0, h1, h2, beta_tmp, hfactor; 787 double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2; 789 double xmom_centroid_store[number_of_elements], ymom_centroid_store[number_of_elements], dk, dv0, dv1, dv2, de[3]; 790 double stage_centroid_store[number_of_elements]; 788 791 789 792 if(extrapolate_velocity_second_order==1){ … … 798 801 ymom_centroid_store[k] = ymom_centroid_values[k]; 799 802 ymom_centroid_values[k] = ymom_centroid_values[k]/dk; 803 804 // Experimental -- replace stage with a weighting of 'stage' and 805 // 'depth', where the weights depend on the froude number **2 806 // stage_centroid_store[k]=stage_centroid_values[k]; 807 // hfactor =(xmom_centroid_values[k]*xmom_centroid_values[k] + 808 // ymom_centroid_values[k]*ymom_centroid_values[k])/20.0; 809 // //(9.8*max(stage_centroid_values[k] - elevation_centroid_values[k],minimum_allowed_height)); 810 //stage_centroid_values[k] -= min(1.0, hfactor)*elevation_centroid_values[k]; 800 811 } 801 812 } … … 1325 1336 for (k=0; k<number_of_elements; k++){ 1326 1337 k3=3*k; 1327 //dv0 = max(stage_edge_values[k3]-elevation_edge_values[k3],0.); 1328 //dv1 = max(stage_edge_values[k3+1]-elevation_edge_values[k3+1],0.); 1329 //dv2 = max(stage_edge_values[k3+2]-elevation_edge_values[k3+2],0.); 1330 1331 // NOTE: dv0 etc may be negative -- we don't prevent this though, 1332 // because if we set them to zero, then the relation between the 1333 // centroid and edge momenta is no longer linear. 1334 dv0 = stage_edge_values[k3]-elevation_edge_values[k3]; 1335 dv1 = stage_edge_values[k3+1]-elevation_edge_values[k3+1]; 1336 dv2 = stage_edge_values[k3+2]-elevation_edge_values[k3+2]; 1337 1338 //Correct centroid and edge values 1338 1339 // Convert energy back to stage 1340 //stage_centroid_values[k] = stage_centroid_store[k]; 1341 1342 //Convert velocity back to momenta 1339 1343 xmom_centroid_values[k] = xmom_centroid_store[k]; 1340 xmom_edge_values[k3]=xmom_edge_values[k3]*dv0;1341 xmom_edge_values[k3+1]=xmom_edge_values[k3+1]*dv1;1342 xmom_edge_values[k3+2]=xmom_edge_values[k3+2]*dv2;1343 1344 1344 ymom_centroid_values[k] = ymom_centroid_store[k]; 1345 ymom_edge_values[k3]=ymom_edge_values[k3]*dv0; 1346 ymom_edge_values[k3+1]=ymom_edge_values[k3+1]*dv1; 1347 ymom_edge_values[k3+2]=ymom_edge_values[k3+2]*dv2; 1345 1346 for (i=0; i<3; i++){ 1347 //stage_edge_values[k3+i] -=(xmom_edge_values[k3+i]*xmom_edge_values[k3+i] 1348 // + ymom_edge_values[k3+i]*ymom_edge_values[k3+i])*0.0; 1349 //stage_edge_values[k3+i] += elevation_edge_values[k3+i]; 1350 //de[i] = max(stage_centroid_values[k]-elevation_centroid_values[k],minimum_allowed_height); 1351 //for(ii=0; ii<1; ii++){ 1352 //hfactor =(xmom_edge_values[k3+i]*xmom_edge_values[k3+i] + 1353 // ymom_edge_values[k3+i]*ymom_edge_values[k3+i])/20.0; // 20 ~= 2*g 1354 //(9.8*de[i]); 1355 //hfactor = 0.0; //max(hfactor/20.0 - 0.2, 0.0); 1356 //stage_edge_values[k3+i] = stage_edge_values[k3+i];//+ min(hfactor, 1.0)*elevation_edge_values[k3+i]; 1357 //de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],1.0e-03); 1358 //} 1359 1360 de[i] = max(stage_edge_values[k3+i]-elevation_edge_values[k3+i],0.0); 1361 xmom_edge_values[k3+i]=xmom_edge_values[k3+i]*de[i]; 1362 ymom_edge_values[k3+i]=ymom_edge_values[k3+i]*de[i]; 1363 } 1364 1348 1365 1349 1366 } -
trunk/anuga_work/development/gareth/tests/channel_floodplain/plotme.py
r8355 r8359 4 4 5 5 # Time-index to plot outputs from 6 index= 11006 index=75 7 7 8 8 #p2 = util.get_output('channel_floodplain1_bal_dev.sww', minimum_allowed_height=0.01) -
trunk/anuga_work/development/gareth/tests/dam_break/plotme.py
r8355 r8359 9 9 p2_st=util.get_centroids(p_st) 10 10 11 p_dev = util.get_output('dam_break_201203 08_153245/dam_break.sww', 0.001)11 p_dev = util.get_output('dam_break_20120313_153101/dam_break.sww', 0.001) 12 12 p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True) 13 13 -
trunk/anuga_work/development/gareth/tests/runup_sinusoid/runup_sinusoid.py
r8358 r8359 18 18 #Setup computational domain 19 19 #--------- 20 points, vertices, boundary = anuga.rectangular_cross(40,40 )20 points, vertices, boundary = anuga.rectangular_cross(40,40, len1=1., len2=1.) 21 21 22 22 domain=Domain(points,vertices,boundary) # Create Domain … … 25 25 domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1}) 26 26 #domain.set_store_vertices_uniquely(True) 27 domain.minimum_allowed_height=0.001 27 28 28 29 #------------------ … … 30 31 #------------------ 31 32 scale_me=1.0 33 34 #domain.minimum_allowed_height=domain.minimum_allowed_height*scale_me # Seems needed to make the algorithms behave 35 32 36 def topography(x,y): 33 37 return (-x/2.0 +0.05*numpy.sin((x+y)*50.0))*scale_me … … 67 71 # Associate boundary tags with boundary objects 68 72 #---------------------------------------------- 69 domain.set_boundary({'left': Br, 'right': B d, 'top': Br, 'bottom':Br})73 domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br}) 70 74 71 75 #------------------------------ … … 75 79 for t in domain.evolve(yieldstep=0.1,finaltime=20.0): 76 80 print domain.timestepping_statistics() 81 xx = domain.quantities['xmomentum'].centroid_values 82 yy = domain.quantities['ymomentum'].centroid_values 83 dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values 84 dd = (dd)*(dd>1.0e-03)+1.0e-06 85 vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5 86 vv = vv*(dd>1.0e-03) 87 print 'Peak velocity is: ', vv.max(), vv.argmax() 77 88 78 89 print 'Finished' -
trunk/anuga_work/development/gareth/tests/shallow_steep_slope/channel_SU_sparse.py
r8354 r8359 57 57 58 58 Br = anuga.Reflective_boundary(domain) # Solid reflective wall 59 Bo= anuga.Dirichlet_boundary([-9.9369037,0.0,0]) # Outflow for steady uniform flow with a depth integrated velocity of 0.1 59 #Bo= anuga.Dirichlet_boundary([-9.9369037,0.0,0]) # Outflow for steady uniform flow with a depth integrated velocity of 0.1 60 61 def constant_water(t): 62 return -9.936903 63 64 #Bo= anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, constant_water) # Outflow for steady uniform flow with a depth integrated velocity of 0.1 65 Bo = anuga.Transmissive_boundary(domain) 60 66 domain.set_boundary({'left': Br, 'right': Bo, 'top': Br, 'bottom': Br}) 61 67 #------------------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.