# Changeset 5162

Ignore:
Timestamp:
Mar 11, 2008, 8:43:22 PM (16 years ago)
Message:

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

Files:
14 edited

Unmodified
Removed
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/domain.py

 r5080 from Numeric import allclose, argmax, zeros, Float from anuga.config import epsilon from anuga.config import beta_euler, beta_rk2 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh self.default_order = n self._order_ = self.default_order if self.default_order == 1: self.set_timestepping_method('euler') #self.set_all_limiters(beta_euler) if self.default_order == 2: self.set_timestepping_method('rk2') #self.set_all_limiters(beta_rk2) elif self._order_ == 2: Q.extrapolate_second_order() Q.limit() #Q.limit() else: raise 'Unknown order' Q.interpolate_from_vertices_to_edges() #Q.interpolate_from_vertices_to_edges()
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

 r4990 #Allocate space for Gradient self.x_gradient = zeros(N, Float) self.y_gradient = zeros(N, Float) self.y_gradient = zeros(N, Float) #Allocate space for Limiter Phi self.phi = zeros(N, Float) #Intialise centroid and edge_values #flux calculations and forcing functions N = len(domain) # number_of_triangles #Allocate space for update fields self.explicit_update = zeros(N, Float ) self.semi_implicit_update = zeros(N, Float ) self.centroid_backup_values = zeros(N, Float) self.centroid_backup_values = zeros(N, Float) self.beta = 1.0 compute_gradients(self) extrapolate_from_gradient(self) def extrapolate_second_order_and_limit(self): #Call correct module function #(either from this module or C-extension) extrapolate_second_order_and_limit(self) extrapolate_second_order_and_limit(self) def bound_vertices_below_by_constant(self, bound): #Call correct module function #(either from this module or C-extension) bound_vertices_below_by_constant(self, bound) def bound_vertices_below_by_quantity(self, quantity): #Call correct module function #(either from this module or C-extension) #check consistency assert self.domain == quantity.domain bound_vertices_below_by_quantity(self, quantity) def backup_centroid_values(self): limit_gradient_by_neighbour,\ extrapolate_from_gradient,\ extrapolate_second_order_and_limit,\ bound_vertices_below_by_constant,\ bound_vertices_below_by_quantity,\ interpolate_from_vertices_to_edges,\ interpolate_from_edges_to_vertices,\
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

 r4978 // // To compile (Python2.X): // use python ../utilities/compile.py // // See the module advection.py // use python ../utilities/compile.py to // compile all C files within a directory // // //------------------------------------------- // Low level routines (called from wrappers) //------------------------------------------ double _compute_fluxes( //Loop //Local Variables double qr,ql; int k3; //Loop through triangles timestep = max_timestep; for (i=0; i<3; i++){ k_i = k3+i; //Quantities inside volume facing neighbour i //Quantities inside triangle facing neighbour i ql = quantity_edge[k_i]; PyObject *domain, *quantity; PyArrayObject * quantity_update,
• ## anuga_core/source/anuga/config.py

 r4837 #timestepping_method = 'rk2'   # 2nd Order TVD scheme # rk2 is a little more stable than euler, so rk2 timestepping # can deal with a larger beta when slope limiting the reconstructed # solution. The large beta is needed if solving problems sensitive # to numerical diffusion, like a small forced wave in an ocean beta_euler = 1.0 beta_rk2   = 1.6 # Option to search for signatures where isolated triangles are # responsible for a small global timestep.
• ## anuga_core/source/anuga/shallow_water/shallow_water_domain.py

 r5089 from anuga.config import alpha_balance from anuga.config import optimise_dry_cells from anuga.config import optimised_gradient_limiter #--------------------- self.minimum_storable_height = minimum_storable_height self.quantities_to_be_stored = ['stage','xmomentum','ymomentum'] # Limiters self.use_old_limiter = True self.optimised_gradient_limiter = optimised_gradient_limiter self.beta_w      = beta self.beta_w_dry  = beta self.quantities['stage'].beta = beta self.beta_uh     = beta self.beta_uh_dry = beta self.quantities['xmomentum'].beta = beta self.beta_vh     = beta self.beta_vh_dry = beta self.quantities['ymomentum'].beta = beta self.beta_h      = beta # Call correct module function # (either from this module or C-extension) distribute_to_vertices_and_edges(self) if self.use_old_limiter: distribute_to_vertices_and_edges(self) else: for name in self.conserved_quantities: Q = self.quantities[name] if self._order_ == 1: Q.extrapolate_first_order() elif self._order_ == 2: Q.extrapolate_second_order_and_limit() if name == 'stage': Q.bound_vertices_below_by_quantity(self.quantities['elevation']) else: raise 'Unknown order' # self.check_integrity() msg = 'Parameter beta_h must be in the interval [0, 1[' assert 0 <= self.beta_h <= 1.0, msg msg = 'Parameter beta_w must be in the interval [0, 1[' assert 0 <= self.beta_w <= 1.0, msg msg = 'Parameter beta_h must be in the interval [0, 2[' assert 0 <= self.beta_h <= 2.0, msg msg = 'Parameter beta_w must be in the interval [0, 2[' assert 0 <= self.beta_w <= 2.0, msg """ from anuga.config import optimised_gradient_limiter # Remove very thin layers of water # Extrapolate all conserved quantities if optimised_gradient_limiter: if domain.optimised_gradient_limiter: # MH090605 if second order, # perform the extrapolation and limiting on Q.extrapolate_first_order() elif domain._order_ == 2: # Experiment #if name == 'stage': #    #print name, 'second' #    Q.extrapolate_second_order() #    Q.limit() #else: #    #print name, 'first' #    Q.extrapolate_first_order() #    #Q.extrapolate_second_order() #    #Q.limit() Q.extrapolate_second_order() Q.limit()
• ## anuga_core/source/anuga/shallow_water/shallow_water_ext.c

 r4978 double dx1,dx2,dy1,dy2,dxv0,dxv1,dxv2,dyv0,dyv1,dyv2,dq0,dq1,dq2,area2; double dqv[3], qmin, qmax, hmin, hmax; double hc, h0, h1, h2, beta_tmp; double hc, h0, h1, h2, beta_tmp, hfactor; h1 = stage_centroid_values[k1] - elevation_centroid_values[k1]; h2 = stage_centroid_values[k2] - elevation_centroid_values[k2]; hmin = min(h0,min(h1,h2)); hmin = min(min(h0,min(h1,h2)),hc); //hfactor = hc/(hc + 1.0); hfactor = 0.0; if (hmin > 0.1 ) { hfactor = (hmin-0.1)/(hmin+0.4); } if (optimise_dry_cells) { // Playing with dry wet interface hmin = qmin; beta_tmp = beta_w; if (hminminimum_allowed_height) beta_tmp = beta_w_dry + (beta_w - beta_w_dry) * hfactor; //printf("min_alled_height = %f\n",minimum_allowed_height); //printf("hmin = %f\n",hmin); //printf("beta_w = %f\n",beta_w); //printf("beta_tmp = %f\n",beta_tmp); // Limit the gradient limit_gradient(dqv,qmin,qmax,beta_tmp); // from the centroid to the min and max find_qmin_and_qmax(dq0,dq1,dq2,&qmin,&qmax); beta_tmp = beta_uh; if (hmin dimensions[0];
• ## anuga_core/source/anuga/shallow_water/test_data_manager.py

 r4924 range = fid.variables['xmomentum_range'][:] assert allclose(range,[0,0.4695096]) or\ allclose(range,[0,0.47790655]) # Old slope limiters assert allclose(range,[0,0.4695096]) or \ allclose(range,[0,0.47790655]) or\ allclose(range,[0,0.46941957]) range = fid.variables['ymomentum_range'][:] #assert allclose(range,[0,0.02174380]) assert allclose(range,[0,0.02174439]) or\ allclose(range,[0,0.02283983]) # Old slope limiters allclose(range,[0,0.02283983]) or\ allclose(range,[0,0.0217342]) # Old slope limiters fid.close() extrema = fid.variables['xmomentum.extrema'][:] assert allclose(extrema,[-0.06062178, 0.47886313]) assert allclose(extrema,[-0.06062178, 0.47873023]) extrema = fid.variables['ymomentum.extrema'][:] assert allclose(extrema,[0.00, 0.06241221]) assert allclose(extrema,[0.00, 0.0625786]) time_interval = fid.variables['extrema.time_interval'][:]
• ## anuga_core/source/anuga/utilities/util_ext.h

 r4978 B = (PyArrayObject*) PyObject_GetAttrString(O, name); if (!B) return NULL; if (!B) { PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array"); return NULL; } //Convert to consecutive array Py_DECREF(B); //FIXME: Is this really needed?? if (!A) return NULL; if (!A) { PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_consecutive_array could not obtain array"); return NULL; } return A; } double get_double(PyObject *O, char *name) { double get_python_double(PyObject *O, char *name) { PyObject *TObject; double tmp; TObject = PyObject_GetAttrString(O, name); if (!TObject) { PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_double could not obtain double from object"); PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_double could not obtain double from object"); return 0.0; } PyObject *get_python_object(PyObject *O, char *name) { PyObject *Oout; Oout = PyObject_GetAttrString(O, name); if (!Oout) { PyErr_SetString(PyExc_RuntimeError, "util_ext.h: get_python_object could not obtain object"); return NULL; } return Oout; }
• ## anuga_core/source/anuga/utilities/xml_tools.py

 r5022 fid = xml try: dom = parse(fid) except Exception, e: # Throw filename into dom exception msg = 'XML file "%s" could not be parsed: ' %fid.name msg += str(e) raise Exception, msg dom = parse(fid) ##     try: ##         dom = parse(fid) ##     except Exception, e: ##         # Throw filename into dom exception ##         msg = 'XML file "%s" could not be parsed: ' %fid.name ##         msg += str(e) ##         raise Exception, msg return dom2object(dom)
• ## anuga_validation/convergence_study/animatesww2d_alt.py

 r4838 def animatesww2d_alt(sww_filename, movie_filename, range): """Plot cross section of model output """ ## def animatesww2d_alt(sww_filename, movie_filename, range): ##     """Plot cross section of model output ##     """ # Read SWW file from Scientific.IO.NetCDF import NetCDFFile fid = NetCDFFile(sww_filename, 'r') ##     # Read SWW file ##     from Scientific.IO.NetCDF import NetCDFFile ##     fid = NetCDFFile(sww_filename, 'r') x = fid.variables['x'] y = fid.variables['y'] volumes = fid.variables['volumes'] ##     x = fid.variables['x'] ##     y = fid.variables['y'] ##     volumes = fid.variables['volumes'] elevation = fid.variables['elevation'] time = fid.variables['time'] stage = fid.variables['stage'] xmomentum = fid.variables['xmomentum'] ymomentum = fid.variables['ymomentum'] ##     elevation = fid.variables['elevation'] ##     time = fid.variables['time'] ##     stage = fid.variables['stage'] ##     xmomentum = fid.variables['xmomentum'] ##     ymomentum = fid.variables['ymomentum'] def animatesww2d(swwfile): # interpolation points are for y = 0 # number of increments in x m = 100 m = 1000 max_x = 100000. x = 0
• ## anuga_validation/convergence_study/convergence_structured.py

 r4992 from anuga.shallow_water import Transmissive_boundary from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary from anuga.geospatial_data.geospatial_data import * #from anuga.geospatial_data.geospatial_data import * from math import cos from Numeric import zeros, Float #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ dx = 1000. dx = 200. dy = dx L = 100000. domain = Domain(points, vertices, boundary) domain.set_timestepping_method('euler') domain.set_timestepping_method('rk2') domain.set_default_order(2) domain.set_name('myexample9') domain.set_datadir('.')                     # Use current directory for output domain.beta_w      = 1.0 domain.beta_w_dry  = 0.2 domain.beta_uh     = 1.0 domain.beta_uh_dry = 0.2 domain.beta_vh     = 1.0 domain.beta_vh_dry = 0.2 domain.beta_h      = 1.0 domain.set_all_limiters(0.9) print domain.beta_w domain.use_old_limiter = False domain.CFL = 1.0 #------------------------------------------------------------------------------ #------------------------------------------------------------------------------ #domain.set_quantity('elevation', topography) # Use function for elevation domain.set_quantity('elevation',-100) domain.set_quantity('elevation',0.0) domain.set_quantity('friction', 0.00) domain.set_quantity('stage', 0.0) h0 = 10.0 h1 = 1.0 def height(x,y): z = zeros(len(x),Float) for i in range(len(x)): if x[i]<=50000.0: z[i] = h0 else: z[i] = h1 return z domain.set_quantity('stage', height) #domain.set_quantity('stage', 0.0) #----------------------------------------------------------------------------- Bw = Time_boundary(domain=domain,     # Time dependent boundary ## Sine wave f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0]) #                   f=lambda t: [(-amplitude*sin((1./300.)*t*2*pi)), 0.0, 0.0]) ## Single wave f=lambda t: [h0, 0.0, 0.0]) ## Sawtooth? #                   f=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0]) #------------------------------------------------------------------------------ for t in domain.evolve(yieldstep = 20.0, finaltime = 10*40*60.): for t in domain.evolve(yieldstep = 50.0, finaltime = 10*40*60.): domain.write_time() vis.update()
• ## anuga_work/development/analytical_solutions/oblique_shock.py

 r5030 from shallow_water import Domain, Constant_height from shallow_water import Transmissive_boundary, Reflective_boundary,\ from anuga.shallow_water import Domain from anuga.shallow_water import Transmissive_boundary, Reflective_boundary,\ Dirichlet_boundary from anuga.visualiser import RealtimeVisualiser from math import sqrt, cos, sin, pi from mesh_factory import oblique from mesh_factory import oblique_cross leny = 30. lenx = 40. n = 50 m = 60 n = 100 m = 120 points, elements, boundary = oblique(m, n, lenx, leny) points, elements, boundary = oblique_cross(m, n, lenx, leny) domain = Domain(points, elements, boundary) # Store output domain.store=True #domain.store=True # Output format domain.format="sww" #NET.CDF binary format #domain.format="sww" #NET.CDF binary format # "dat" for ASCII #Initial condition h = 0.5 domain.set_quantity('stage', Constant_height(x_slope, h) ) domain.set_quantity('stage', expression = 'elevation + %d'%h ) #--------------------------------- # Setup visualization #--------------------------------- vis = RealtimeVisualiser(domain) vis.render_quantity_height("elevation", dynamic=False) vis.render_quantity_height("stage", zScale=10, dynamic=True) vis.colour_height_quantity('stage', (0.75, 0.5, 0.5)) vis.start() ###################### for t in domain.evolve(yieldstep = 0.5, finaltime = 50): domain.write_time() vis.update() print 'That took %.2f seconds' %(time.time()-t0) vis.evolveFinished() vis.join() #FIXME: Compute average water depth on either side of shock and compare
• ## anuga_work/development/near_shore_PMD/run_beach.py

 r5154 domain.set_boundary( {'wall': Br, 'wave': Bwp_velocity} ) #------------------------------------------------------------------------- # Evolve system through time for t in domain.evolve(yieldstep, finaltime): domain.write_time() print 'That took %.2f seconds' %(time.time()-t0) print 'finished'
Note: See TracChangeset for help on using the changeset viewer.