Changeset 4957


Ignore:
Timestamp:
Jan 18, 2008, 5:30:13 PM (17 years ago)
Author:
steve
Message:

Removed Conserved_quantity from domain, now just using Quantity.

In the process of adding a new limiter limit_gradient_by_neighbour to Quantity

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r4897 r4957  
    104104        if verbose: print 'Initialising Domain'
    105105        from Numeric import zeros, Float, Int, ones
    106         from quantity import Quantity, Conserved_quantity
     106        from quantity import Quantity
    107107
    108108        # List of quantity names entering
    … …  
    125125        #FIXME: remove later - maybe OK, though....
    126126        for name in self.conserved_quantities:
    127             self.quantities[name] = Conserved_quantity(self)
     127            self.quantities[name] = Quantity(self)
    128128        for name in self.other_quantities:
    129129            self.quantities[name] = Quantity(self)

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

    r4897 r4957  
    6060        self.edge_values = zeros((N, 3), Float)
    6161
     62        #Allocate space for Gradient
     63        self.x_gradient = zeros(N, Float)
     64        self.y_gradient = zeros(N, Float)       
     65
    6266        #Intialise centroid and edge_values
    6367        self.interpolate()
     68
     69        #Allocate space for boundary values
     70        L = len(domain.boundary)
     71        self.boundary_values = zeros(L, Float)
     72
     73        #Allocate space for updates of conserved quantities by
     74        #flux calculations and forcing functions
     75
     76        N = len(domain) # number_of_triangles
     77        self.explicit_update = zeros(N, Float )
     78        self.semi_implicit_update = zeros(N, Float )
     79        self.centroid_backup_values = zeros(N, Float)       
    6480
    6581
    … …  
    13301346            qe[:,i] = qc
    13311347
     1348        self.x_gradient *= 0.0
     1349        self.y_gradient *= 0.0
     1350
    13321351
    13331352    def get_integral(self):
    … …  
    13421361        return integral
    13431362
    1344 
    1345 
    1346 
    1347 class Conserved_quantity(Quantity):
    1348     """Class conserved quantity adds to Quantity:
    1349 
    1350     boundary values, storage and method for updating, and
    1351     methods for (second order) extrapolation from centroid to vertices inluding
    1352     gradients and limiters
    1353     """
    1354 
    1355     def __init__(self, domain, vertex_values=None):
    1356         Quantity.__init__(self, domain, vertex_values)
    1357 
    1358         from Numeric import zeros, Float
    1359 
    1360         #Allocate space for boundary values
    1361         L = len(domain.boundary)
    1362         self.boundary_values = zeros(L, Float)
    1363 
    1364         #Allocate space for updates of conserved quantities by
    1365         #flux calculations and forcing functions
    1366 
    1367         N = len(domain) # number_of_triangles
    1368         self.explicit_update = zeros(N, Float )
    1369         self.semi_implicit_update = zeros(N, Float )
    1370         self.centroid_backup_values = zeros(N, Float)       
    1371 
    1372        
     1363    def get_gradients(self):
     1364        """Provide gradients. Use compute_gradients first
     1365        """
     1366
     1367        return self.x_gradient, self.y_gradient
    13731368
    13741369
    … …  
    14061401        #Call correct module function
    14071402        #(either from this module or C-extension)
    1408         extrapolate_second_order(self)
     1403        compute_gradients(self)
     1404        extrapolate_from_gradient(self)
    14091405
    14101406    def backup_centroid_values(self):
    … …  
    14181414        saxpy_centroid_values(self,a,b)
    14191415   
     1416
     1417
     1418class Conserved_quantity(Quantity):
     1419    """Class conserved quantity being removed, use Quantity
     1420
     1421    """
     1422
     1423    def __init__(self, domain, vertex_values=None):
     1424        #Quantity.__init__(self, domain, vertex_values)
     1425
     1426        msg = 'ERROR: Use Quantity instead of Conserved_quantity'
     1427
     1428        raise Exception, msg
     1429
    14201430
    14211431
    … …  
    14331443         limit_edges_by_all_neighbours,\
    14341444         limit_edges_by_neighbour,\
    1435          extrapolate_second_order,\
     1445         limit_gradient_by_neighbour,\
     1446         extrapolate_from_gradient,\
    14361447         interpolate_from_vertices_to_edges,\
    14371448         interpolate_from_edges_to_vertices,\

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

    r4897 r4957  
    9393
    9494
    95 int _extrapolate(int N,
     95int _extrapolate_from_gradient(int N,
    9696                 double* centroids,
    9797                 double* centroid_values,
    … …  
    258258                     double* vertex_values,
    259259                     double* edge_values,
     260                     long*   neighbours) {
     261
     262        int i, k, k2, k3, k6;
     263        long n;
     264        double qmin, qmax, qn, qc;
     265        double dq, dqa[3], phi, r;
     266
     267        for (k=0; k<N; k++){
     268                k6 = 6*k;
     269                k3 = 3*k;
     270                k2 = 2*k;
     271
     272                qc = centroid_values[k];
     273                phi = 1.0;
     274
     275                for (i=0; i<3; i++) {
     276                    dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values
     277                    dqa[i] = dq;                      //Save dq for use in updating vertex values
     278
     279                    n = neighbours[k3+i];
     280                    if (n >= 0) {
     281                        qn = centroid_values[n]; //Neighbour's centroid value
     282
     283                        qmin = min(qc, qn);
     284                        qmax = max(qc, qn);
     285
     286                        r = 1.0;
     287     
     288                        if (dq > 0.0) r = (qmax - qc)/dq;
     289                        if (dq < 0.0) r = (qmin - qc)/dq;     
     290                   
     291                        phi = min( min(r*beta, 1.0), phi);   
     292                    }
     293                }
     294
     295
     296                //Update edge and vertex values using phi limiter
     297                edge_values[k3+0] = qc + phi*dqa[0];
     298                edge_values[k3+1] = qc + phi*dqa[1];
     299                edge_values[k3+2] = qc + phi*dqa[2];
     300               
     301                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
     302                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
     303                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     304
     305        }
     306
     307        return 0;
     308}
     309
     310
     311int _limit_gradient_by_neighbour(int N, double beta,
     312                     double* centroid_values,
     313                     double* vertex_values,
     314                     double* edge_values,
     315                     double* x_gradient,
     316                     double* y_gradient,
    260317                     long*   neighbours) {
    261318
    … …  
    752809
    753810
    754 PyObject *compute_gradients(PyObject *self, PyObject *args) {
    755   //"""Compute gradients of triangle surfaces defined by centroids of
    756   //neighbouring volumes.
    757   //If one edge is on the boundary, use own centroid as neighbour centroid.
    758   //If two or more are on the boundary, fall back to first order scheme.
    759   //"""
    760 
    761 
    762         PyObject *quantity, *domain, *R;
    763         PyArrayObject
    764                 *centroids,            //Coordinates at centroids
    765                 *centroid_values,      //Values at centroids
    766                 *number_of_boundaries, //Number of boundaries for each triangle
    767                 *surrogate_neighbours, //True neighbours or - if one missing - self
    768                 *a, *b;                //Return values
    769 
    770         int dimensions[1], N, err;
    771 
    772         // Convert Python arguments to C
    773         if (!PyArg_ParseTuple(args, "O", &quantity)) {
    774           PyErr_SetString(PyExc_RuntimeError,
    775                           "quantity_ext.c: compute_gradients could not parse input");   
    776           return NULL;
    777         }
    778 
    779         domain = PyObject_GetAttrString(quantity, "domain");
    780         if (!domain) {
    781           PyErr_SetString(PyExc_RuntimeError,
    782                           "compute_gradients could not obtain domain object from quantity");
    783           return NULL;
    784         }
    785 
    786         // Get pertinent variables
    787 
    788         centroids = get_consecutive_array(domain, "centroid_coordinates");
    789         centroid_values = get_consecutive_array(quantity, "centroid_values");
    790         surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
    791         number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
    792 
    793         N = centroid_values -> dimensions[0];
    794 
    795         // Release
    796         Py_DECREF(domain);
    797 
    798         // Allocate space for return vectors a and b (don't DECREF)
    799         dimensions[0] = N;
    800         a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    801         b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    802 
    803 
    804 
    805         err = _compute_gradients(N,
    806                         (double*) centroids -> data,
    807                         (double*) centroid_values -> data,
    808                         (long*) number_of_boundaries -> data,
    809                         (long*) surrogate_neighbours -> data,
    810                         (double*) a -> data,
    811                         (double*) b -> data);
    812 
    813         if (err != 0) {
    814           PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    815           return NULL;
    816         }
    817 
    818         // Release
    819         Py_DECREF(centroids);
    820         Py_DECREF(centroid_values);
    821         Py_DECREF(number_of_boundaries);
    822         Py_DECREF(surrogate_neighbours);
    823 
    824         // Build result, release and return
    825         R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
    826         Py_DECREF(a);
    827         Py_DECREF(b);
    828         return R;
    829 }
    830 
    831 
    832 
    833 PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
     811/* PyObject *compute_gradients(PyObject *self, PyObject *args) { */
     812/*   //"""Compute gradients of triangle surfaces defined by centroids of */
     813/*   //neighbouring volumes. */
     814/*   //If one edge is on the boundary, use own centroid as neighbour centroid. */
     815/*   //If two or more are on the boundary, fall back to first order scheme. */
     816/*   //""" */
     817
     818
     819/*      PyObject *quantity, *domain, *R; */
     820/*      PyArrayObject */
     821/*              *centroids,            //Coordinates at centroids */
     822/*              *centroid_values,      //Values at centroids */
     823/*              *number_of_boundaries, //Number of boundaries for each triangle */
     824/*              *surrogate_neighbours, //True neighbours or - if one missing - self */
     825/*              *a, *b;                //Return values */
     826
     827/*      int dimensions[1], N, err; */
     828
     829/*      // Convert Python arguments to C */
     830/*      if (!PyArg_ParseTuple(args, "O", &quantity)) { */
     831/*        PyErr_SetString(PyExc_RuntimeError,  */
     832/*                        "quantity_ext.c: compute_gradients could not parse input");    */
     833/*        return NULL; */
     834/*      } */
     835
     836/*      domain = PyObject_GetAttrString(quantity, "domain"); */
     837/*      if (!domain) { */
     838/*        PyErr_SetString(PyExc_RuntimeError,  */
     839/*                        "compute_gradients could not obtain domain object from quantity"); */
     840/*        return NULL; */
     841/*      } */
     842
     843/*      // Get pertinent variables */
     844
     845/*      centroids = get_consecutive_array(domain, "centroid_coordinates"); */
     846/*      centroid_values = get_consecutive_array(quantity, "centroid_values"); */
     847/*      surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours"); */
     848/*      number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries"); */
     849
     850/*      N = centroid_values -> dimensions[0]; */
     851
     852/*      // Release */
     853/*      Py_DECREF(domain); */
     854
     855/*      // Allocate space for return vectors a and b (don't DECREF) */
     856/*      dimensions[0] = N; */
     857/*      a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
     858/*      b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); */
     859
     860
     861
     862/*      err = _compute_gradients(N, */
     863/*                      (double*) centroids -> data, */
     864/*                      (double*) centroid_values -> data, */
     865/*                      (long*) number_of_boundaries -> data, */
     866/*                      (long*) surrogate_neighbours -> data, */
     867/*                      (double*) a -> data, */
     868/*                      (double*) b -> data); */
     869
     870/*      if (err != 0) { */
     871/*        PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); */
     872/*        return NULL; */
     873/*      } */
     874
     875/*      // Release */
     876/*      Py_DECREF(centroids); */
     877/*      Py_DECREF(centroid_values); */
     878/*      Py_DECREF(number_of_boundaries); */
     879/*      Py_DECREF(surrogate_neighbours); */
     880
     881/*      // Build result, release and return */
     882/*      R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); */
     883/*      Py_DECREF(a); */
     884/*      Py_DECREF(b); */
     885/*      return R; */
     886/* } */
     887
     888
     889
     890PyObject *extrapolate_from_gradient(PyObject *self, PyObject *args) {
    834891
    835892        PyObject *quantity, *domain;
    … …  
    842899            *number_of_boundaries, //Number of boundaries for each triangle
    843900            *surrogate_neighbours, //True neighbours or - if one missing - self
    844             *a, *b;                //Gradients
     901            *x_gradient,           //x gradient
     902            *y_gradient;           //y gradient
    845903
    846904        //int N, err;
    847         int dimensions[1], N, err;
     905        //int dimensions[1];
     906        int N, err;
    848907        //double *a, *b;  //Gradients
    849908
    … …  
    851910        if (!PyArg_ParseTuple(args, "O", &quantity)) {
    852911          PyErr_SetString(PyExc_RuntimeError,
    853                           "extrapolate_second_order could not parse input");   
     912                          "extrapolate_gradient could not parse input");       
    854913          return NULL;
    855914        }
    … …  
    858917        if (!domain) {
    859918          PyErr_SetString(PyExc_RuntimeError,
    860                           "extrapolate_second_order could not obtain domain object from quantity");     
     919                          "extrapolate_gradient could not obtain domain object from quantity");
    861920          return NULL;
    862921        }
    … …  
    870929        vertex_values        = get_consecutive_array(quantity, "vertex_values");
    871930        edge_values          = get_consecutive_array(quantity, "edge_values");
     931        x_gradient           = get_consecutive_array(quantity, "x_gradient");
     932        y_gradient           = get_consecutive_array(quantity, "y_gradient");
    872933
    873934        N = centroid_values -> dimensions[0];
    … …  
    877938
    878939        //Allocate space for return vectors a and b (don't DECREF)
    879         dimensions[0] = N;
    880         a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    881         b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     940        //dimensions[0] = N;
     941        //a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     942        //b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
    882943
    883944        //FIXME: Odd that I couldn't use normal arrays
    … …  
    889950
    890951
    891         err = _compute_gradients(N,
    892                         (double*) centroids -> data,
    893                         (double*) centroid_values -> data,
    894                         (long*) number_of_boundaries -> data,
    895                         (long*) surrogate_neighbours -> data,
    896                         (double*) a -> data,
    897                         (double*) b -> data);
    898 
    899         if (err != 0) {
    900           PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
    901           return NULL;
    902         }
    903 
    904         err = _extrapolate(N,
     952/*      err = _compute_gradients(N, */
     953/*                      (double*) centroids -> data, */
     954/*                      (double*) centroid_values -> data, */
     955/*                      (long*) number_of_boundaries -> data, */
     956/*                      (long*) surrogate_neighbours -> data, */
     957/*                      (double*) x_gradient -> data, */
     958/*                      (double*) y_gradient -> data); */
     959
     960/*      if (err != 0) { */
     961/*        PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); */
     962/*        return NULL; */
     963/*      } */
     964
     965        err = _extrapolate_from_gradient(N,
    905966                        (double*) centroids -> data,
    906967                        (double*) centroid_values -> data,
    … …  
    908969                        (double*) vertex_values -> data,
    909970                        (double*) edge_values -> data,
    910                         (double*) a -> data,
    911                         (double*) b -> data);
     971                        (double*) x_gradient -> data,
     972                        (double*) y_gradient -> data);
    912973
    913974
    … …  
    928989        Py_DECREF(vertex_values);
    929990        Py_DECREF(edge_values);
    930         Py_DECREF(a);
    931         Py_DECREF(b);
     991        Py_DECREF(x_gradient);
     992        Py_DECREF(y_gradient);
     993
     994        return Py_BuildValue("");
     995}
     996
     997
     998
     999PyObject *compute_gradients(PyObject *self, PyObject *args) {
     1000
     1001        PyObject *quantity, *domain;
     1002        PyArrayObject
     1003            *centroids,            //Coordinates at centroids
     1004            *centroid_values,      //Values at centroids
     1005            *vertex_coordinates,   //Coordinates at vertices
     1006            *vertex_values,        //Values at vertices
     1007            *edge_values,          //Values at edges
     1008            *number_of_boundaries, //Number of boundaries for each triangle
     1009            *surrogate_neighbours, //True neighbours or - if one missing - self
     1010            *x_gradient,           //x gradient
     1011            *y_gradient;           //y gradient
     1012
     1013        //int N, err;
     1014        //int dimensions[1];
     1015        int N, err;
     1016        //double *a, *b;  //Gradients
     1017
     1018        // Convert Python arguments to C
     1019        if (!PyArg_ParseTuple(args, "O", &quantity)) {
     1020          PyErr_SetString(PyExc_RuntimeError,
     1021                          "compute_gradients could not parse input");   
     1022          return NULL;
     1023        }
     1024
     1025        domain = PyObject_GetAttrString(quantity, "domain");
     1026        if (!domain) {
     1027          PyErr_SetString(PyExc_RuntimeError,
     1028                          "compute_gradients could not obtain domain object from quantity");   
     1029          return NULL;
     1030        }
     1031
     1032        // Get pertinent variables
     1033        centroids            = get_consecutive_array(domain,   "centroid_coordinates");
     1034        centroid_values      = get_consecutive_array(quantity, "centroid_values");
     1035        surrogate_neighbours = get_consecutive_array(domain,   "surrogate_neighbours");
     1036        number_of_boundaries = get_consecutive_array(domain,   "number_of_boundaries");
     1037        vertex_coordinates   = get_consecutive_array(domain,   "vertex_coordinates");
     1038        vertex_values        = get_consecutive_array(quantity, "vertex_values");
     1039        edge_values          = get_consecutive_array(quantity, "edge_values");
     1040        x_gradient           = get_consecutive_array(quantity, "x_gradient");
     1041        y_gradient           = get_consecutive_array(quantity, "y_gradient");
     1042
     1043        N = centroid_values -> dimensions[0];
     1044
     1045        // Release
     1046        Py_DECREF(domain);
     1047
     1048        //Allocate space for return vectors a and b (don't DECREF)
     1049        //dimensions[0] = N;
     1050        //a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     1051        //b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
     1052
     1053        //FIXME: Odd that I couldn't use normal arrays
     1054        //Allocate space for return vectors a and b (don't DECREF)
     1055        //a = (double*) malloc(N * sizeof(double));
     1056        //if (!a) return NULL;
     1057        //b = (double*) malloc(N * sizeof(double));
     1058        //if (!b) return NULL;
     1059
     1060
     1061        err = _compute_gradients(N,
     1062                        (double*) centroids -> data,
     1063                        (double*) centroid_values -> data,
     1064                        (long*) number_of_boundaries -> data,
     1065                        (long*) surrogate_neighbours -> data,
     1066                        (double*) x_gradient -> data,
     1067                        (double*) y_gradient -> data);
     1068
     1069        if (err != 0) {
     1070          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
     1071          return NULL;
     1072        }
     1073
     1074/*      err = _extrapolate(N, */
     1075/*                      (double*) centroids -> data, */
     1076/*                      (double*) centroid_values -> data, */
     1077/*                      (double*) vertex_coordinates -> data, */
     1078/*                      (double*) vertex_values -> data, */
     1079/*                      (double*) edge_values -> data, */
     1080/*                      (double*) x_gradient -> data, */
     1081/*                      (double*) y_gradient -> data); */
     1082
     1083
     1084/*      if (err != 0) { */
     1085/*        PyErr_SetString(PyExc_RuntimeError, */
     1086/*                        "Internal function _extrapolate failed"); */
     1087/*        return NULL; */
     1088/*      } */
     1089
     1090
     1091
     1092        // Release
     1093        Py_DECREF(centroids);
     1094        Py_DECREF(centroid_values);
     1095        Py_DECREF(number_of_boundaries);
     1096        Py_DECREF(surrogate_neighbours);
     1097        Py_DECREF(vertex_coordinates);
     1098        Py_DECREF(vertex_values);
     1099        Py_DECREF(edge_values);
     1100        Py_DECREF(x_gradient);
     1101        Py_DECREF(y_gradient);
    9321102
    9331103        return Py_BuildValue("");
    … …  
    12821452
    12831453
     1454PyObject *limit_gradient_by_neighbour(PyObject *self, PyObject *args) {
     1455  //Limit slopes for each volume to eliminate artificial variance
     1456  //introduced by e.g. second order extrapolator
     1457
     1458  //This is an unsophisticated limiter as it does not take into
     1459  //account dependencies among quantities.
     1460
     1461  //precondition:
     1462  //  vertex values are estimated from gradient
     1463  //postcondition:
     1464  //  vertex and edge values are updated
     1465  //
     1466
     1467        PyObject *quantity, *domain, *Tmp;
     1468        PyArrayObject
     1469            *vertex_values,   //Conserved quantities at vertices
     1470            *centroid_values, //Conserved quantities at centroids
     1471            *edge_values,     //Conserved quantities at edges
     1472            *x_gradient,
     1473            *y_gradient,
     1474            *neighbours;
     1475
     1476        double beta_w; //Safety factor
     1477        int N, err;
     1478
     1479
     1480        // Convert Python arguments to C
     1481        if (!PyArg_ParseTuple(args, "O", &quantity)) {
     1482          PyErr_SetString(PyExc_RuntimeError,
     1483                          "quantity_ext.c: limit_gradient_by_neighbour could not parse input");
     1484          return NULL;
     1485        }
     1486
     1487        domain = PyObject_GetAttrString(quantity, "domain");
     1488        if (!domain) {
     1489          PyErr_SetString(PyExc_RuntimeError,
     1490                          "quantity_ext.c: limit_gradient_by_neighbour could not obtain domain object from quantity");                 
     1491         
     1492          return NULL;
     1493        }
     1494
     1495        // Get safety factor beta_w
     1496        Tmp = PyObject_GetAttrString(domain, "beta_w");
     1497        if (!Tmp) {
     1498          PyErr_SetString(PyExc_RuntimeError,
     1499                          "quantity_ext.c: limit_gradient_by_neighbour could not obtain beta_w object from domain");                   
     1500         
     1501          return NULL;
     1502        }       
     1503
     1504
     1505        // Get pertinent variables
     1506        neighbours       = get_consecutive_array(domain, "neighbours");
     1507        centroid_values  = get_consecutive_array(quantity, "centroid_values");
     1508        vertex_values    = get_consecutive_array(quantity, "vertex_values");
     1509        edge_values      = get_consecutive_array(quantity, "edge_values");
     1510        x_gradient       = get_consecutive_array(quantity, "x_gradient");
     1511        y_gradient       = get_consecutive_array(quantity, "y_gradient");
     1512
     1513        beta_w           = PyFloat_AsDouble(Tmp);
     1514
     1515
     1516        N = centroid_values -> dimensions[0];
     1517
     1518        err = _limit_gradient_by_neighbour(N, beta_w,
     1519                                        (double*) centroid_values -> data,
     1520                                        (double*) vertex_values -> data,
     1521                                        (double*) edge_values -> data,
     1522                                        (double*) x_gradient -> data,
     1523                                        (double*) y_gradient -> data,
     1524                                        (long*)   neighbours -> data);
     1525       
     1526        if (err != 0) {
     1527          PyErr_SetString(PyExc_RuntimeError,
     1528                          "Internal function _limit_gradient_by_neighbour failed");
     1529          return NULL;
     1530        }       
     1531
     1532
     1533        // Release
     1534        Py_DECREF(neighbours);
     1535        Py_DECREF(centroid_values);
     1536        Py_DECREF(vertex_values);
     1537        Py_DECREF(edge_values);
     1538        Py_DECREF(x_gradient);
     1539        Py_DECREF(y_gradient);
     1540        Py_DECREF(Tmp);
     1541
     1542
     1543        return Py_BuildValue("");
     1544}
     1545
     1546
    12841547// Method table for python module
    12851548static struct PyMethodDef MethodTable[] = {
    … …  
    12881551        {"limit_edges_by_all_neighbours", limit_edges_by_all_neighbours, METH_VARARGS, "Print out"},
    12891552        {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"},
     1553        {"limit_gradient_by_neighbour", limit_gradient_by_neighbour, METH_VARARGS, "Print out"},
    12901554        {"update", update, METH_VARARGS, "Print out"},
    12911555        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
    12921556        {"saxpy_centroid_values", saxpy_centroid_values, METH_VARARGS, "Print out"},
    12931557        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
    1294         {"extrapolate_second_order", extrapolate_second_order,
     1558        {"extrapolate_from_gradient", extrapolate_from_gradient,
    12951559                METH_VARARGS, "Print out"},
    12961560        {"interpolate_from_vertices_to_edges",

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

    r4897 r4957  
    9999
    100100    def test_interpolation2(self):
    101         quantity = Conserved_quantity(self.mesh4,
     101        quantity = Quantity(self.mesh4,
    102102                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    103103        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroid
    … …  
    125125
    126126    def test_get_extrema_1(self):
    127         quantity = Conserved_quantity(self.mesh4,
     127        quantity = Quantity(self.mesh4,
    128128                                      [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    129129        assert allclose(quantity.centroid_values, [2., 5., 3., 0.]) #Centroids
    … …  
    231231
    232232    def test_boundary_allocation(self):
    233         quantity = Conserved_quantity(self.mesh4,
     233        quantity = Quantity(self.mesh4,
    234234                            [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])
    235235
    … …  
    12891289
    12901290
    1291     def test_gradient(self):
    1292         quantity = Conserved_quantity(self.mesh4)
     1291    def test_compute_gradient(self):
     1292        quantity = Quantity(self.mesh4)
    12931293
    12941294        #Set up for a gradient of (2,0) at mid triangle
    … …  
    12981298
    12991299        #Gradients
    1300         a, b = quantity.compute_gradients()
    1301 
     1300        quantity.compute_gradients()
     1301
     1302        a = quantity.x_gradient
     1303        b = quantity.y_gradient
    13021304        #print self.mesh4.centroid_coordinates
    13031305        #print a, b
    … …  
    13461348        assert allclose(quantity.vertex_values[2,2], 8)
    13471349
     1350    def test_get_gradients(self):
     1351        quantity = Quantity(self.mesh4)
     1352
     1353        #Set up for a gradient of (2,0) at mid triangle
     1354        quantity.set_values([2.0, 4.0, 6.0, 2.0],
     1355                            location = 'centroids')
     1356
     1357
     1358        #Gradients
     1359        quantity.compute_gradients()
     1360
     1361        a, b = quantity.get_gradients()
     1362        #print self.mesh4.centroid_coordinates
     1363        #print a, b
     1364
     1365        #The central triangle (1)
     1366        #(using standard gradient based on neigbours controid values)
     1367        assert allclose(a[1], 2.0)
     1368        assert allclose(b[1], 0.0)
     1369
     1370
     1371        #Left triangle (0) using two point gradient
     1372        #q0 = q1 + a*(x0-x1) + b*(y0-y1)  <=>
     1373        #2  = 4  + a*(-2/3)  + b*(-2/3)
     1374        assert allclose(a[0] + b[0], 3)
     1375        #From orthogonality (a*(y0-y1) + b*(x0-x1) == 0)
     1376        assert allclose(a[0] - b[0], 0)
     1377
     1378
     1379        #Right triangle (2) using two point gradient
     1380        #q2 = q1 + a*(x2-x1) + b*(y2-y1)  <=>
     1381        #6  = 4  + a*(4/3)  + b*(-2/3)
     1382        assert allclose(2*a[2] - b[2], 3)
     1383        #From orthogonality (a*(y1-y2) + b*(x2-x1) == 0)
     1384        assert allclose(a[2] + 2*b[2], 0)
     1385
     1386
     1387        #Top triangle (3) using two point gradient
     1388        #q3 = q1 + a*(x3-x1) + b*(y3-y1)  <=>
     1389        #2  = 4  + a*(-2/3)  + b*(4/3)
     1390        assert allclose(a[3] - 2*b[3], 3)
     1391        #From orthogonality (a*(y1-y3) + b*(x3-x1) == 0)
     1392        assert allclose(2*a[3] + b[3], 0)
     1393
    13481394
    13491395    def test_second_order_extrapolation2(self):
    1350         quantity = Conserved_quantity(self.mesh4)
     1396        quantity = Quantity(self.mesh4)
    13511397
    13521398        #Set up for a gradient of (3,1), f(x) = 3x+y
    … …  
    13551401
    13561402        #Gradients
    1357         a, b = quantity.compute_gradients()
    1358 
     1403        quantity.compute_gradients()
     1404
     1405        a = quantity.x_gradient
     1406        b = quantity.y_gradient
     1407       
    13591408        #print a, b
    13601409
    … …  
    13741423
    13751424    def test_backup_saxpy_centroid_values(self):
    1376         quantity = Conserved_quantity(self.mesh4)
     1425        quantity = Quantity(self.mesh4)
    13771426
    13781427        #Set up for a gradient of (3,1), f(x) = 3x+y
    … …  
    13971446
    13981447    def test_first_order_extrapolator(self):
    1399         quantity = Conserved_quantity(self.mesh4)
     1448        quantity = Quantity(self.mesh4)
    14001449
    14011450        #Test centroids
    … …  
    14061455        quantity.extrapolate_first_order()
    14071456
     1457        #Check that gradient is zero
     1458        a,b = quantity.get_gradients()
     1459        assert allclose(a, [0,0,0,0])
     1460        assert allclose(b, [0,0,0,0])
     1461
    14081462        #Check vertices but not edge values
    14091463        assert allclose(quantity.vertex_values,
    … …  
    14121466
    14131467    def test_second_order_extrapolator(self):
    1414         quantity = Conserved_quantity(self.mesh4)
     1468        quantity = Quantity(self.mesh4)
    14151469
    14161470        #Set up for a gradient of (3,0) at mid triangle
    … …  
    14461500
    14471501    def test_limit_vertices_by_all_neighbours(self):
    1448         quantity = Conserved_quantity(self.mesh4)
     1502        quantity = Quantity(self.mesh4)
    14491503
    14501504        #Create a deliberate overshoot (e.g. from gradient computation)
    … …  
    14761530
    14771531    def test_limit_edges_by_all_neighbours(self):
    1478         quantity = Conserved_quantity(self.mesh4)
     1532        quantity = Quantity(self.mesh4)
    14791533
    14801534        #Create a deliberate overshoot (e.g. from gradient computation)
    … …  
    15051559
    15061560    def test_limit_edges_by_neighbour(self):
    1507         quantity = Conserved_quantity(self.mesh4)
     1561        quantity = Quantity(self.mesh4)
    15081562
    15091563        #Create a deliberate overshoot (e.g. from gradient computation)
    … …  
    15351589        """Taken from test_shallow_water
    15361590        """
    1537         quantity = Conserved_quantity(self.mesh4)
     1591        quantity = Quantity(self.mesh4)
    15381592        quantity.domain.beta_w = 0.9
    15391593       
    … …  
    15691623
    15701624    def test_distribute_first_order(self):
    1571         quantity = Conserved_quantity(self.mesh4)
     1625        quantity = Quantity(self.mesh4)
    15721626
    15731627        #Test centroids
    … …  
    16171671
    16181672    def test_distribute_second_order(self):
    1619         quantity = Conserved_quantity(self.mesh4)
     1673        quantity = Quantity(self.mesh4)
    16201674
    16211675        #Test centroids
    … …  
    16311685
    16321686    def test_update_explicit(self):
    1633         quantity = Conserved_quantity(self.mesh4)
     1687        quantity = Quantity(self.mesh4)
    16341688
    16351689        #Test centroids
    … …  
    16471701
    16481702    def test_update_semi_implicit(self):
    1649         quantity = Conserved_quantity(self.mesh4)
     1703        quantity = Quantity(self.mesh4)
    16501704
    16511705        #Test centroids
    … …  
    16681722
    16691723    def test_both_updates(self):
    1670         quantity = Conserved_quantity(self.mesh4)
     1724        quantity = Quantity(self.mesh4)
    16711725
    16721726        #Test centroids

  • anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

    r4868 r4957  
    23042304        domain.set_quantity('stage', stage, location='centroids')
    23052305
    2306         a, b = domain.quantities['stage'].compute_gradients()
     2306        domain.quantities['stage'].compute_gradients()
     2307
     2308        a, b = domain.quantities['stage'].get_gradients()
     2309               
    23072310        assert allclose(a[1], 3.33333334)
    23082311        assert allclose(b[1], 0.0)
    … …  
    23482351        #print domain.quantities['stage'].centroid_values
    23492352
    2350         a, b = domain.quantities['stage'].compute_gradients()
     2353        domain.quantities['stage'].compute_gradients()
     2354        a, b = domain.quantities['stage'].get_gradients()
    23512355        assert allclose(a[1], 25.18518519)
    23522356        assert allclose(b[1], 3.33333333)
Note: See TracChangeset for help on using the changeset viewer.