# Changeset 4886

Ignore:
Timestamp:
Dec 12, 2007, 5:38:37 PM (16 years ago)
Message:

Working towards an edge limiter (at least in quantity)

Location:
anuga_core/source/anuga
Files:
4 edited

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

 r4883 #(either from this module or C-extension) interpolate_from_vertices_to_edges(self) def interpolate_from_edges_to_vertices(self): #Call correct module function #(either from this module or C-extension) interpolate_from_edges_to_vertices(self) ##     def limit_by_vertex(self): ##         #Call correct module function ##         #(either from this module or C-extension) ##         limit_by_vertex() ##     def limit_by_edge(self): ##         #Call correct module function ##         #(either from this module or C-extension) ##         limit_by_edge() def limit_by_vertex(self): #Call correct module function #(either from this module or C-extension) limit_by_vertex(self) def limit_by_edge(self): #Call correct module function #(either from this module or C-extension) limit_by_edge(self) def extrapolate_second_order(self): compute_gradients,\ limit_old,\ limit_by_vertex,\ limit_by_edge,\ extrapolate_second_order,\ interpolate_from_vertices_to_edges,\ interpolate_from_edges_to_vertices,\ update else:
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c

 r4883 y2 = vertex_coordinates[k6 + 5]; // Extrapolate // Extrapolate to Vertices vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y); vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y); vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y); // Extrapolate to Edges (midpoints) edge_values[k3+0] = 0.5*(vertex_values[k3 + 1]+vertex_values[k3 + 2]); edge_values[k3+1] = 0.5*(vertex_values[k3 + 2]+vertex_values[k3 + 0]); int _interpolate(int N, int _limit_by_vertex(int N, double beta, double* centroid_values, double* vertex_values, double* edge_values, long*   neighbours) { int i, k, k2, k3, k6; long n; double qmin, qmax, qn, qc; double dq, dqa[3], phi, r; for (k=0; k= 0) { qn = centroid_values[n]; //Neighbour's centroid value qmin = min(qmin, qn); qmax = max(qmax, qn); } } phi = 1.0; for (i=0; i<3; i++) { r = 1.0; dq = vertex_values[k3+i] - qc;    //Delta between vertex and centroid values dqa[i] = dq;                      //Save dq for use in updating vertex values if (dq > 0.0) r = (qmax - qc)/dq; if (dq < 0.0) r = (qmin - qc)/dq; phi = min( min(r*beta, 1.0), phi); } //Update vertex and edge values using phi limiter vertex_values[k3+0] = qc + phi*dqa[0]; vertex_values[k3+1] = qc + phi*dqa[1]; vertex_values[k3+2] = qc + phi*dqa[2]; edge_values[k3+0] = 0.5*(vertex_values[k3+1] + vertex_values[k3+2]); edge_values[k3+1] = 0.5*(vertex_values[k3+2] + vertex_values[k3+0]); edge_values[k3+2] = 0.5*(vertex_values[k3+0] + vertex_values[k3+1]); } return 0; } int _limit_by_edge(int N, double beta, double* centroid_values, double* vertex_values, double* edge_values, long*   neighbours) { int i, k, k2, k3, k6; long n; double qmin, qmax, qn, qc; double dq, dqa[3], phi, r; for (k=0; k= 0) { qn = centroid_values[n]; //Neighbour's centroid value qmin = min(qmin, qn); qmax = max(qmax, qn); } } phi = 1.0; for (i=0; i<3; i++) { r = 1.0; dq = edge_values[k3+i] - qc;     //Delta between edge and centroid values dqa[i] = dq;                      //Save dq for use in updating vertex values if (dq > 0.0) r = (qmax - qc)/dq; if (dq < 0.0) r = (qmin - qc)/dq; phi = min( min(r*beta, 1.0), phi); } //Update vertex and edge values using phi limiter edge_values[k3+0] = qc + phi*dqa[0]; edge_values[k3+1] = qc + phi*dqa[1]; edge_values[k3+2] = qc + phi*dqa[2]; vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2] - edge_values[k3+0]; vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0] - edge_values[k3+1]; vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1] - edge_values[k3+2]; } return 0; } int _interpolate_from_vertices_to_edges(int N, double* vertex_values, double* edge_values) { edge_values[k3 + 1] = 0.5*(q0+q2); edge_values[k3 + 2] = 0.5*(q0+q1); } return 0; } int _interpolate_from_edges_to_vertices(int N, double* vertex_values, double* edge_values) { int k, k3; double e0, e1, e2; for (k=0; k dimensions[0]; err = _interpolate(N, err = _interpolate_from_vertices_to_edges(N, (double*) vertex_values -> data, (double*) edge_values -> data); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed"); return NULL; } // Release and return Py_DECREF(vertex_values); Py_DECREF(edge_values); return Py_BuildValue(""); } PyObject *interpolate_from_edges_to_vertices(PyObject *self, PyObject *args) { // //Compute vertex values from edge values using linear interpolation // PyObject *quantity; PyArrayObject *vertex_values, *edge_values; int N, err; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "O", &quantity)) { PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: interpolate_from_edges_to_vertices could not parse input"); return NULL; } vertex_values = get_consecutive_array(quantity, "vertex_values"); edge_values = get_consecutive_array(quantity, "edge_values"); N = vertex_values -> dimensions[0]; err = _interpolate_from_edges_to_vertices(N, (double*) vertex_values -> data, (double*) edge_values -> data); PyObject *limit_by_vertex(PyObject *self, PyObject *args) { //Limit slopes for each volume to eliminate artificial variance //introduced by e.g. second order extrapolator //This is an unsophisticated limiter as it does not take into //account dependencies among quantities. //precondition: //  vertex values are estimated from gradient //postcondition: //  vertex and edge values are updated // PyObject *quantity, *domain, *Tmp; PyArrayObject *vertex_values,   //Conserved quantities at vertices *centroid_values, //Conserved quantities at centroids *edge_values,     //Conserved quantities at edges *neighbours; double beta_w; //Safety factor int N, err; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "O", &quantity)) { PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_by_vertex could not parse input"); return NULL; } domain = PyObject_GetAttrString(quantity, "domain"); if (!domain) { PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_by_vertex could not obtain domain object from quantity"); return NULL; } // Get safety factor beta_w Tmp = PyObject_GetAttrString(domain, "beta_w"); if (!Tmp) { PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain"); return NULL; } // Get pertinent variables neighbours       = get_consecutive_array(domain, "neighbours"); centroid_values  = get_consecutive_array(quantity, "centroid_values"); vertex_values    = get_consecutive_array(quantity, "vertex_values"); edge_values      = get_consecutive_array(quantity, "edge_values"); beta_w           = PyFloat_AsDouble(Tmp); N = centroid_values -> dimensions[0]; err = _limit_by_vertex(N, beta_w, (double*) centroid_values -> data, (double*) vertex_values -> data, (double*) edge_values -> data, (long*)   neighbours -> data); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "Internal function _limit_by_vertex failed"); return NULL; } // Release Py_DECREF(neighbours); Py_DECREF(centroid_values); Py_DECREF(vertex_values); Py_DECREF(edge_values); Py_DECREF(Tmp); return Py_BuildValue(""); } PyObject *limit_by_edge(PyObject *self, PyObject *args) { //Limit slopes for each volume to eliminate artificial variance //introduced by e.g. second order extrapolator //This is an unsophisticated limiter as it does not take into //account dependencies among quantities. //precondition: //  vertex values are estimated from gradient //postcondition: //  vertex and edge values are updated // PyObject *quantity, *domain, *Tmp; PyArrayObject *vertex_values,   //Conserved quantities at vertices *centroid_values, //Conserved quantities at centroids *edge_values,     //Conserved quantities at edges *neighbours; double beta_w; //Safety factor int N, err; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "O", &quantity)) { PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_by_edge could not parse input"); return NULL; } domain = PyObject_GetAttrString(quantity, "domain"); if (!domain) { PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_by_edge could not obtain domain object from quantity"); return NULL; } // Get safety factor beta_w Tmp = PyObject_GetAttrString(domain, "beta_w"); if (!Tmp) { PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_by_edge could not obtain beta_w object from domain"); return NULL; } // Get pertinent variables neighbours       = get_consecutive_array(domain, "neighbours"); centroid_values  = get_consecutive_array(quantity, "centroid_values"); vertex_values    = get_consecutive_array(quantity, "vertex_values"); edge_values      = get_consecutive_array(quantity, "edge_values"); beta_w           = PyFloat_AsDouble(Tmp); N = centroid_values -> dimensions[0]; err = _limit_by_edge(N, beta_w, (double*) centroid_values -> data, (double*) vertex_values -> data, (double*) edge_values -> data, (long*)   neighbours -> data); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "Internal function _limit_by_vertex failed"); return NULL; } // Release Py_DECREF(neighbours); Py_DECREF(centroid_values); Py_DECREF(vertex_values); Py_DECREF(edge_values); Py_DECREF(Tmp); return Py_BuildValue(""); } // Method table for python module static struct PyMethodDef MethodTable[] = { {"limit_old", limit_old, METH_VARARGS, "Print out"}, {"limit_by_vertex", limit_by_vertex, METH_VARARGS, "Print out"}, {"limit_by_edge", limit_by_edge, METH_VARARGS, "Print out"}, {"update", update, METH_VARARGS, "Print out"}, {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"}, interpolate_from_vertices_to_edges, METH_VARARGS, "Print out"}, {"interpolate_from_edges_to_vertices", interpolate_from_edges_to_vertices, METH_VARARGS, "Print out"}, {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"}, {NULL, NULL, 0, NULL}   // sentinel
• ## anuga_core/source/anuga/abstract_2d_finite_volumes/test_quantity.py

 r4883 def test_limiter(self): def test_vertex_limiter(self): quantity = Conserved_quantity(self.mesh4) #Limit quantity.limit() quantity.limit_by_vertex() #Assert that central triangle is limited by neighbours def test_edge_limiter(self): quantity = Conserved_quantity(self.mesh4) #Create a deliberate overshoot (e.g. from gradient computation) quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]]) #Limit quantity.limit_by_edge() #Assert that central triangle is limited by neighbours assert quantity.edge_values[1,0] <= quantity.centroid_values[2] assert quantity.edge_values[1,0] >= quantity.centroid_values[0] assert quantity.edge_values[1,1] <= quantity.centroid_values[2] assert quantity.edge_values[1,1] >= quantity.centroid_values[0] assert quantity.edge_values[1,2] <= quantity.centroid_values[2] assert quantity.edge_values[1,2] <= quantity.centroid_values[0] #Assert that quantities are conserved from Numeric import sum for k in range(quantity.centroid_values.shape[0]): assert allclose (quantity.centroid_values[k], sum(quantity.vertex_values[k,:])/3) def test_limiter2(self): """Taken from test_shallow_water #Extrapolate #Extrapolate from centroid to vertices and edges quantity.extrapolate_first_order() #Interpolate quantity.interpolate_from_vertices_to_edges() #quantity.interpolate_from_vertices_to_edges() assert allclose(quantity.vertex_values, assert allclose(quantity.edge_values, [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) def test_interpolate_from_vertices_to_edges(self): quantity = Quantity(self.mesh4) quantity.vertex_values = array([[1,0,2], [1,2,4], [4,2,5], [3,1,4]],Float) quantity.interpolate_from_vertices_to_edges() assert allclose(quantity.edge_values, [[1., 1.5, 0.5], [3., 2.5, 1.5], [3.5, 4.5, 3.], [2.5, 3.5, 2]]) def test_interpolate_from_edges_to_vertices(self): quantity = Quantity(self.mesh4) quantity.edge_values = array([[1., 1.5, 0.5], [3., 2.5, 1.5], [3.5, 4.5, 3.], [2.5, 3.5, 2]],Float) quantity.interpolate_from_edges_to_vertices() assert allclose(quantity.vertex_values, [[1,0,2], [1,2,4], [4,2,5], [3,1,4]])
• ## anuga_core/source/anuga/utilities/util_ext.h

 r4883 void _limit_by_edge(int N, double beta, double* qc, double* qv, double* qmin, double* qmax) { //N are the number of elements int k, i, k3; double dq, dqa[3], phi, r; //printf("INSIDE\n"); for (k=0; k 0.0) r = (qmax[k] - qc[k])/dq; if (dq < 0.0) r = (qmin[k] - qc[k])/dq; phi = min( min(r*beta, 1.0), phi); } //Then update using phi limiter for (i=0; i<3; i++) { qv[k3+i] = qc[k] + phi*dqa[i]; } } } void _limit_by_vertex(int N, double beta, double* qc, double* qv, double* qmin, double* qmax) { //N are the number of elements int k, i, k3; double dq, dqa[3], phi, r; //printf("INSIDE\n"); for (k=0; k 0.0) r = (qmax[k] - qc[k])/dq; if (dq < 0.0) r = (qmin[k] - qc[k])/dq; phi = min( min(r*beta, 1.0), phi); } //Then update using phi limiter for (i=0; i<3; i++) { qv[k3+i] = qc[k] + phi*dqa[i]; } } }
Note: See TracChangeset for help on using the changeset viewer.