Changeset 4897


Ignore:
Timestamp:
Dec 27, 2007, 6:35:08 PM (17 years ago)
Author:
steve
Message:

Working on limiters

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

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

    r4836 r4897  
    380380            Depth = domain.create_quantity_from_expression('stage-elevation')
    381381
    382             exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5')       
     382            exp = '(xmomentum*xmomentum + ymomentum*ymomentum)**0.5'       
    383383            Absolute_momentum = domain.create_quantity_from_expression(exp)
    384384
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

    r4886 r4897  
    250250                   polygon = None,
    251251                   indices = None,
     252                   smooth = False,
    252253                   verbose = False,
    253254                   use_cache = False):
     
    389390            self.set_values_from_constant(numeric,
    390391                                          location, indices, verbose)
    391            
     392
     393
    392394            self.extrapolate_first_order()
     395
     396            if smooth:
     397                self.smooth_vertex_values()
     398
     399               
    393400            return
    394401       
     
    11811188
    11821189
    1183     def smooth_vertex_values(self, value_array='field_values',
    1184                              precision = None):
    1185         """ Smooths field_values or conserved_quantities data.
    1186         TODO: be able to smooth individual fields
    1187         NOTE:  This function does not have a test.
    1188         FIXME: NOT DONE - do we need it?
    1189         FIXME: this function isn't called by anything.
    1190                Maybe it should be removed..-DSG
    1191         """
    1192 
    1193         from Numeric import concatenate, zeros, Float, Int, array, reshape
    1194 
    1195 
    1196         A,V = self.get_vertex_values(xy=False,
    1197                                      value_array=value_array,
    1198                                      smooth = True,
    1199                                      precision = precision)
    1200 
    1201         #Set some field values
    1202         for volume in self:
    1203             for i,v in enumerate(volume.vertices):
    1204                 if value_array == 'field_values':
    1205                     volume.set_field_values('vertex', i, A[v,:])
    1206                 elif value_array == 'conserved_quantities':
    1207                     volume.set_conserved_quantities('vertex', i, A[v,:])
    1208 
    1209         if value_array == 'field_values':
    1210             self.precompute()
    1211         elif value_array == 'conserved_quantities':
    1212             Volume.interpolate_conserved_quantities()
     1190    def smooth_vertex_values(self):
     1191        """ Smooths vertex values.
     1192        """
     1193
     1194        A,V = self.get_vertex_values(xy=False, smooth = True)
     1195        self.set_vertex_values(A)
    12131196
    12141197
     
    12541237        if smooth is None:
    12551238            # Take default from domain
    1256             smooth = self.domain.smooth
     1239            try:
     1240                smooth = self.domain.smooth
     1241            except:
     1242                smooth = False
    12571243
    12581244        if precision is None:
     
    13921378        return update(self, timestep)
    13931379
    1394 
    13951380    def compute_gradients(self):
    13961381        #Call correct module function
     
    13981383        return compute_gradients(self)
    13991384
    1400 
    14011385    def limit(self):
    14021386        #Call correct module depending on whether
     
    14041388        limit_old(self)
    14051389
    1406 
    1407     def limit_by_vertex(self):
     1390    def limit_vertices_by_all_neighbours(self):
    14081391        #Call correct module function
    14091392        #(either from this module or C-extension)
    1410         limit_by_vertex(self)
    1411 
    1412 
    1413     def limit_by_edge(self):
     1393        limit_vertices_by_all_neighbours(self)
     1394
     1395    def limit_edges_by_all_neighbours(self):
    14141396        #Call correct module function
    14151397        #(either from this module or C-extension)
    1416         limit_by_edge(self)       
     1398        limit_edges_by_all_neighbours(self)
     1399
     1400    def limit_edges_by_neighbour(self):
     1401        #Call correct module function
     1402        #(either from this module or C-extension)
     1403        limit_edges_by_neighbour(self)               
    14171404
    14181405    def extrapolate_second_order(self):
     
    14201407        #(either from this module or C-extension)
    14211408        extrapolate_second_order(self)
    1422 
    14231409
    14241410    def backup_centroid_values(self):
     
    14441430         compute_gradients,\
    14451431         limit_old,\
    1446          limit_by_vertex,\
    1447          limit_by_edge,\
     1432         limit_vertices_by_all_neighbours,\
     1433         limit_edges_by_all_neighbours,\
     1434         limit_edges_by_neighbour,\
    14481435         extrapolate_second_order,\
    14491436         interpolate_from_vertices_to_edges,\
  • anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

    r4886 r4897  
    137137
    138138
    139 int _limit_by_vertex(int N, double beta,
     139int _limit_vertices_by_all_neighbours(int N, double beta,
    140140                     double* centroid_values,
    141141                     double* vertex_values,
     
    195195}
    196196
    197 int _limit_by_edge(int N, double beta,
     197int _limit_edges_by_all_neighbours(int N, double beta,
    198198                     double* centroid_values,
    199199                     double* vertex_values,
     
    240240   
    241241                //Update vertex and edge values using phi limiter
     242                edge_values[k3+0] = qc + phi*dqa[0];
     243                edge_values[k3+1] = qc + phi*dqa[1];
     244                edge_values[k3+2] = qc + phi*dqa[2];
     245               
     246                vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0];
     247                vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1];
     248                vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2];
     249
     250        }
     251
     252        return 0;
     253}
     254
     255
     256int _limit_edges_by_neighbour(int N, double beta,
     257                     double* centroid_values,
     258                     double* vertex_values,
     259                     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
    242297                edge_values[k3+0] = qc + phi*dqa[0];
    243298                edge_values[k3+1] = qc + phi*dqa[1];
     
    9741029
    9751030
    976 PyObject *limit_by_vertex(PyObject *self, PyObject *args) {
     1031PyObject *limit_vertices_by_all_neighbours(PyObject *self, PyObject *args) {
    9771032  //Limit slopes for each volume to eliminate artificial variance
    9781033  //introduced by e.g. second order extrapolator
     
    10331088        N = centroid_values -> dimensions[0];
    10341089
    1035         err = _limit_by_vertex(N, beta_w,
    1036                                (double*) centroid_values -> data,
    1037                                (double*) vertex_values -> data,
    1038                                (double*) edge_values -> data,
    1039                                (long*)   neighbours -> data);
    1040 
     1090        err = _limit_vertices_by_all_neighbours(N, beta_w,
     1091                                                (double*) centroid_values -> data,
     1092                                                (double*) vertex_values -> data,
     1093                                                (double*) edge_values -> data,
     1094                                                (long*)   neighbours -> data);
     1095       
    10411096        if (err != 0) {
    10421097          PyErr_SetString(PyExc_RuntimeError,
     
    10591114
    10601115
    1061 PyObject *limit_by_edge(PyObject *self, PyObject *args) {
     1116PyObject *limit_edges_by_all_neighbours(PyObject *self, PyObject *args) {
    10621117  //Limit slopes for each volume to eliminate artificial variance
    10631118  //introduced by e.g. second order extrapolator
     
    10861141        if (!PyArg_ParseTuple(args, "O", &quantity)) {
    10871142          PyErr_SetString(PyExc_RuntimeError,
    1088                           "quantity_ext.c: limit_by_edge could not parse input");
     1143                          "quantity_ext.c: limit_edges_by_all_neighbours could not parse input");
    10891144          return NULL;
    10901145        }
     
    10931148        if (!domain) {
    10941149          PyErr_SetString(PyExc_RuntimeError,
    1095                           "quantity_ext.c: limit_by_edge could not obtain domain object from quantity");                       
     1150                          "quantity_ext.c: limit_edges_by_all_neighbours could not obtain domain object from quantity");                       
    10961151         
    10971152          return NULL;
     
    11021157        if (!Tmp) {
    11031158          PyErr_SetString(PyExc_RuntimeError,
    1104                           "quantity_ext.c: limit_by_edge could not obtain beta_w object from domain");                 
     1159                          "quantity_ext.c: limit_edges_by_all_neighbours could not obtain beta_w object from domain");                 
    11051160         
    11061161          return NULL;
     
    11181173        N = centroid_values -> dimensions[0];
    11191174
    1120         err = _limit_by_edge(N, beta_w,
    1121                                (double*) centroid_values -> data,
    1122                                (double*) vertex_values -> data,
    1123                                (double*) edge_values -> data,
    1124                                (long*)   neighbours -> data);
     1175        err = _limit_edges_by_all_neighbours(N, beta_w,
     1176                                             (double*) centroid_values -> data,
     1177                                             (double*) vertex_values -> data,
     1178                                             (double*) edge_values -> data,
     1179                                             (long*)   neighbours -> data);
    11251180
    11261181        if (err != 0) {
     
    11431198
    11441199
     1200PyObject *limit_edges_by_neighbour(PyObject *self, PyObject *args) {
     1201  //Limit slopes for each volume to eliminate artificial variance
     1202  //introduced by e.g. second order extrapolator
     1203
     1204  //This is an unsophisticated limiter as it does not take into
     1205  //account dependencies among quantities.
     1206
     1207  //precondition:
     1208  //  vertex values are estimated from gradient
     1209  //postcondition:
     1210  //  vertex and edge values are updated
     1211  //
     1212
     1213        PyObject *quantity, *domain, *Tmp;
     1214        PyArrayObject
     1215            *vertex_values,   //Conserved quantities at vertices
     1216            *centroid_values, //Conserved quantities at centroids
     1217            *edge_values,     //Conserved quantities at edges
     1218            *neighbours;
     1219
     1220        double beta_w; //Safety factor
     1221        int N, err;
     1222
     1223
     1224        // Convert Python arguments to C
     1225        if (!PyArg_ParseTuple(args, "O", &quantity)) {
     1226          PyErr_SetString(PyExc_RuntimeError,
     1227                          "quantity_ext.c: limit_edges_by_neighbour could not parse input");
     1228          return NULL;
     1229        }
     1230
     1231        domain = PyObject_GetAttrString(quantity, "domain");
     1232        if (!domain) {
     1233          PyErr_SetString(PyExc_RuntimeError,
     1234                          "quantity_ext.c: limit_edges_by_neighbour could not obtain domain object from quantity");                     
     1235         
     1236          return NULL;
     1237        }
     1238
     1239        // Get safety factor beta_w
     1240        Tmp = PyObject_GetAttrString(domain, "beta_w");
     1241        if (!Tmp) {
     1242          PyErr_SetString(PyExc_RuntimeError,
     1243                          "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain");                       
     1244         
     1245          return NULL;
     1246        }       
     1247
     1248
     1249        // Get pertinent variables
     1250        neighbours       = get_consecutive_array(domain, "neighbours");
     1251        centroid_values  = get_consecutive_array(quantity, "centroid_values");
     1252        vertex_values    = get_consecutive_array(quantity, "vertex_values");
     1253        edge_values      = get_consecutive_array(quantity, "edge_values");
     1254        beta_w           = PyFloat_AsDouble(Tmp);
     1255
     1256
     1257        N = centroid_values -> dimensions[0];
     1258
     1259        err = _limit_edges_by_neighbour(N, beta_w,
     1260                                        (double*) centroid_values -> data,
     1261                                        (double*) vertex_values -> data,
     1262                                        (double*) edge_values -> data,
     1263                                        (long*)   neighbours -> data);
     1264       
     1265        if (err != 0) {
     1266          PyErr_SetString(PyExc_RuntimeError,
     1267                          "Internal function _limit_by_vertex failed");
     1268          return NULL;
     1269        }       
     1270
     1271
     1272        // Release
     1273        Py_DECREF(neighbours);
     1274        Py_DECREF(centroid_values);
     1275        Py_DECREF(vertex_values);
     1276        Py_DECREF(edge_values);
     1277        Py_DECREF(Tmp);
     1278
     1279
     1280        return Py_BuildValue("");
     1281}
     1282
    11451283
    11461284// Method table for python module
    11471285static struct PyMethodDef MethodTable[] = {
    11481286        {"limit_old", limit_old, METH_VARARGS, "Print out"},
    1149         {"limit_by_vertex", limit_by_vertex, METH_VARARGS, "Print out"},
    1150         {"limit_by_edge", limit_by_edge, METH_VARARGS, "Print out"},
     1287        {"limit_vertices_by_all_neighbours", limit_vertices_by_all_neighbours, METH_VARARGS, "Print out"},
     1288        {"limit_edges_by_all_neighbours", limit_edges_by_all_neighbours, METH_VARARGS, "Print out"},
     1289        {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"},
    11511290        {"update", update, METH_VARARGS, "Print out"},
    11521291        {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
  • anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

    r4886 r4897  
    14451445
    14461446
    1447     def test_vertex_limiter(self):
     1447    def test_limit_vertices_by_all_neighbours(self):
    14481448        quantity = Conserved_quantity(self.mesh4)
    14491449
     
    14531453
    14541454        #Limit
    1455         quantity.limit_by_vertex()
     1455        quantity.limit_vertices_by_all_neighbours()
    14561456
    14571457        #Assert that central triangle is limited by neighbours
     
    14751475
    14761476
    1477     def test_edge_limiter(self):
     1477    def test_limit_edges_by_all_neighbours(self):
    14781478        quantity = Conserved_quantity(self.mesh4)
    14791479
     
    14831483
    14841484        #Limit
    1485         quantity.limit_by_edge()
     1485        quantity.limit_edges_by_all_neighbours()
    14861486
    14871487        #Assert that central triangle is limited by neighbours
     
    14931493
    14941494        assert quantity.edge_values[1,2] <= quantity.centroid_values[2]
    1495         assert quantity.edge_values[1,2] <= quantity.centroid_values[0]
     1495        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
    14961496
    14971497
     
    15031503                             sum(quantity.vertex_values[k,:])/3)
    15041504
     1505
     1506    def test_limit_edges_by_neighbour(self):
     1507        quantity = Conserved_quantity(self.mesh4)
     1508
     1509        #Create a deliberate overshoot (e.g. from gradient computation)
     1510        quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
     1511
     1512
     1513        #Limit
     1514        quantity.limit_edges_by_neighbour()
     1515
     1516        #Assert that central triangle is limited by neighbours
     1517        assert quantity.edge_values[1,0] <= quantity.centroid_values[3]
     1518        assert quantity.edge_values[1,0] >= quantity.centroid_values[1]
     1519
     1520        assert quantity.edge_values[1,1] <= quantity.centroid_values[2]
     1521        assert quantity.edge_values[1,1] >= quantity.centroid_values[1]
     1522
     1523        assert quantity.edge_values[1,2] <= quantity.centroid_values[1]
     1524        assert quantity.edge_values[1,2] >= quantity.centroid_values[0]
     1525
     1526
     1527
     1528        #Assert that quantities are conserved
     1529        from Numeric import sum
     1530        for k in range(quantity.centroid_values.shape[0]):
     1531            assert allclose (quantity.centroid_values[k],
     1532                             sum(quantity.vertex_values[k,:])/3)
    15051533
    15061534    def test_limiter2(self):
     
    22072235        assert allclose(subset, answer)
    22082236
     2237    def test_smooth_vertex_values(self):
     2238        """
     2239        get values based on triangle lists.
     2240        """
     2241        from mesh_factory import rectangular
     2242        from shallow_water import Domain
     2243        from Numeric import zeros, Float
     2244
     2245        #Create basic mesh
     2246        points, vertices, boundary = rectangular(2, 2)
     2247
     2248        #print "points",points
     2249        #print "vertices",vertices
     2250        #print "boundary",boundary
     2251
     2252        #Create shallow water domain
     2253        domain = Domain(points, vertices, boundary)
     2254        #print "domain.number_of_elements ",domain.number_of_elements
     2255        quantity = Quantity(domain,[[0,0,0],[1,1,1],[2,2,2],[3,3,3],
     2256                                    [4,4,4],[5,5,5],[6,6,6],[7,7,7]])
     2257
     2258        #print "quantity.get_values(location = 'unique vertices')", \
     2259        #      quantity.get_values(location = 'unique vertices')
     2260
     2261        #print "quantity.get_values(location = 'unique vertices')", \
     2262        #      quantity.get_values(indices=[0,1,2,3,4,5,6,7], \
     2263        #                          location = 'unique vertices')
     2264
     2265        #print quantity.get_values(location = 'unique vertices')
     2266        #print quantity.domain.number_of_triangles_per_node
     2267        #print quantity.vertex_values
     2268       
     2269        #answer = [0.5, 2, 3, 3, 3.5, 4, 4, 5, 6.5]
     2270        #assert allclose(answer,
     2271        #                quantity.get_values(location = 'unique vertices'))
     2272
     2273        quantity.smooth_vertex_values()
     2274
     2275        #print quantity.vertex_values
     2276
     2277
     2278        answer_vertex_values = [[3,3.5,0.5],[2,0.5,3.5],[3.5,4,2],[3,2,4],
     2279                                [4,5,3],[3.5,3,5],[5,6.5,3.5],[4,3.5,6.5]]
     2280       
     2281        assert allclose(answer_vertex_values,
     2282                        quantity.vertex_values)
     2283        #print "quantity.centroid_values",quantity.centroid_values
     2284        #print "quantity.get_values(location = 'centroids') ",\
     2285        #      quantity.get_values(location = 'centroids')
    22092286
    22102287
  • anuga_core/source/anuga/visualiser/realtime.py

    r3966 r4897  
    1 from Numeric import Float, zeros
     1from Numeric import Float, zeros, shape
    22from Tkinter import Button, E, Tk, W
    33from threading import Event
     
    4545        verticies = self.source.get_vertex_coordinates()
    4646        N_vert = len(verticies)
     47
    4748        # Also build vert_index - a list of the x & y values of each vertex
    4849        self.vert_index = zeros((N_vert,2), Float)
     
    5051            self.vtk_cells.InsertNextCell(3)
    5152            for v in range(3):
    52                 self.vert_index[triangles[n][v]] = verticies[n * 3 + v]
    53                 self.vtk_cells.InsertCellPoint(triangles[n][v])
     53                self.vert_index[n * 3 + v] = verticies[n * 3 + v]
     54                self.vtk_cells.InsertCellPoint(n * 3 + v)
    5455
    5556    def update_height_quantity(self, quantityName, dynamic=True):
     
    5960        vertex_values, _ = self.source.get_quantity(quantityName).get_vertex_values(xy=False, smooth=False)
    6061
    61         for n in range(len(triangles)):
    62             for v in range(3):
    63                 #qty_index[triangles[n][v]] = self.source.get_quantity(quantityName).vertex_values[n][v]
    64                 qty_index[triangles[n][v]] = vertex_values[n * 3 + v]
     62        for n in range(N_vert):
     63            qty_index[n] = vertex_values[n]
    6564
    6665        points = vtkPoints()
     
    124123        Visualiser.redraw(self)
    125124
    126     def update(self):
     125    def update(self,pause=False):
    127126        """Sync the visualiser to the domain. Call this in the evolve loop."""
     127           
    128128        if self.running:
    129129            self.sync_redrawReady.set()
     
    132132            self.sync_unpaused.wait()
    133133
     134        if pause and self.running:
     135            if self.sync_unpaused.isSet():
     136                self.sync_unpaused.clear()
     137                self.tk_pauseResume.config(text="Resume")
     138               
     139                self.sync_redrawReady.set()
     140                self.sync_idle.wait()
     141                self.sync_idle.clear()
     142                self.sync_unpaused.wait()
     143           
     144
     145
    134146    def evolveFinished(self):
    135147        """Stop the visualiser from waiting on signals from the evolve loop.
Note: See TracChangeset for help on using the changeset viewer.