Changeset 5563


Ignore:
Timestamp:
Jul 23, 2008, 4:38:10 PM (11 years ago)
Author:
steve
Message:

Working version of comp_flux_ext

Location:
anuga_work/development/anuga_1d
Files:
1 added
4 edited

Legend:

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

    r5562 r5563  
    55const double pi = 3.14159265358979;
    66
    7 double max(double a, double b) {
    8         double z;
    9         z=(a>b)?a:b;
    10         return z;}
    11        
    12 double min(double a, double b) {
    13         double z;
    14         z=(a<b)?a:b;
    15         return z;}
     7// Shared code snippets
     8#include "util_ext.h"
     9
     10
     11/* double max(double a, double b) { */
     12/*      double z; */
     13/*      z=(a>b)?a:b; */
     14/*      return z;} */
     15       
     16/* double min(double a, double b) { */
     17/*      double z; */
     18/*      z=(a<b)?a:b; */
     19/*      return z;} */
    1620       
    1721
     
    154158                               
    155159                                normal = normals[ki];
    156                                 _flux_function(ql, qr, zl, zr, normal, g, epsilon, flux, &max_speed);
     160                                _flux_function(ql, qr, zl, zr, normal, g, epsilon, edgeflux, &max_speed);
    157161                                flux[0] -= edgeflux[0];
    158162                                flux[1] -= edgeflux[1];
     
    189193//=========================================================================
    190194PyObject *compute_fluxes_ext(PyObject *self, PyObject *args) {
    191   PyArrayObject *neighbours,
     195 
     196    PyArrayObject
     197        *neighbours,
    192198        *neighbour_vertices,
    193     *normals,
     199        *normals,
    194200        *areas,
    195     *stage_edge_values,
    196     *xmom_edge_values,
    197     *bed_edge_values,
    198     *stage_boundary_values,
    199     *xmom_boundary_values,
     201        *stage_edge_values,
     202        *xmom_edge_values,
     203        *bed_edge_values,
     204        *stage_boundary_values,
     205        *xmom_boundary_values,
    200206        *stage_explicit_update,
    201207        *xmom_explicit_update,
     
    231237  // the explicit update arrays
    232238  timestep = _compute_fluxes_ext(timestep,
    233                                      epsilon,
    234                                      g,
    235                                      (long*) neighbours -> data,
    236                                      (long*) neighbour_vertices -> data,
    237                                      (double*) normals -> data,
    238                                      (double*) areas -> data,
    239                                      (double*) stage_edge_values -> data,
    240                                      (double*) xmom_edge_values -> data,
    241                                      (double*) bed_edge_values -> data,
    242                                      (double*) stage_boundary_values -> data,
    243                                      (double*) xmom_boundary_values -> data,
    244                                          (double*) stage_explicit_update -> data,
    245                                          (double*) xmom_explicit_update -> data,
    246                                          number_of_elements,
    247                                          (double*) max_speed_array -> data);
     239                                 epsilon,
     240                                 g,
     241                                 (long*) neighbours -> data,
     242                                 (long*) neighbour_vertices -> data,
     243                                 (double*) normals -> data,
     244                                 (double*) areas -> data,
     245                                 (double*) stage_edge_values -> data,
     246                                 (double*) xmom_edge_values -> data,
     247                                 (double*) bed_edge_values -> data,
     248                                 (double*) stage_boundary_values -> data,
     249                                 (double*) xmom_boundary_values -> data,
     250                                 (double*) stage_explicit_update -> data,
     251                                 (double*) xmom_explicit_update -> data,
     252                                 number_of_elements,
     253                                 (double*) max_speed_array -> data);
     254
    248255  // Return updated flux timestep
    249256  return Py_BuildValue("d", timestep);
    250257}
    251258
     259
     260PyObject *compute_fluxes_ext_short(PyObject *self, PyObject *args) {
     261 
     262   PyObject
     263        *domain,
     264        *stage,
     265        *xmom,
     266        *bed;
     267
     268    PyArrayObject
     269        *neighbours,
     270        *neighbour_vertices,
     271        *normals,
     272        *areas,
     273        *stage_vertex_values,
     274        *xmom_vertex_values,
     275        *bed_vertex_values,
     276        *stage_boundary_values,
     277        *xmom_boundary_values,
     278        *stage_explicit_update,
     279        *xmom_explicit_update,
     280        *max_speed_array;
     281   
     282  double timestep, epsilon, g;
     283  int number_of_elements;
     284
     285   
     286  // Convert Python arguments to C
     287  if (!PyArg_ParseTuple(args, "dOOOO",
     288                        &timestep,
     289                        &domain,
     290                        &stage,
     291                        &xmom,
     292                        &bed)) {
     293      PyErr_SetString(PyExc_RuntimeError, "comp_flux_ext.c: compute_fluxes_ext_short could not parse input");
     294      return NULL;
     295  }
     296
     297
     298    epsilon           = get_python_double(domain,"epsilon");
     299    g                 = get_python_double(domain,"g");
     300 
     301   
     302    neighbours        = get_consecutive_array(domain, "neighbours");
     303    neighbour_vertices= get_consecutive_array(domain, "neighbour_vertices");
     304    normals           = get_consecutive_array(domain, "normals");
     305    areas             = get_consecutive_array(domain, "areas");   
     306    max_speed_array   = get_consecutive_array(domain, "max_speed_array");
     307   
     308    stage_vertex_values = get_consecutive_array(stage, "vertex_values");   
     309    xmom_vertex_values  = get_consecutive_array(xmom, "vertex_values");   
     310    bed_vertex_values   = get_consecutive_array(bed, "vertex_values");   
     311
     312    stage_boundary_values = get_consecutive_array(stage, "boundary_values");   
     313    xmom_boundary_values  = get_consecutive_array(xmom, "boundary_values");   
     314
     315
     316    stage_explicit_update = get_consecutive_array(stage, "explicit_update");   
     317    xmom_explicit_update  = get_consecutive_array(xmom, "explicit_update");   
     318
     319
     320
     321    number_of_elements = stage_vertex_values -> dimensions[0];
     322
     323
     324 
     325    // Call underlying flux computation routine and update
     326    // the explicit update arrays
     327    timestep = _compute_fluxes_ext(timestep,
     328                                 epsilon,
     329                                 g,
     330                                 (long*) neighbours -> data,
     331                                 (long*) neighbour_vertices -> data,
     332                                 (double*) normals -> data,
     333                                 (double*) areas -> data,
     334                                 (double*) stage_vertex_values -> data,
     335                                 (double*) xmom_vertex_values -> data,
     336                                 (double*) bed_vertex_values -> data,
     337                                 (double*) stage_boundary_values -> data,
     338                                 (double*) xmom_boundary_values -> data,
     339                                 (double*) stage_explicit_update -> data,
     340                                 (double*) xmom_explicit_update -> data,
     341                                 number_of_elements,
     342                                 (double*) max_speed_array -> data);
     343
     344
     345  Py_DECREF(neighbours);
     346  Py_DECREF(neighbour_vertices);
     347  Py_DECREF(normals);
     348  Py_DECREF(areas);
     349  Py_DECREF(stage_vertex_values);
     350  Py_DECREF(xmom_vertex_values);
     351  Py_DECREF(bed_vertex_values);
     352  Py_DECREF(stage_boundary_values);
     353  Py_DECREF(xmom_boundary_values);
     354  Py_DECREF(stage_explicit_update);
     355  Py_DECREF(xmom_explicit_update);
     356  Py_DECREF(max_speed_array);
     357
     358
     359
     360
     361  // Return updated flux timestep
     362  return Py_BuildValue("d", timestep);
     363}
     364
    252365 
    253366
     
    259372static struct PyMethodDef MethodTable[] = {
    260373  {"compute_fluxes_ext", compute_fluxes_ext, METH_VARARGS, "Print out"},
     374  {"compute_fluxes_ext_short", compute_fluxes_ext_short, METH_VARARGS, "Print out"},
    261375  {NULL, NULL}
    262376};
  • anuga_work/development/anuga_1d/domain.py

    r5536 r5563  
    4747        self.centroids = zeros(N, Float)
    4848        self.areas     = zeros(N, Float)
     49
     50        self.max_speed_array = zeros(N, Float)
    4951
    5052        self.normals = zeros((N, 2), Float)
  • anuga_work/development/anuga_1d/shallow_water_domain.py

    r5536 r5563  
    9090
    9191        self.quantities_to_be_stored = ['stage','xmomentum']
     92
     93
    9294
    9395
     
    155157        #Call correct module function
    156158        #(either from this module or C-extension)
    157         compute_fluxes(self)
     159        compute_fluxes_C_short(self)
    158160
    159161    def compute_timestep(self):
     
    572574    #print domain.quantities['stage'].centroid_values
    573575
     576
     577# Compute flux definition
     578def compute_fluxes_C_long(domain):
     579        from Numeric import zeros, Float
     580        import sys
     581
     582       
     583        timestep = float(sys.maxint)
     584        #print 'timestep=',timestep
     585        #print 'The type of timestep is',type(timestep)
     586       
     587        epsilon = domain.epsilon
     588        #print 'epsilon=',epsilon
     589        #print 'The type of epsilon is',type(epsilon)
     590       
     591        g = domain.g
     592        #print 'g=',g
     593        #print 'The type of g is',type(g)
     594       
     595        neighbours = domain.neighbours
     596        #print 'neighbours=',neighbours
     597        #print 'The type of neighbours is',type(neighbours)
     598       
     599        neighbour_vertices = domain.neighbour_vertices
     600        #print 'neighbour_vertices=',neighbour_vertices
     601        #print 'The type of neighbour_vertices is',type(neighbour_vertices)
     602       
     603        normals = domain.normals
     604        #print 'normals=',normals
     605        #print 'The type of normals is',type(normals)
     606       
     607        areas = domain.areas
     608        #print 'areas=',areas
     609        #print 'The type of areas is',type(areas)
     610       
     611        stage_edge_values = domain.quantities['stage'].vertex_values
     612        #print 'stage_edge_values=',stage_edge_values
     613        #print 'The type of stage_edge_values is',type(stage_edge_values)
     614       
     615        xmom_edge_values = domain.quantities['xmomentum'].vertex_values
     616        #print 'xmom_edge_values=',xmom_edge_values
     617        #print 'The type of xmom_edge_values is',type(xmom_edge_values)
     618       
     619        bed_edge_values = domain.quantities['elevation'].vertex_values
     620        #print 'bed_edge_values=',bed_edge_values
     621        #print 'The type of bed_edge_values is',type(bed_edge_values)
     622       
     623        stage_boundary_values = domain.quantities['stage'].boundary_values
     624        #print 'stage_boundary_values=',stage_boundary_values
     625        #print 'The type of stage_boundary_values is',type(stage_boundary_values)
     626       
     627        xmom_boundary_values = domain.quantities['xmomentum'].boundary_values
     628        #print 'xmom_boundary_values=',xmom_boundary_values
     629        #print 'The type of xmom_boundary_values is',type(xmom_boundary_values)
     630       
     631        stage_explicit_update = domain.quantities['stage'].explicit_update
     632        #print 'stage_explicit_update=',stage_explicit_update
     633        #print 'The type of stage_explicit_update is',type(stage_explicit_update)
     634       
     635        xmom_explicit_update = domain.quantities['xmomentum'].explicit_update
     636        #print 'xmom_explicit_update=',xmom_explicit_update
     637        #print 'The type of xmom_explicit_update is',type(xmom_explicit_update)
     638       
     639        number_of_elements = len(stage_edge_values)
     640        #print 'number_of_elements=',number_of_elements
     641        #print 'The type of number_of_elements is',type(number_of_elements)
     642       
     643        max_speed_array = domain.max_speed_array
     644        #print 'max_speed_array=',max_speed_array
     645        #print 'The type of max_speed_array is',type(max_speed_array)
     646       
     647       
     648        from comp_flux_ext import compute_fluxes_ext
     649       
     650        domain.timestep = compute_fluxes_ext(timestep,
     651                                  epsilon,
     652                                  g,
     653                                  neighbours,
     654                                  neighbour_vertices,
     655                                  normals,
     656                                  areas,
     657                                  stage_edge_values,
     658                                  xmom_edge_values,
     659                                  bed_edge_values,
     660                                  stage_boundary_values,
     661                                  xmom_boundary_values,
     662                                  stage_explicit_update,
     663                                  xmom_explicit_update,
     664                                  number_of_elements,
     665                                  max_speed_array)
     666
     667
     668# Compute flux definition
     669def compute_fluxes_C_short(domain):
     670        from Numeric import zeros, Float
     671        import sys
     672
     673       
     674        timestep = float(sys.maxint)
     675
     676        stage = domain.quantities['stage']
     677        xmom = domain.quantities['xmomentum']
     678        bed = domain.quantities['elevation']
     679
     680
     681        from comp_flux_ext import compute_fluxes_ext_short
     682
     683        domain.timestep = compute_fluxes_ext_short(timestep,domain,stage,xmom,bed)
     684
     685
     686
    574687####################################
     688
     689
     690
     691
     692
    575693
    576694# Module functions for gradient limiting (distribute_to_vertices_and_edges)
  • anuga_work/development/anuga_1d/test_shallow_water.py

    r5538 r5563  
    4141        assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
    4242        assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
     43
     44
     45    def test_local_flux_function(self):
     46        normal = 1.0
     47        ql = array([1.0, 2.0],Float)
     48        qr = array([1.0, 2.0],Float)
     49        zl = 0.0
     50        zr = 0.0
     51
     52        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
     53
     54        assert allclose(array([2.0, 8.9],Float), edgeflux)
     55        assert allclose(5.1304951685, maxspeed)
     56
     57        normal = -1.0
     58        ql = array([1.0, 2.0],Float)
     59        qr = array([1.0, 2.0],Float)
     60        zl = 0.0
     61        zr = 0.0
     62
     63        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
     64
     65        assert allclose(array([-2.0, -8.9],Float), edgeflux)
     66        assert allclose(5.1304951685, maxspeed)
     67
     68    def test_domain_flux_function(self):
     69        normal = 1.0
     70        ql = array([1.0, 2.0],Float)
     71        qr = array([1.0, 2.0],Float)
     72        zl = 0.0
     73        zr = 0.0
     74
     75        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
     76
     77        #print edgeflux
     78
     79        from shallow_water_domain import flux_function as domain_flux_function
     80
     81        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
     82
     83        #print domainedgeflux
     84
     85        assert allclose(domainedgeflux, edgeflux)
     86        assert allclose(domainmaxspeed, maxspeed)
     87
    4388
    4489
     
    128173                #m = domain.neighbour_edges[k,i]
    129174                m = domain.neighbour_vertices[k,i]
     175                #print i, ' ' , m
    130176                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
    131177                qr[0] = stage[n, m]
     
    169215
    170216
     217def flux_function(normal, ql, qr, zl, zr):
     218    """Compute fluxes between volumes for the shallow water wave equation
     219    cast in terms of w = h+z using the 'central scheme' as described in
     220
     221    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
     222    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
     223    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
     224
     225    The implemented formula is given in equation (3.15) on page 714
     226
     227    Conserved quantities w, uh, are stored as elements 0 and 1
     228    in the numerical vectors ql an qr.
     229
     230    Bed elevations zl and zr.
     231    """
     232
     233    from config import g, epsilon
     234    from math import sqrt
     235    from Numeric import array
     236
     237    #print 'ql',ql
     238
     239    #Align momentums with x-axis
     240    #q_left  = rotate(ql, normal, direction = 1)
     241    #q_right = rotate(qr, normal, direction = 1)
     242    q_left = ql
     243    q_left[1] = q_left[1]*normal
     244    q_right = qr
     245    q_right[1] = q_right[1]*normal
     246
     247    #z = (zl+zr)/2 #Take average of field values
     248    z = 0.5*(zl+zr) #Take average of field values
     249
     250    w_left  = q_left[0]   #w=h+z
     251    h_left  = w_left-z
     252    uh_left = q_left[1]
     253
     254    if h_left < epsilon:
     255        u_left = 0.0  #Could have been negative
     256        h_left = 0.0
     257    else:
     258        u_left  = uh_left/h_left
     259
     260
     261    w_right  = q_right[0]  #w=h+z
     262    h_right  = w_right-z
     263    uh_right = q_right[1]
     264
     265
     266    if h_right < epsilon:
     267        u_right = 0.0  #Could have been negative
     268        h_right = 0.0
     269    else:
     270        u_right  = uh_right/h_right
     271
     272    #vh_left  = q_left[2]
     273    #vh_right = q_right[2]
     274
     275    #print h_right
     276    #print u_right
     277    #print h_left
     278    #print u_right
     279
     280    soundspeed_left  = sqrt(g*h_left)
     281    soundspeed_right = sqrt(g*h_right)
     282
     283    #Maximal wave speed
     284    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
     285   
     286    #Minimal wave speed
     287    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
     288   
     289    #Flux computation
     290
     291    #flux_left  = array([u_left*h_left,
     292    #                    u_left*uh_left + 0.5*g*h_left**2])
     293    #flux_right = array([u_right*h_right,
     294    #                    u_right*uh_right + 0.5*g*h_right**2])
     295    flux_left  = array([u_left*h_left,
     296                        u_left*uh_left + 0.5*g*h_left*h_left])
     297    flux_right = array([u_right*h_right,
     298                        u_right*uh_right + 0.5*g*h_right*h_right])
     299
     300    denom = s_max-s_min
     301    if denom == 0.0:
     302        edgeflux = array([0.0, 0.0])
     303        max_speed = 0.0
     304    else:
     305        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
     306        edgeflux += s_max*s_min*(q_right-q_left)/denom
     307       
     308        edgeflux[1] = edgeflux[1]*normal
     309
     310        max_speed = max(abs(s_max), abs(s_min))
     311
     312    return edgeflux, max_speed
     313
    171314
    172315#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.