Changeset 266
- Timestamp:
- Sep 1, 2004, 6:31:42 PM (21 years ago)
- Location:
- inundation/ga/storm_surge/pyvolution
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/pyvolution/shallow_water.py
r263 r266 474 474 475 475 def balance_deep_and_shallow(domain): 476 """Compute linear combination between stage as computed by gradient-limiters and 477 stage computed as constant height above bed. 478 The former takes precedence when heights are large compared to the bed slope while the latter 479 takes precedence when heights are relatively small. 480 Anything in between is computed as a balanced linear combination in order to avoid numerical 481 disturbances which would otherwise appear as a result of hard switching between modes. 476 """Compute linear combination between stage as computed by 477 gradient-limiters and stage computed as constant height above bed. 478 The former takes precedence when heights are large compared to the 479 bed slope while the latter takes precedence when heights are 480 relatively small. Anything in between is computed as a balanced 481 linear combination in order to avoid numerical disturbances which 482 would otherwise appear as a result of hard switching between 483 modes. 482 484 """ 483 485 … … 553 555 ymomv[k,:] = (1-alpha)*ymomc[k] + alpha*ymomv[k,:]; 554 556 557 558 def balance_deep_and_shallow_c(domain): 559 """Wrapper for C implementation 560 """ 561 562 #Shortcuts 563 wc = domain.quantities['level'].centroid_values 564 zc = domain.quantities['elevation'].centroid_values 565 hc = wc - zc 566 567 wv = domain.quantities['level'].vertex_values 568 zv = domain.quantities['elevation'].vertex_values 569 hv = wv-zv 570 571 #Momentums at centroids 572 xmomc = domain.quantities['xmomentum'].centroid_values 573 ymomc = domain.quantities['ymomentum'].centroid_values 574 575 #Momentums at vertices 576 xmomv = domain.quantities['xmomentum'].vertex_values 577 ymomv = domain.quantities['ymomentum'].vertex_values 578 579 580 581 from shallow_water_ext import balance_deep_and_shallow 582 balance_deep_and_shallow(domain.number_of_elements, 583 wc, zc, wv, zv, hv, 584 xmomc, ymomc, xmomv, ymomv) 585 555 586 556 587 … … 954 985 gravity = gravity_c 955 986 manning_friction = manning_friction_c 956 987 #balance_deep_and_shallow = balance_deep_and_shallow_c 957 988 958 989 #distribute_to_vertices_and_edges = distribute_to_vertices_and_edges_c -
inundation/ga/storm_surge/pyvolution/shallow_water_ext.c
r260 r266 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. 204 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 178 235 179 236 /////////////////////////////////////////////////////////////////// … … 482 539 483 540 541 PyObject *balance_deep_and_shallow(PyObject *self, PyObject *args) { 542 // 543 // balance_deep_and_shallow(domain.number_of_elements, 544 // wc, zc, wv, zv, hv, 545 // xmomc, ymomc, xmomv, ymomv) 546 547 548 PyArrayObject 549 *wc, //Level at centroids 550 *zc, //Elevation at centroids 551 *wv, //Level at vertices 552 *zv, //Elevation at vertices 553 *hv, //Heights at vertices 554 *xmomc, //Momentums at centroids and vertices 555 *ymomc, 556 *xmomv, 557 *ymomv; 558 559 int N; //, err; 560 561 // Convert Python arguments to C 562 if (!PyArg_ParseTuple(args, "dOOOOOOOOO", &N, &wc, &zc, &wv, &zv, &hv, 563 &xmomc, &ymomc, &xmomv, &ymomv)) 564 return NULL; 565 566 567 /* 568 _balance_deep_and_shallow(N, 569 (double*) wc -> data, 570 (double*) zc -> data, 571 (double*) wv -> data, 572 (double*) zv -> data, 573 (double*) hv -> data, 574 (double*) xmomc -> data, 575 (double*) ymomc -> data, 576 (double*) xmomv -> data, 577 (double*) ymomv -> data); 578 */ 579 580 return Py_BuildValue(""); 581 } 484 582 485 583 … … 499 597 {"gravity", gravity, METH_VARARGS, "Print out"}, 500 598 {"manning_friction", manning_friction, METH_VARARGS, "Print out"}, 599 {"balance_deep_and_shallow", balance_deep_and_shallow, 600 METH_VARARGS, "Print out"}, 501 601 //{"distribute_to_vertices_and_edges", 502 602 // distribute_to_vertices_and_edges, METH_VARARGS},
Note: See TracChangeset
for help on using the changeset viewer.