Changeset 4883


Ignore:
Timestamp:
Dec 11, 2007, 6:48:46 PM (16 years ago)
Author:
steve
Message:

Starting to develop edge limiting and vertex limiting in
quantity (and later in shallow_water_domain

Location:
anuga_core
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • anuga_core/documentation/user_manual/demos/runup.py

    r4551 r4883  
    1010#------------------------------------------------------------------------------
    1111
    12 from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
     12from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
    1313from anuga.shallow_water import Domain
    1414from anuga.shallow_water import Reflective_boundary
    … …  
    2222#------------------------------------------------------------------------------
    2323
    24 points, vertices, boundary = rectangular(10, 10) # Basic mesh
     24points, vertices, boundary = rectangular_cross(10, 10,len1=1.0, len2= 1.0 ) # Basic mesh
    2525
    2626domain = Domain(points, vertices, boundary) # Create domain
    … …  
    3737
    3838def topography(x,y):
    39     return -x/2                             # linear bed slope
     39    return -x/2                              # linear bed slope
    4040    #return x*(-(2.0-x)*.5)                  # curved bed slope
    4141
    … …  
    5454Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
    5555Bw = Time_boundary(domain=domain,     # Time dependent boundary 
    56                    f=lambda t: [(.1*sin(t*2*pi)-0.3) * exp(-t), 0.0, 0.0])
     56                   f=lambda t: [(.1*sin(t*2*pi)-0.3) , 0.0, 0.0])
    5757
    5858# Associate boundary tags with boundary objects
    5959domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
    6060
     61#===============================================================================
     62from anuga.visualiser import RealtimeVisualiser
     63vis = RealtimeVisualiser(domain)
     64vis.render_quantity_height("elevation", zScale=1.0, dynamic=False)
     65vis.render_quantity_height("stage", dynamic=True)
     66vis.colour_height_quantity('stage', (lambda q:q['stage'], -0.5, 0.5))
     67vis.start()
     68#===============================================================================
    6169
     70import time
     71t0 = time.time()
    6272#------------------------------------------------------------------------------
    6373# Evolve system through time
    … …  
    6676for t in domain.evolve(yieldstep = 0.1, finaltime = 10.0):
    6777    domain.write_time()
    68    
     78    vis.update()
    6979
     80vis.evolveFinished()
     81print 'That took %.2f seconds' %(time.time()-t0)
    7082

  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r4872 r4883  
    13271327    def extrapolate_first_order(self):
    13281328        """Extrapolate conserved quantities from centroid to
    1329         vertices for each volume using
     1329        vertices and edges for each volume using
    13301330        first order scheme.
    13311331        """
    … …  
    13331333        qc = self.centroid_values
    13341334        qv = self.vertex_values
     1335        qe = self.edge_values
    13351336
    13361337        for i in range(3):
    13371338            qv[:,i] = qc
     1339            qe[:,i] = qc
    13381340
    13391341
    … …  
    13911393        return compute_gradients(self)
    13921394
     1395
    13931396    def limit(self):
    1394         #Call correct module function
    1395         #(either from this module or C-extension)
    1396         limit(self)
     1397        #Call correct module depending on whether
     1398        #basing limit calculations on edges or vertices
     1399        limit_old(self)
     1400
     1401
     1402##     def limit_by_vertex(self):
     1403##         #Call correct module function
     1404##         #(either from this module or C-extension)
     1405##         limit_by_vertex()
     1406
     1407
     1408##     def limit_by_edge(self):
     1409##         #Call correct module function
     1410##         #(either from this module or C-extension)
     1411##         limit_by_edge()       
    13971412
    13981413    def extrapolate_second_order(self):
    … …  
    14001415        #(either from this module or C-extension)
    14011416        extrapolate_second_order(self)
     1417
    14021418
    14031419    def backup_centroid_values(self):
    … …  
    14171433    # Underlying C implementations can be accessed
    14181434
    1419     from quantity_ext import average_vertex_values, backup_centroid_values, \
    1420          saxpy_centroid_values
    1421 
    1422     from quantity_ext import compute_gradients, limit,\
    1423     extrapolate_second_order, interpolate_from_vertices_to_edges, update   
     1435    from quantity_ext import \
     1436         average_vertex_values,\
     1437         backup_centroid_values,\
     1438         saxpy_centroid_values,\
     1439         compute_gradients,\
     1440         limit_old,\
     1441         extrapolate_second_order,\
     1442         interpolate_from_vertices_to_edges,\
     1443         update   
    14241444else:
    14251445    msg = 'C implementations could not be accessed by %s.\n ' %__file__

  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r4769 r4883  
    9898                 double* vertex_coordinates,
    9999                 double* vertex_values,
     100                 double* edge_values,
    100101                 double* a,
    101102                 double* b) {
    … …  
    125126                vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
    126127                vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
     128
     129
     130                edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]);
     131                edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]);
     132                edge_values[k3+2] = 0.5*(vertex_values[k3 + 0]+vertex_values[k3 + 1]);
    127133
    128134        }
    … …  
    595601        PyObject *quantity, *domain;
    596602        PyArrayObject
    597                 *centroids,            //Coordinates at centroids
    598                 *centroid_values,      //Values at centroids
    599                 *vertex_coordinates,   //Coordinates at vertices
    600                 *vertex_values,        //Values at vertices
    601                 *number_of_boundaries, //Number of boundaries for each triangle
    602                 *surrogate_neighbours, //True neighbours or - if one missing - self
    603                 *a, *b;                //Gradients
     603            *centroids,            //Coordinates at centroids
     604            *centroid_values,      //Values at centroids
     605            *vertex_coordinates,   //Coordinates at vertices
     606            *vertex_values,        //Values at vertices
     607            *edge_values,          //Values at edges
     608            *number_of_boundaries, //Number of boundaries for each triangle
     609            *surrogate_neighbours, //True neighbours or - if one missing - self
     610            *a, *b;                //Gradients
    604611
    605612        //int N, err;
    … …  
    622629
    623630        // Get pertinent variables
    624         centroids = get_consecutive_array(domain, "centroid_coordinates");
    625         centroid_values = get_consecutive_array(quantity, "centroid_values");
    626         surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
    627         number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
    628         vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
    629         vertex_values = get_consecutive_array(quantity, "vertex_values");
     631        centroids            = get_consecutive_array(domain,   "centroid_coordinates");
     632        centroid_values      = get_consecutive_array(quantity, "centroid_values");
     633        surrogate_neighbours = get_consecutive_array(domain,   "surrogate_neighbours");
     634        number_of_boundaries = get_consecutive_array(domain,   "number_of_boundaries");
     635        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
     636        vertex_values        = get_consecutive_array(quantity, "vertex_values");
     637        edge_values          = get_consecutive_array(quantity, "edge_values");
    630638
    631639        N = centroid_values -> dimensions[0];
    … …  
    665673                        (double*) vertex_coordinates -> data,
    666674                        (double*) vertex_values -> data,
     675                        (double*) edge_values -> data,
    667676                        (double*) a -> data,
    668677                        (double*) b -> data);
    … …  
    684693        Py_DECREF(vertex_coordinates);
    685694        Py_DECREF(vertex_values);
     695        Py_DECREF(edge_values);
    686696        Py_DECREF(a);
    687697        Py_DECREF(b);
    … …  
    692702
    693703
    694 PyObject *limit(PyObject *self, PyObject *args) {
     704PyObject *limit_old(PyObject *self, PyObject *args) {
    695705  //Limit slopes for each volume to eliminate artificial variance
    696706  //introduced by e.g. second order extrapolator
    … …  
    707717        PyObject *quantity, *domain, *Tmp;
    708718        PyArrayObject
    709                 *qv, //Conserved quantities at vertices
    710                 *qc, //Conserved quantities at centroids
    711                 *neighbours;
     719            *qv, //Conserved quantities at vertices
     720            *qc, //Conserved quantities at centroids
     721            *neighbours;
    712722
    713723        int k, i, n, N, k3;
    … …  
    718728        if (!PyArg_ParseTuple(args, "O", &quantity)) {
    719729          PyErr_SetString(PyExc_RuntimeError,
    720                           "quantity_ext.c: limit could not parse input");
     730                          "quantity_ext.c: limit_old could not parse input");
    721731          return NULL;
    722732        }
    … …  
    725735        if (!domain) {
    726736          PyErr_SetString(PyExc_RuntimeError,
    727                           "quantity_ext.c: limit could not obtain domain object from quantity");                       
     737                          "quantity_ext.c: limit_old could not obtain domain object from quantity");                   
    728738         
    729739          return NULL;
    … …  
    737747        if (!Tmp) {
    738748          PyErr_SetString(PyExc_RuntimeError,
    739                           "quantity_ext.c: limit could not obtain beta_w object from domain");                 
     749                          "quantity_ext.c: limit_old could not obtain beta_w object from domain");                     
    740750         
    741751          return NULL;
    … …  
    777787
    778788        // Call underlying routine
    779         _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
     789        _limit_old(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
    780790
    781791        free(qmin);
    … …  
    788798// Method table for python module
    789799static struct PyMethodDef MethodTable[] = {
    790         {"limit", limit, METH_VARARGS, "Print out"},
     800        {"limit_old", limit_old, METH_VARARGS, "Print out"},
    791801        {"update", update, METH_VARARGS, "Print out"},
    792802        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},

  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r4779 r4883  
    107107
    108108        #print quantity.vertex_values
    109         #assert allclose(quantity.vertex_values, [[2., 2., 2.],
    110         #                                         [3.+2./3, 6.+2./3, 4.+2./3],
    111         #                                         [7.5, 0.5, 1.],
    112         #                                         [-5, -2.5, 7.5]])
    113 
    114         assert allclose(quantity.vertex_values[1,:],[3.+2./3, 6.+2./3, 4.+2./3])
    115         #FIXME: Work out the others
    116 
     109        assert allclose(quantity.vertex_values, [[3.5, -1.0, 3.5],
     110                                                 [3.+2./3, 6.+2./3, 4.+2./3],
     111                                                 [4.6, 3.4, 1.],
     112                                                 [-5.0, 1.0, 4.0]])
    117113
    118114        #print quantity.edge_values
    119         assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
    120                                                [5., 5., 5.],
    121                                                [4.5, 4.5, 0.],
    122                                                [3.0, -1.5, -1.5]])
     115        assert allclose(quantity.edge_values, [[1.25, 3.5, 1.25],
     116                                               [5. + 2/3.0, 4.0 + 1.0/6, 5.0 + 1.0/6],
     117                                               [2.2, 2.8, 4.0],
     118                                               [2.5, -0.5, -2.0]])
     119
     120
     121        #assert allclose(quantity.edge_values, [[2.5, 2.0, 1.5],
     122        #                                       [5., 5., 5.],
     123        #                                       [4.5, 4.5, 0.],
     124        #                                       [3.0, -1.5, -1.5]])
    123125
    124126    def test_get_extrema_1(self):

  • anuga_core/source/anuga/shallow_water/shallow_water_ext.c

    r4834 r4883  
    24162416
    24172417  // Call underlying standard routine
    2418   _limit(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
     2418  _limit_old(N, beta_h, (double*) hc -> data, (double*) hvbar -> data, hmin, hmax);
    24192419
    24202420  // // //Py_DECREF(domain); //FIXME: Necessary?

  • anuga_core/source/anuga/utilities/util_ext.h

    r2508 r4883  
    132132
    133133
    134 void _limit(int N, double beta, double* qc, double* qv,
     134void _limit_old(int N, double beta, double* qc, double* qv,
     135            double* qmin, double* qmax) {
     136
     137  //N are the number of elements
     138  int k, i, k3;
     139  double dq, dqa[3], phi, r;
     140 
     141  //printf("INSIDE\n");
     142  for (k=0; k<N; k++) {
     143    k3 = k*3;
     144   
     145    //Find the gradient limiter (phi) across vertices 
     146    phi = 1.0;
     147    for (i=0; i<3; i++) {   
     148      r = 1.0;
     149     
     150      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
     151      dqa[i] = dq;              //Save dq for use in the next loop
     152     
     153      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
     154      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;     
     155 
     156 
     157      phi = min( min(r*beta, 1.0), phi);   
     158    }
     159   
     160    //Then update using phi limiter
     161    for (i=0; i<3; i++) {   
     162      qv[k3+i] = qc[k] + phi*dqa[i];
     163    }
     164  }
     165}
     166
     167
     168void _limit_by_edge(int N, double beta, double* qc, double* qv,
     169            double* qmin, double* qmax) {
     170
     171  //N are the number of elements
     172  int k, i, k3;
     173  double dq, dqa[3], phi, r;
     174 
     175  //printf("INSIDE\n");
     176  for (k=0; k<N; k++) {
     177    k3 = k*3;
     178   
     179    //Find the gradient limiter (phi) across vertices 
     180    phi = 1.0;
     181    for (i=0; i<3; i++) {   
     182      r = 1.0;
     183     
     184      dq = qv[k3+i] - qc[k];    //Delta between vertex and centroid values
     185      dqa[i] = dq;              //Save dq for use in the next loop
     186     
     187      if (dq > 0.0) r = (qmax[k] - qc[k])/dq;
     188      if (dq < 0.0) r = (qmin[k] - qc[k])/dq;     
     189 
     190 
     191      phi = min( min(r*beta, 1.0), phi);   
     192    }
     193   
     194    //Then update using phi limiter
     195    for (i=0; i<3; i++) {   
     196      qv[k3+i] = qc[k] + phi*dqa[i];
     197    }
     198  }
     199}
     200
     201void _limit_by_vertex(int N, double beta, double* qc, double* qv,
    135202            double* qmin, double* qmax) {
    136203
Note: See TracChangeset for help on using the changeset viewer.