Changeset 267
- Timestamp:
- Sep 2, 2004, 3:00:07 AM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/doit.py
r240 r267 1 1 import os 2 3 #Clean up 4 5 for file in os.listdir('.'): 6 if file[-1] == '~': 7 os.remove(file) 8 2 9 3 10 os.system('python compile.py') -
inundation/ga/storm_surge/pyvolution/shallow_water.py
r266 r267 511 511 hmin = min( hv[k,:] ) 512 512 513 print dz, hmin 514 513 515 #Create alpha in [0,1], where alpha==0 means using shallow 514 516 #first order scheme and alpha==1 means using the stage w as … … 552 554 # xmomc and ymomc (shallow) and momentum 553 555 # from extrapolator xmomv and ymomv (deep). 554 xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]; 555 ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]; 556 xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:] 557 ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:] 558 556 559 557 560 … … 580 583 581 584 from shallow_water_ext import balance_deep_and_shallow 582 balance_deep_and_shallow(domain.number_of_elements, 583 wc, zc, wv, zv, hv, 585 balance_deep_and_shallow(wc, zc, hc, wv, zv, hv, 584 586 xmomc, ymomc, xmomv, ymomv) 585 587 … … 985 987 gravity = gravity_c 986 988 manning_friction = manning_friction_c 987 #balance_deep_and_shallow = balance_deep_and_shallow_c989 balance_deep_and_shallow = balance_deep_and_shallow_c 988 990 989 991 #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c -
inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
r266 r267 176 176 177 177 178 /* 179 int _balance_deep_and_shallow() 180 181 #Computed linear combination between constant levels and and 182 #levels parallel to the bed elevation. 183 for k in range(domain.number_of_elements): 184 #Compute maximal variation in bed elevation 185 # This quantitiy is 186 # dz = max_i abs(z_i - z_c) 187 # and it is independent of dimension 188 # In the 1d case zc = (z0+z1)/2 189 # In the 2d case zc = (z0+z1+z2)/3 190 191 dz = max(abs(zv[k,0]-zc[k]), 192 abs(zv[k,1]-zc[k]), 193 abs(zv[k,2]-zc[k])) 194 195 196 hmin = min( hv[k,:] ) 197 198 #Create alpha in [0,1], where alpha==0 means using shallow 199 #first order scheme and alpha==1 means using the stage w as 200 #computed by the gradient limiter (1st or 2nd order) 201 # 202 #If hmin > dz/2 then alpha = 1 and the bed will have no effect 203 #If hmin < 0 then alpha = 0 reverting to constant height above bed. 178 179 int _balance_deep_and_shallow(int N, 180 double* wc, 181 double* zc, 182 double* hc, 183 double* wv, 184 double* zv, 185 double* hv, 186 double* xmomc, 187 double* ymomc, 188 double* xmomv, 189 double* ymomv) { 190 191 int k, k3, i; 192 double dz, hmin, alpha; 193 194 //Compute linear combination between constant levels and and 195 //levels parallel to the bed elevation. 196 197 for (k=0; k<N; k++) { 198 // Compute maximal variation in bed elevation 199 // This quantitiy is 200 // dz = max_i abs(z_i - z_c) 201 // and it is independent of dimension 202 // In the 1d case zc = (z0+z1)/2 203 // In the 2d case zc = (z0+z1+z2)/3 204 205 k3 = 3*k; 206 207 //FIXME: Try with this one precomputed 208 dz = 0.0; 209 hmin = hv[k3]; 210 for (i=0; i<3; i++) { 211 dz = max(dz, fabs(zv[k3+i]-zc[k])); 212 hmin = min(hmin, hv[k3+i]); 213 } 214 215 216 //Create alpha in [0,1], where alpha==0 means using shallow 217 //first order scheme and alpha==1 means using the stage w as 218 //computed by the gradient limiter (1st or 2nd order) 219 // 220 //If hmin > dz/2 then alpha = 1 and the bed will have no effect 221 //If hmin < 0 then alpha = 0 reverting to constant height above bed. 222 223 if (dz > 0.0) 224 alpha = max( min( 2*hmin/dz, 1.0), 0.0 ); 225 else 226 alpha = 1.0; //Flat bed 227 228 229 //Weighted balance between stage parallel to bed elevation 230 //(wvi = zvi + hc) and stage as computed by 1st or 2nd 231 //order gradient limiter 232 //(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids 233 // 234 //It follows that the updated wvi is 235 // wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 236 // zvi + hc + alpha*(hvi - hc) 237 // 238 //Note that hvi = zc+hc-zvi in the first order case (constant). 239 240 if (alpha < 1) { 241 for (i=0; i<3; i++) { 242 wv[k3+i] = zv[k3+i] + hc[k] + alpha*(hv[k3+i]-hc[k]); 243 204 244 205 if dz > 0.0: 206 alpha = max( min( 2*hmin/dz, 1.0), 0.0 ) 207 else: 208 #Flat bed 209 alpha = 1.0 210 211 212 #Weighted balance between stage parallel to bed elevation 213 #(wvi = zvi + hc) and stage as computed by 1st or 2nd 214 #order gradient limiter 215 #(wvi = zvi + hvi) where i=0,1,2 denotes the vertex ids 216 # 217 #It follows that the updated wvi is 218 # wvi := (1-alpha)*(zvi+hc) + alpha*(zvi+hvi) = 219 # zvi + hc + alpha*(hvi - hc) 220 # 221 #Note that hvi = zc+hc-zvi in the first order case (constant). 222 223 if alpha < 1: 224 for i in range(3): 225 wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) 226 227 228 # Update momentum as a linear combination of 229 # xmomc and ymomc (shallow) and momentum 230 # from extrapolator xmomv and ymomv (deep). 231 xmomv[k,:] = (1-alpha)*xmomc[k] + alpha*xmomv[k,:]; 232 ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]; 233 */ 234 235 245 //Update momentum as a linear combination of 246 //xmomc and ymomc (shallow) and momentum 247 //from extrapolator xmomv and ymomv (deep). 248 xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i]; 249 ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i]; 250 } 251 } 252 } 253 return 0; 254 } 255 236 256 /////////////////////////////////////////////////////////////////// 237 257 // Gateways to Python … … 542 562 // 543 563 // balance_deep_and_shallow(domain.number_of_elements, 544 // wc, zc, wv, zv, hv,564 // wc, zc, hc, wv, zv, hv, 545 565 // xmomc, ymomc, xmomv, ymomv) 546 566 … … 549 569 *wc, //Level at centroids 550 570 *zc, //Elevation at centroids 571 *hc, //Height at centroids 551 572 *wv, //Level at vertices 552 573 *zv, //Elevation at vertices … … 560 581 561 582 // Convert Python arguments to C 562 if (!PyArg_ParseTuple(args, "dOOOOOOOOO", &N, &wc, &zc, &wv, &zv, &hv, 583 if (!PyArg_ParseTuple(args, "OOOOOOOOOO", 584 &wc, &zc, &hc, 585 &wv, &zv, &hv, 563 586 &xmomc, &ymomc, &xmomv, &ymomv)) 564 587 return NULL; 565 588 566 567 /*589 N = wc -> dimensions[0]; 590 568 591 _balance_deep_and_shallow(N, 569 592 (double*) wc -> data, 570 593 (double*) zc -> data, 594 (double*) hc -> data, 571 595 (double*) wv -> data, 572 596 (double*) zv -> data, … … 576 600 (double*) xmomv -> data, 577 601 (double*) ymomv -> data); 578 */602 579 603 580 604 return Py_BuildValue("");
Note: See TracChangeset
for help on using the changeset viewer.