# Changeset 5587

Ignore:
Timestamp:
Jul 30, 2008, 5:03:47 PM (15 years ago)
Message:

Added in minimum height

Location:
anuga_work/development/anuga_1d
Files:
7 edited

Unmodified
Removed
• ## anuga_work/development/anuga_1d/comp_flux_ext.c

 r5564 // Function to obtain speed from momentum and depth. // This is used by flux functions // Input parameters uh and h may be modified by this function. double _compute_speed(double *uh, double *h, double epsilon, double h0) { double u; if (*h < epsilon) { *h = 0.0;  //Could have been negative u = 0.0; } else { u = *uh/(*h + h0/ *h); } // Adjust momentum to be consistent with speed *uh = u * *h; return u; } //Innermost flux function (using w=z+h) int _flux_function(double *q_left, double *q_right, double z_left, double z_right, double normals, double g, double epsilon, double normals, double g, double epsilon, double h0, double *edgeflux, double *max_speed) { double s_max, s_min, denom; //printf("h0 = %f \n",h0); ql[0] = q_left[0]; ql[1] = q_left[1]; h_left = w_left-z; uh_left = ql[1]; if (h_left < epsilon) { u_left = 0.0; h_left = 0.0; } else { u_left = uh_left/h_left; } u_left = _compute_speed(&uh_left, &h_left, epsilon, h0); w_right = qr[0]; h_right = w_right-z; uh_right = qr[1]; if (h_right < epsilon) { u_right = 0.0; h_right = 0.0; } else { u_right = uh_right/h_right; } u_right = _compute_speed(&uh_right, &h_right, epsilon, h0); soundspeed_left = sqrt(g*h_left); // Computational function for flux computation double _compute_fluxes_ext(double timestep, double epsilon, double g, long* neighbours, long* neighbour_vertices, double* normals, double* areas, double* stage_edge_values, double* xmom_edge_values, double* bed_edge_values, double* stage_boundary_values, double* xmom_boundary_values, double* stage_explicit_update, double* xmom_explicit_update, int number_of_elements, double* max_speed_array) { double flux[2], ql[2], qr[2], edgeflux[2]; double epsilon, double g, double h0, long* neighbours, long* neighbour_vertices, double* normals, double* areas, double* stage_edge_values, double* xmom_edge_values, double* bed_edge_values, double* stage_boundary_values, double* xmom_boundary_values, double* stage_explicit_update, double* xmom_explicit_update, int number_of_elements, double* max_speed_array) { double flux[2], ql[2], qr[2], edgeflux[2]; double zl, zr, max_speed, normal; int k, i, ki, n, m, nm=0; normal = normals[ki]; _flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed); _flux_function(ql, qr, zl, zr, normal, g, epsilon, h0, edgeflux, &max_speed); flux[0] -= edgeflux[0]; flux[1] -= edgeflux[1]; *max_speed_array; double timestep, epsilon, g; double timestep, epsilon, g, h0; int number_of_elements; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "dddOOOOOOOOOOOiO", if (!PyArg_ParseTuple(args, "ddddOOOOOOOOOOOiO", ×tep, &epsilon, &g, &h0, &neighbours, &neighbour_vertices, epsilon, g, h0, (long*) neighbours -> data, (long*) neighbour_vertices -> data, *max_speed_array; double timestep, epsilon, g; double timestep, epsilon, g, h0; int number_of_elements; epsilon           = get_python_double(domain,"epsilon"); g                 = get_python_double(domain,"g"); h0                = get_python_double(domain,"h0"); // the explicit update arrays timestep = _compute_fluxes_ext(timestep, epsilon, g, (long*) neighbours -> data, (long*) neighbour_vertices -> data, (double*) normals -> data, (double*) areas -> data, (double*) stage_vertex_values -> data, (double*) xmom_vertex_values -> data, (double*) bed_vertex_values -> data, (double*) stage_boundary_values -> data, (double*) xmom_boundary_values -> data, (double*) stage_explicit_update -> data, (double*) xmom_explicit_update -> data, number_of_elements, (double*) max_speed_array -> data); epsilon, g, h0, (long*) neighbours -> data, (long*) neighbour_vertices -> data, (double*) normals -> data, (double*) areas -> data, (double*) stage_vertex_values -> data, (double*) xmom_vertex_values -> data, (double*) bed_vertex_values -> data, (double*) stage_boundary_values -> data, (double*) xmom_boundary_values -> data, (double*) stage_explicit_update -> data, (double*) xmom_explicit_update -> data, number_of_elements, (double*) max_speed_array -> data);
• ## anuga_work/development/anuga_1d/config.py

 r5535 """Module where global model parameters are set """Module where global model parameters are set for anuga_1d """ epsilon = 1.0e-12 h0 = 1.0e-6 default_boundary_tag = 'exterior' beta_w = 0.9 beta_w = 1.5 beta_h = 0.2 CFL = 1.0  #FIXME (ole): Is this in use yet??
• ## anuga_work/development/anuga_1d/dam_h_sudi.py

 r5535 import os from math import sqrt from shallow_water_h import * #from shallow_water_h import * from shallow_water_domain import * from Numeric import zeros, Float from analytic_dam_sudi import AnalyticDam L=2000.0 N=100 N=400 cell_len=L/N domain=Domain(points) domain.default_order = 1 domain.default_order = 2 domain.default_time_order = 1 #domain.cfl = 1.0 #domain.limiter = "minmod" domain.cfl = 1.0 domain.limiter = "vanleer" def height(x): y=zeros(len(x), Float) return y domain.set_quantity('height', height) domain.set_quantity('stage',height) #('height', height) domain.order=domain.default_order print "domain order", domain.order t0=time.time() yieldstep=30.0 finaltime=30.0 print "integral", domain.quantities['height'].get_integral() finaltime=20.0 print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral() for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime): domain.write_time() print "integral", domain.quantities['height'].get_integral() print "integral", domain.quantities['stage'].get_integral() #['height'].get_integral() if t>0.0: HeightQ=domain.quantities['height'].vertex_values HeightQ=domain.quantities['stage'].vertex_values #['height'].vertex_values MomentumQ=domain.quantities['xmomentum'].vertex_values h, uh=analytical_sol(X.flat, domain.time)
• ## anuga_work/development/anuga_1d/dry_dam_sudi.py

 r5565 # Calculate Analytical Solution at time t > 0 u3 = 2.0/3.0*(sqrt(g*h1)+x[i]/t) h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t)) h3 = 4.0/(9.0*g)*(sqrt(g*h1)-x[i]/(2.0*t))*(sqrt(g*h1)-x[i]/(2.0*t)) u3_ = 2.0/3.0*((x[i]+L/2.0)/t-sqrt(g*h1)) h3_ = 1.0/(9.0*g)*((x[i]+L/2.0)/t+2*sqrt(g*h1))*((x[i]+L/2.0)/t+2*sqrt(g*h1)) if ( x[i] <= -t*sqrt(g*h1) ): if ( x[i] <= -1*L/2.0+2*(-sqrt(g*h1)*t)): u[i] = 0.0 h[i] = h0 elif ( x[i] <= -1*L/2.0-(-sqrt(g*h1)*t)): u[i] = u3_ h[i] = h3_ elif ( x[i] <= -t*sqrt(g*h1) ): u[i] = 0.0 h[i] = h1 L = 2000.0     # Length of channel (m) N = 100        # Number of compuational cells N = 800        # Number of computational cells cell_len = L/N # Origin = 0.0 yieldstep = 1.0 L = 2000.0     # Length of channel (m) number_of_cells = [100]#,200,500,1000,2000,5000,10000,20000] number_of_cells = [200]#,200,500,1000,2000,5000,10000,20000] h_error = zeros(len(number_of_cells),Float) uh_error = zeros(len(number_of_cells),Float) domain.set_boundary({'exterior': Reflective_boundary(domain)}) domain.default_order = 2 domain.default_time_order = 1 domain.default_time_order = 2 domain.cfl = 1.0 domain.limiter = "vanleer" domain.limiter = "minmod" t0 = time.time() for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): pass domain.write_time() N = float(N)
• ## anuga_work/development/anuga_1d/run_profile_dry.py

 r5565 S = pstats.Stats(FN) s = S.sort_stats('cumulative').print_stats(20) s = S.sort_stats('cumulative').print_stats(30) print s
• ## anuga_work/development/anuga_1d/shallow_water_domain.py

 r5563 tagged_elements) from config import minimum_allowed_height, g from config import minimum_allowed_height, g, h0 self.minimum_allowed_height = minimum_allowed_height self.g = g self.h0 = h0 #forcing terms not included in 1d domain ?WHy? """ from config import g, epsilon from config import g, epsilon, h0 from math import sqrt from Numeric import array h_left = 0.0 else: u_left  = uh_left/h_left u_left  = uh_left/(h_left +  h0/h_left) uh_left = u_left*h_left h_right  = w_right-z uh_right = q_right[1] if h_right < epsilon: h_right = 0.0 else: u_right  = uh_right/h_right u_right  = uh_right/(h_right + h0/h_right) uh_right = u_right*h_right #vh_left  = q_left[2]
• ## anuga_work/development/anuga_1d/test_shallow_water.py

 r5563 zr = 0.0 #This assumes h0 = 1.0e-3!! edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr) assert allclose(array([2.0, 8.9],Float), edgeflux) assert allclose(5.1304951685, maxspeed) #print maxspeed #print edgeflux assert allclose(array([1.998002, 8.89201198],Float), edgeflux) assert allclose(5.1284971665, maxspeed) normal = -1.0 edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr) assert allclose(array([-2.0, -8.9],Float), edgeflux) assert allclose(5.1304951685, maxspeed) #print maxspeed #print edgeflux assert allclose(array([-1.998002, -8.89201198],Float), edgeflux) assert allclose(5.1284971665, maxspeed) def test_domain_flux_function(self): """ from config import g, epsilon from config import g, epsilon, h0 from math import sqrt from Numeric import array h_left = 0.0 else: u_left  = uh_left/h_left u_left  = uh_left/(h_left +  h0/h_left) uh_left = u_left*h_left w_right  = q_right[0]  #w=h+z h_right = 0.0 else: u_right  = uh_right/h_right u_right  = uh_right/(h_right + h0/h_right) uh_right = u_right*h_right #vh_left  = q_left[2] #vh_right = q_right[2]
Note: See TracChangeset for help on using the changeset viewer.