Changeset 6694
- Timestamp:
- Apr 1, 2009, 8:41:37 PM (16 years ago)
- Location:
- anuga_work/development/anuga_1d
- Files:
-
- 1 added
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/comp_flux_ext_wellbalanced.c
r5832 r6694 3 3 #include "math.h" 4 4 #include <stdio.h> 5 #include "util_ext.h" 5 6 const double pi = 3.14159265358979; 6 7 7 8 // Shared code snippets9 #include "util_ext.h"10 11 12 //double max(double a, double b) {13 // double z;14 // z=(a>b)?a:b;15 // return z;}16 17 //double min(double a, double b) {18 // double z;19 // z=(a<b)?a:b;20 // return z;}21 22 23 24 25 26 8 //Innermost flux function (using w=z+h) 27 int _flux_function_wellbalanced(double *q_left, double *q_right, 28 double normals, double g, double epsilon, 29 double *edgeflux, double *max_speed) { 30 9 int _flux_function_wellbalanced(double *q_left, 10 double *q_right, 11 double normals, 12 double g, 13 double epsilon, 14 double *edgeflux, 15 double *max_speed) { 31 16 int i; 32 17 double z_left, z_right; 33 18 double flux_left[2], flux_right[2]; 34 19 double z_star, w_left, h_left, h_left_star, uh_left, soundspeed_left, u_left; 35 double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right; //h_right,20 double w_right, h_right, h_right_star, uh_right, soundspeed_right, u_right; 36 21 double s_max, s_min, denom; 37 22 38 w_left = q_left[0];23 w_left = q_left[0]; 39 24 uh_left = q_left[1]; 40 25 uh_left = uh_left*normals; 41 z_left = q_left[2]; 42 43 w_right = q_right[0]; 26 z_left = q_left[2]; 27 h_left = q_left[3]; 28 u_left = q_left[4]; 29 30 w_right = q_right[0]; 44 31 uh_right = q_right[1]; 45 32 uh_right = uh_right*normals; 46 z_right = q_right[2]; 33 z_right = q_right[2]; 34 h_right = q_right[3]; 35 u_right = q_right[4]; 47 36 48 37 if (z_left-z_right != 0.0 ) printf("z_l - z_r = %f \n",z_left-z_right); 49 38 50 //z = (z_left+z_right)/2.0;51 39 z_star = max(z_left, z_right); //equation(7) 52 40 53 // Compute speeds in x-direction 54 h_left = w_left-z_left; //This is used for computing the edgeflux[1]. 55 h_left_star = max(0, w_left-z_star); //(8) 56 41 // Compute speeds in x-direction. 42 h_left_star = max(0.0, w_left-z_star); //equation(8) 43 57 44 if (h_left_star < epsilon) { 58 45 u_left = 0.0; h_left_star = 0.0; … … 61 48 u_left = uh_left/h_left_star; 62 49 } 63 64 h_right = w_right-z_right; 65 h_right_star = max(0, w_right-z_star); //(9) 50 51 h_right_star = max(0.0, w_right-z_star); //equation(9) 66 52 67 53 if (h_right_star < epsilon) { … … 92 78 // Flux computation 93 79 denom = s_max-s_min; 94 if (denom < epsilon ) {80 if (denom < epsilon && z_right > z_left) { 95 81 for (i=0; i<2; i++) edgeflux[i] = 0.0; 82 // Maximal wavespeed 83 edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star ); 84 edgeflux[1] *= normals; 96 85 *max_speed = 0.0; 97 } else { 86 } else if (denom < epsilon){ 87 for (i=0; i<2; i++) edgeflux[i] = 0.0; 88 // Maximal wavespeed 89 *max_speed = 0.0; 90 }else { 98 91 edgeflux[0] = s_max*flux_left[0] - s_min*flux_right[0]; 99 92 edgeflux[0] += s_max*s_min*(h_right_star-h_left_star); //s_max*s_min*(qr[0]-ql[0]); … … 102 95 edgeflux[1] += s_max*s_min*(flux_right[0]-flux_left[0]); //s_max*s_min*(qr[1]-ql[1]); 103 96 edgeflux[1] /= denom; 97 edgeflux[1] = edgeflux[1] + (0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star ); 104 98 edgeflux[1] *= normals; 105 106 //edgeflux[1] = edgeflux[1] + 0.5*g*0.5*(h_left+h_right)*0.5*(h_left+h_right) - 0.5*g*h_left_star*h_left_star; 107 //edgeflux[1] = edgeflux[1] + 0.5*g*h_right*h_right - 0.5*g*h_right_star*h_right_star; 108 edgeflux[1] = edgeflux[1] + 0.5*g*h_left*h_left - 0.5*g*h_left_star*h_left_star; //what about the h_right and h_right_star?????????????????????? 109 // Maximal wavespeed 110 *max_speed = max(fabs(s_max), fabs(s_min)); 99 // Maximal wavespeed 100 *max_speed = max(fabs(s_max), fabs(s_min)); 111 101 } 112 //printf("max_speed = %g \n",*max_speed);102 113 103 return 0; 114 104 } … … 121 111 double epsilon, 122 112 double g, 113 123 114 long* neighbours, 124 115 long* neighbour_vertices, 125 116 double* normals, 126 117 double* areas, 118 127 119 double* stage_edge_values, 128 120 double* xmom_edge_values, 129 121 double* bed_edge_values, 122 double* height_edge_values, 123 double* vel_edge_values, 124 130 125 double* stage_boundary_values, 131 126 double* xmom_boundary_values, 127 double* bed_boundary_values, 128 double* height_boundary_values, 129 double* vel_boundary_values, 130 132 131 double* stage_explicit_update, 133 132 double* xmom_explicit_update, 133 double* bed_explicit_update, 134 double* height_explicit_update, 135 double* vel_explicit_update, 136 134 137 int number_of_elements, 135 138 double* max_speed_array) { 136 139 137 double flux[2], ql[ 3], qr[3], edgeflux[2];140 double flux[2], ql[5], qr[5], edgeflux[2]; 138 141 double max_speed, normal; 139 142 int k, i, ki, n, m, nm=0; … … 149 152 ql[1] = xmom_edge_values[ki]; 150 153 ql[2] = bed_edge_values[ki]; 151 //ql[3] = velocity_edge_values[ki];152 //ql[4] = height_edge_values[ki];154 ql[3] = height_edge_values[ki]; 155 ql[4] = vel_edge_values[ki]; 153 156 154 157 n = neighbours[ki]; 158 printf("neighbours[ki] = %li \n",neighbours[ki]); 155 159 if (n<0) { 156 160 m = -n-1; 157 161 qr[0] = stage_boundary_values[m]; 158 162 qr[1] = xmom_boundary_values[m]; 159 qr[2] = ql[2]; 163 qr[2] = ql[2]; 164 qr[3] = height_edge_values[m]; 165 qr[4] = vel_edge_values[m]; 160 166 161 167 } else { … … 164 170 qr[0] = stage_edge_values[nm]; 165 171 qr[1] = xmom_edge_values[nm]; 166 qr[2] = bed_edge_values[nm]; 172 qr[2] = bed_edge_values[nm]; 173 qr[3] = height_edge_values[nm]; 174 qr[4] = vel_edge_values[nm]; 167 175 } 168 176 … … 175 183 if (max_speed > epsilon) { 176 184 // Original CFL calculation 177 178 timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. ????????????????????????????????????????????? 185 timestep = min(timestep, 0.5*areas[k]/max_speed); //Here, CFL=1.0 is assumed. 179 186 if (n>=0) { 180 timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. ?????????????????????????????????????????????187 timestep = min(timestep, 0.5*areas[n]/max_speed); //Here, CFL=1.0 is assumed. 181 188 } 182 189 } … … 190 197 max_speed_array[k]=max_speed; 191 198 } 192 //printf("timestep = %f \n",timestep);199 193 200 return timestep; 194 201 } … … 209 216 *normals, 210 217 *areas, 218 211 219 *stage_edge_values, 212 220 *xmom_edge_values, 213 221 *bed_edge_values, 222 *height_edge_values, 223 *vel_edge_values, 224 214 225 *stage_boundary_values, 215 226 *xmom_boundary_values, 227 *bed_boundary_values, 228 *height_boundary_values, 229 *vel_boundary_values, 230 216 231 *stage_explicit_update, 217 232 *xmom_explicit_update, 233 *bed_explicit_update, 234 *height_explicit_update, 235 *vel_explicit_update, 236 218 237 *max_speed_array; 219 238 … … 222 241 223 242 // Convert Python arguments to C 224 if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOO iO",243 if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOOOOOOOOOiO", 225 244 ×tep, 226 245 &epsilon, 227 246 &g, 247 228 248 &neighbours, 229 249 &neighbour_vertices, 230 250 &normals, 231 251 &areas, 252 232 253 &stage_edge_values, 233 254 &xmom_edge_values, 234 255 &bed_edge_values, 256 &height_edge_values, 257 &vel_edge_values, 258 235 259 &stage_boundary_values, 236 260 &xmom_boundary_values, 261 &bed_boundary_values, 262 &height_boundary_values, 263 &vel_boundary_values, 264 237 265 &stage_explicit_update, 238 266 &xmom_explicit_update, 267 &bed_explicit_update, 268 &height_explicit_update, 269 &vel_explicit_update, 270 239 271 &number_of_elements, 240 272 &max_speed_array)) { 241 PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext .c: compute_fluxes_extcould not parse input");273 PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext_wellbalanced.c: compute_fluxes_ext_wellbalanced could not parse input"); 242 274 return NULL; 243 275 } … … 249 281 epsilon, 250 282 g, 283 251 284 (long*) neighbours -> data, 252 285 (long*) neighbour_vertices -> data, 253 286 (double*) normals -> data, 254 287 (double*) areas -> data, 288 255 289 (double*) stage_edge_values -> data, 256 290 (double*) xmom_edge_values -> data, 257 291 (double*) bed_edge_values -> data, 292 (double*) height_edge_values -> data, 293 (double*) vel_edge_values -> data, 294 258 295 (double*) stage_boundary_values -> data, 259 296 (double*) xmom_boundary_values -> data, 297 (double*) bed_boundary_values -> data, 298 (double*) height_boundary_values -> data, 299 (double*) vel_boundary_values -> data, 300 260 301 (double*) stage_explicit_update -> data, 261 302 (double*) xmom_explicit_update -> data, 303 (double*) bed_explicit_update -> data, 304 (double*) height_explicit_update -> data, 305 (double*) vel_explicit_update -> data, 306 262 307 number_of_elements, 263 308 (double*) max_speed_array -> data); … … 266 311 return Py_BuildValue("d", timestep); 267 312 } 268 269 270 313 271 314 … … 285 328 import_array(); 286 329 } 287 288 -
anuga_work/development/anuga_1d/config.py
r5827 r6694 109 109 110 110 #Specific to shallow water W.E. 111 minimum_allowed_height =1.0e-6 #Water depth below which it is considered to be 0111 h_min = minimum_allowed_height = 1.0e-6 #1.0e-6 #Water depth below which it is considered to be 0 112 112 maximum_allowed_speed = 1000.0 -
anuga_work/development/anuga_1d/domain.py
r6453 r6694 916 916 # Update ghosts 917 917 self.update_ghosts() 918 918 919 919 # Initial update of vertex and edge values 920 self.distribute_to_vertices_and_edges() 921 920 self.distribute_to_vertices_and_edges() 921 922 922 # Update extrema if necessary (for reporting) 923 923 self.update_extrema() … … 1502 1502 Q = self.quantities[name] 1503 1503 Q.boundary_values[i] = q[j] 1504 #print 'Q=',Q 1504 1505 1505 1506 def update_timestep(self, yieldstep, finaltime): -
anuga_work/development/anuga_1d/shallow_water_domain_suggestion2.py
r5827 r6694 53 53 54 54 conserved_quantities = ['stage', 'xmomentum'] 55 other_quantities = ['elevation', 'friction', 'height', 'velocity'] 55 evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity'] 56 other_quantities = ['friction'] 56 57 Generic_Domain.__init__(self, coordinates, boundary, 57 conserved_quantities, other_quantities,58 conserved_quantities, evolved_quantities, other_quantities, 58 59 tagged_elements) 59 60 … … 91 92 self.quantities_to_be_stored = ['stage','xmomentum'] 92 93 93 self.__doc__ = 'shallow_water_domain' 94 self.__doc__ = 'shallow_water_domain_suggestion2' 95 self.check_integrity() 94 96 95 97 … … 148 150 msg = 'Second conserved quantity must be "xmomentum"' 149 151 assert self.conserved_quantities[1] == 'xmomentum', msg 152 153 msg = 'First evolved quantity must be "stage"' 154 assert self.evolved_quantities[0] == 'stage', msg 155 msg = 'Second evolved quantity must be "xmomentum"' 156 assert self.evolved_quantities[1] == 'xmomentum', msg 157 msg = 'Third evolved quantity must be "elevation"' 158 assert self.evolved_quantities[2] == 'elevation', msg 159 msg = 'Fourth evolved quantity must be "height"' 160 assert self.evolved_quantities[3] == 'height', msg 161 msg = 'Fifth evolved quantity must be "velocity"' 162 assert self.evolved_quantities[4] == 'velocity', msg 150 163 151 164 def extrapolate_second_order_sw(self): … … 168 181 distribute_to_vertices_and_edges(self) 169 182 170 def evolve(self, yieldstep = None, finaltime = None, duration = None, 171 skip_initial_step = False): 172 """Specialisation of basic evolve method from parent class 173 """ 174 175 #Call check integrity here rather than from user scripts 176 #self.check_integrity() 177 178 #msg = 'Parameter beta_h must be in the interval [0, 1)' 179 #assert 0 <= self.beta_h < 1.0, msg 180 #msg = 'Parameter beta_w must be in the interval [0, 1)' 181 #assert 0 <= self.beta_w < 1.0, msg 182 183 184 #Initial update of vertex and edge values before any storage 185 #and or visualisation 186 187 #self.distribute_to_vertices_and_edges() ????????????????????????????????? 188 189 #Initialise real time viz if requested 190 #if self.visualise is True and self.time == 0.0: 191 # if self.visualiser is None: 192 # self.initialise_visualiser() 193 # 194 # self.visualiser.update_timer() 195 # self.visualiser.setup_all() 196 197 #Store model data, e.g. for visualisation 198 #if self.store is True and self.time == 0.0: 199 # self.initialise_storage() 200 # #print 'Storing results in ' + self.writer.filename 201 #else: 202 # pass 203 # #print 'Results will not be stored.' 204 # #print 'To store results set domain.store = True' 205 # #FIXME: Diagnostic output should be controlled by 206 # # a 'verbose' flag living in domain (or in a parent class) 207 183 def evolve(self, yieldstep = None, finaltime = None, duration = None, skip_initial_step = False): 208 184 #Call basic machinery from parent class 209 for t in Generic_Domain.evolve(self, yieldstep, finaltime,duration, 210 skip_initial_step): 211 #Real time viz 212 # if self.visualise is True: 213 # self.visualiser.update_all() 214 # self.visualiser.update_timer() 215 216 217 #Store model data, e.g. for subsequent visualisation 218 # if self.store is True: 219 # self.store_timestep(self.quantities_to_be_stored) 220 221 #FIXME: Could maybe be taken from specified list 222 #of 'store every step' quantities 223 224 #Pass control on to outer loop for more specific actions 185 for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration, skip_initial_step): 225 186 yield(t) 226 187 … … 352 313 353 314 354 #We have got h and u at vertex, then the following is the calculation of fluxes!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!315 #We have got h and u at vertex, then the following is the calculation of fluxes! 355 316 soundspeed_left = sqrt(g*h_left) 356 317 soundspeed_right = sqrt(g*h_right) … … 656 617 import sys 657 618 from Numeric import zeros, Float 658 from config import epsilon, h0659 619 660 620 N = domain.number_of_elements 661 621 timestep = float(sys.maxint) 662 #epsilon = domain.epsilon622 epsilon = domain.epsilon 663 623 g = domain.g 664 624 neighbours = domain.neighbours … … 667 627 areas = domain.areas 668 628 629 Stage = domain.quantities['stage'] 630 Xmom = domain.quantities['xmomentum'] 631 Bed = domain.quantities['elevation'] 632 Height = domain.quantities['height'] 633 Velocity = domain.quantities['velocity'] 634 635 636 stage_boundary_values = Stage.boundary_values 637 xmom_boundary_values = Xmom.boundary_values 638 bed_boundary_values = Bed.boundary_values 639 height_boundary_values= Height.boundary_values 640 vel_boundary_values = Velocity.boundary_values 641 642 stage_explicit_update = Stage.explicit_update 643 xmom_explicit_update = Xmom.explicit_update 644 bed_explicit_values = Bed.explicit_update 645 height_explicit_values= Height.explicit_update 646 vel_explicit_values = Velocity.explicit_update 647 648 max_speed_array = domain.max_speed_array 649 650 domain.distribute_to_vertices_and_edges() 651 domain.update_boundary() 652 653 stage_V = Stage.vertex_values 654 xmom_V = Xmom.vertex_values 655 bed_V = Bed.vertex_values 656 height_V= Height.vertex_values 657 vel_V = Velocity.vertex_values 658 659 number_of_elements = len(stage_V) 660 661 #from config import epsilon, h0, h_min 662 # ##Check this: 663 #for i in range(N): 664 # height_V[i] = stage_V[i] - bed_V[i] 665 # if height_V[i] <= h_min: #epsilon : 666 # height_V[i] = 0.0 667 # stage_V[i] = bed_V[i] 668 # xmom_V[i] = 0.0 669 # vel_V[i] = 0.0 670 # else: 671 # vel_V[i] = xmom_V[i]/(height_V[i] + h0/height_V[i]) 672 # ##End of part to be checked. 673 674 675 #print 'neighbours=',neighbours 676 #print 'neighbour_vertices',neighbour_vertices 677 #print 'normals=',normals 678 #print 'areas=',areas 679 from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced 680 domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep, 681 epsilon, 682 g, 683 684 neighbours, 685 neighbour_vertices, 686 normals, 687 areas, 688 689 stage_V, 690 xmom_V, 691 bed_V, 692 height_V, 693 vel_V, 694 695 stage_boundary_values, 696 xmom_boundary_values, 697 bed_boundary_values, 698 height_boundary_values, 699 vel_boundary_values, 700 701 stage_explicit_update, 702 xmom_explicit_update, 703 bed_explicit_values, 704 height_explicit_values, 705 vel_explicit_values, 706 707 number_of_elements, 708 max_speed_array) 709 710 711 # ################################### 712 713 714 715 716 717 718 # Module functions for gradient limiting (distribute_to_vertices_and_edges) 719 720 def distribute_to_vertices_and_edges(domain): 721 """Distribution from centroids to vertices specific to the 722 shallow water wave 723 equation. 724 725 It will ensure that h (w-z) is always non-negative even in the 726 presence of steep bed-slopes by taking a weighted average between shallow 727 and deep cases. 728 729 In addition, all conserved quantities get distributed as per either a 730 constant (order==1) or a piecewise linear function (order==2). 731 732 FIXME: more explanation about removal of artificial variability etc 733 734 Precondition: 735 All quantities defined at centroids and bed elevation defined at 736 vertices. 737 738 Postcondition 739 Conserved quantities defined at vertices 740 741 """ 742 743 #from config import optimised_gradient_limiter 744 745 #Remove very thin layers of water 746 #protect_against_infinitesimal_and_negative_heights(domain) 747 748 import sys 749 from Numeric import zeros, Float 750 from config import epsilon, h0, h_min 751 752 N = domain.number_of_elements 753 754 #Shortcuts 669 755 Stage = domain.quantities['stage'] 670 Xmom 671 Bed 756 Xmom = domain.quantities['xmomentum'] 757 Bed = domain.quantities['elevation'] 672 758 Height = domain.quantities['height'] 673 759 Velocity = domain.quantities['velocity'] 674 760 761 #Arrays 675 762 w_C = Stage.centroid_values 676 763 uh_C = Xmom.centroid_values 677 z_C = Bed.centroid_values 764 z_C = Bed.centroid_values 678 765 h_C = Height.centroid_values 679 766 u_C = Velocity.centroid_values … … 681 768 for i in range(N): 682 769 h_C[i] = w_C[i] - z_C[i] 683 if h_C[i] < 770 if h_C[i] <= h_min: #epsilon : 684 771 h_C[i] = 0.0 685 772 w_C[i] = z_C[i] 686 773 uh_C[i] = 0.0 687 688 for i in range(len(h_C)): 689 if h_C[i] < epsilon: 690 u_C[i] = 0.0 #Could have been negative 691 h_C[i] = 0.0 774 u_C[i] = 0.0 692 775 else: 693 694 695 for name in [' height', 'velocity']:776 u_C[i] = uh_C[i]/(h_C[i] + h0/h_C[i]) 777 778 for name in ['stage', 'height', 'velocity']: 696 779 Q = domain.quantities[name] 697 780 if domain.order == 1: … … 705 788 raise 'Unknown order' 706 789 707 w_V = domain.quantities['stage'].vertex_values 708 h_V = domain.quantities['height'].vertex_values 709 u_V = domain.quantities['velocity'].vertex_values 710 uh_V = domain.quantities['xmomentum'].vertex_values 711 z_V = w_V - h_V 712 790 stage_V = domain.quantities['stage'].vertex_values 791 bed_V = domain.quantities['elevation'].vertex_values 792 h_V = domain.quantities['height'].vertex_values 793 u_V = domain.quantities['velocity'].vertex_values 794 xmom_V = domain.quantities['xmomentum'].vertex_values 795 796 bed_V[:] = stage_V - h_V 797 xmom_V[:] = u_V * h_V 713 798 714 stage_boundary_values = Stage.boundary_values 715 xmom_boundary_values = Xmom.boundary_values 716 stage_explicit_update = Stage.explicit_update 717 xmom_explicit_update = Xmom.explicit_update 718 max_speed_array = domain.max_speed_array 719 #number_of_elements = len(w_V) 799 return 800 720 801 721 #domain.distribute_to_vertices_and_edges()722 #domain.update_boundary()723 #stage_V = Stage.vertex_values724 #xmom_V = Xmom.vertex_values725 #bed_V = Bed.vertex_values726 727 728 from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced #from comp_flux_ext import compute_fluxes_ext729 730 domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,731 epsilon,732 g,733 neighbours,734 neighbour_vertices,735 normals,736 areas,737 w_V,738 uh_V,739 z_V,740 stage_boundary_values,741 xmom_boundary_values,742 stage_explicit_update,743 xmom_explicit_update,744 N,745 max_speed_array)746 747 # ###################################748 749 750 751 752 753 754 # Module functions for gradient limiting (distribute_to_vertices_and_edges)755 756 def distribute_to_vertices_and_edges(domain):757 """Distribution from centroids to vertices specific to the758 shallow water wave759 equation.760 761 It will ensure that h (w-z) is always non-negative even in the762 presence of steep bed-slopes by taking a weighted average between shallow763 and deep cases.764 765 In addition, all conserved quantities get distributed as per either a766 constant (order==1) or a piecewise linear function (order==2).767 768 FIXME: more explanation about removal of artificial variability etc769 770 Precondition:771 All quantities defined at centroids and bed elevation defined at772 vertices.773 774 Postcondition775 Conserved quantities defined at vertices776 777 """778 779 #from config import optimised_gradient_limiter780 781 #Remove very thin layers of water782 protect_against_infinitesimal_and_negative_heights(domain)783 784 import sys785 from Numeric import zeros, Float786 787 N = domain.number_of_elements788 789 #Shortcuts790 Stage = domain.quantities['stage']791 Xmom = domain.quantities['xmomentum']792 Bed = domain.quantities['elevation']793 Height = domain.quantities['height']794 Velocity = domain.quantities['velocity']795 796 #Arrays797 w_C = Stage.centroid_values798 uh_C = Xmom.centroid_values799 z_C = Bed.centroid_values800 h_C = Height.centroid_values801 u_C = Velocity.centroid_values802 803 #print id(h_C)804 for i in range(N):805 h_C[i] = w_C[i] - z_C[i] #This is h at centroids: OK806 807 from config import epsilon, h0808 809 for i in range(len(h_C)):810 if h_C[i] < epsilon:811 u_C[i] = 0.0 #Could have been negative812 h_C[i] = 0.0813 else:814 u_C[i] = uh_C[i]/(h_C[i] + h0/h_C[i])815 816 for name in ['stage', 'velocity', 'height' ]:817 Q = domain.quantities[name]818 if domain.order == 1:819 Q.extrapolate_first_order()820 elif domain.order == 2:821 #print "add extrapolate second order to shallow water"822 #if name != 'height':823 Q.extrapolate_second_order()824 #Q.limit()825 else:826 raise 'Unknown order'827 828 stage = domain.quantities['stage'].vertex_values #This is w at vertices: OK829 h_V = domain.quantities['height'].vertex_values830 bed = stage - h_V831 u_V = domain.quantities['velocity'].vertex_values832 xmom_V = domain.quantities['xmomentum'].vertex_values833 834 #h_V[:] = stage - bed835 xmom_V[:] = u_V * h_V836 return stage, xmom_V, bed, h_V, u_V837 802 # 838 803 def protect_against_infinitesimal_and_negative_heights(domain): … … 1064 1029 #self.xmom = domain.quantities['xmomentum'].edge_values 1065 1030 #self.ymom = domain.quantities['ymomentum'].edge_values 1066 self.normals = domain.normals 1067 self.stage = domain.quantities['stage'].vertex_values 1068 self.xmom = domain.quantities['xmomentum'].vertex_values 1031 self.normals = domain.normals 1032 self.stage = domain.quantities['stage'].vertex_values 1033 self.xmom = domain.quantities['xmomentum'].vertex_values 1034 self.bed = domain.quantities['elevation'].vertex_values 1035 self.height = domain.quantities['height'].vertex_values 1036 self.velocity = domain.quantities['velocity'].vertex_values 1069 1037 1070 1038 from Numeric import zeros, Float 1071 1039 #self.conserved_quantities = zeros(3, Float) 1072 self.conserved_quantities = zeros(2, Float) 1040 #self.conserved_quantities = zeros(2, Float) 1041 self.evolved_quantities = zeros(5, Float) 1073 1042 1074 1043 def __repr__(self): 1075 1044 return 'Reflective_boundary' 1076 1045 1077 1046 """ 1078 1047 def evaluate(self, vol_id, edge_id): 1079 """Reflective boundaries reverses the outward momentum1080 of the volume they serve.1081 """1082 1083 q = self. conserved_quantities1048 #Reflective boundaries reverses the outward momentum 1049 #of the volume they serve. 1050 1051 #q = self.conserved_quantities 1052 q = self.evolved_quantities 1084 1053 q[0] = self.stage[vol_id, edge_id] 1085 1054 q[1] = self.xmom[vol_id, edge_id] … … 1101 1070 1102 1071 return q 1072 """ 1073 def evaluate(self, vol_id, edge_id): 1074 """Reflective boundaries reverses the outward momentum 1075 of the volume they serve. 1076 """ 1077 q = self.evolved_quantities 1078 q[0] = self.stage[vol_id, edge_id] 1079 q[1] = -self.xmom[vol_id, edge_id] 1080 q[2] = self.bed[vol_id, edge_id] 1081 q[3] = self.height[vol_id, edge_id] 1082 q[4] = -self.velocity[vol_id, edge_id] 1083 return q 1084 1085 1086 1103 1087 1104 1088 class Dirichlet_boundary(Boundary):
Note: See TracChangeset
for help on using the changeset viewer.