# Changeset 587

Ignore:
Timestamp:
Nov 17, 2004, 11:40:28 PM (19 years ago)
Message:

testing sparse

Location:
inundation/ga/storm_surge
Files:
5 edited

Unmodified
Removed
• ## inundation/ga/storm_surge/analytical solutions/oblique_test_chris.py

 r564 leny = 30. lenx = 40. n = 50 m = 60 points, elements, boundary = oblique(m, n, lenx, leny) # Order of solver domain.default_order=1 domain.default_order=2 # Store output D = Dirichlet_boundary([1.0, 8.57, 0.0]) domain.set_boundary({'left': D, 'right': T, 'top': T, 'bottom': R}) domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R}) ######################
• ## inundation/ga/storm_surge/pyvolution/quantity_ext.c

 r515 int N, err; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "Od", &quantity, ×tep))
• ## inundation/ga/storm_surge/pyvolution/sparse.py

 r586 print 'FIXME: Only Numeric types implemented so far' ##        R = csr_mv(self,B) #Assume numeric types from now on return R def csr_mv_python(self, B): from Numeric import zeros, Float if len(B.shape) == 1: #Vector assert B.shape[0] == self.N, 'Mismatching dimensions' R = zeros(self.M, Float) #Result #Multiply nonzero elements for i in range(self.M): for ckey in range(self.row_ptr[i],self.row_ptr[i+1]): j = self.colind[ckey] R[i] += self.data[ckey]*B[j] return R else: raise ValueError, 'Dimension too high: d=%d' %len(B.shape) return R import compile if compile.can_use_C_extension('sparse_ext.c'): #Replace python version with c implementations from sparse_ext import csr_mv if __name__ == '__main__':
• ## inundation/ga/storm_surge/pyvolution/sparse_ext.c

 r586 #include "math.h" //Shared code snippets #include "util_ext.h" #include "quantity_ext.h" int _compute_csr_mult(int M, double* data, int* colind, int* row_ptr, double* R) { int _csr_mv(int M, double* data, int* colind, int* row_ptr, double* x, double* y) { int i, j, ckey; for (i=0; i dimensions[0])-1; //Allocate space for return vectors y (don't DECREF) dimensions[0] = M; y = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); err = _csr_mv(M, (double*) data -> data, (int*)    colind -> data, (int*)    row_ptr -> data, (double*) x -> data, (double*) y -> data); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "matrix vector mult could not be calculated"); return NULL; } return 0; //Release Py_DECREF(data); Py_DECREF(colind); Py_DECREF(row_ptr); Py_DECREF(x); //Build result, release and return R = Py_BuildValue("O", PyArray_Return(y)); Py_DECREF(y); return R; } int _interpolate(int N, double* vertex_values, double* edge_values) { int k, k3; double q0, q1, q2; for (k=0; k dimensions[0]; err = _update(N, timestep, (double*) centroid_values -> data, (double*) explicit_update -> data, (double*) semi_implicit_update -> data); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "Zero division in semi implicit update - call Stephen :)"); return NULL; } //Release and return Py_DECREF(centroid_values); Py_DECREF(explicit_update); Py_DECREF(semi_implicit_update); return Py_BuildValue(""); } PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) { PyObject *quantity; PyArrayObject *vertex_values, *edge_values; int N, err; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "O", &quantity)) return NULL; vertex_values = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values"); if (!vertex_values) return NULL; edge_values = (PyArrayObject*) PyObject_GetAttrString(quantity, "edge_values"); if (!edge_values) return NULL; N = vertex_values -> dimensions[0]; err = _interpolate(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 *compute_gradients(PyObject *self, PyObject *args) { PyObject *quantity, *domain, *R; PyArrayObject *centroids,            //Coordinates at centroids *centroid_values,      //Values at centroids *number_of_boundaries, //Number of boundaries for each triangle *surrogate_neighbours, //True neighbours or - if one missing - self *a, *b;                //Return values int dimensions[1], N, err; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "O", &quantity)) return NULL; domain = PyObject_GetAttrString(quantity, "domain"); if (!domain) return NULL; //Get pertinent variables centroids = (PyArrayObject*) PyObject_GetAttrString(domain, "centroid_coordinates"); if (!centroids) return NULL; centroid_values = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); if (!centroid_values) return NULL; surrogate_neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "surrogate_neighbours"); if (!surrogate_neighbours) return NULL; number_of_boundaries = (PyArrayObject*) PyObject_GetAttrString(domain, "number_of_boundaries"); if (!number_of_boundaries) return NULL; N = centroid_values -> dimensions[0]; //Release Py_DECREF(domain); //Allocate space for return vectors a and b (don't DECREF) dimensions[0] = N; a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); err = _compute_gradients(N, (double*) centroids -> data, (double*) centroid_values -> data, (int*) number_of_boundaries -> data, (int*) surrogate_neighbours -> data, (double*) a -> data, (double*) b -> data); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); return NULL; } //Release Py_DECREF(centroids); Py_DECREF(centroid_values); Py_DECREF(number_of_boundaries); Py_DECREF(surrogate_neighbours); //Build result, release and return R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b)); Py_DECREF(a); Py_DECREF(b); return R; } PyObject *extrapolate_second_order(PyObject *self, PyObject *args) { PyObject *quantity, *domain; PyArrayObject *centroids,            //Coordinates at centroids *centroid_values,      //Values at centroids *vertex_coordinates,   //Coordinates at vertices *vertex_values,        //Values at vertices *number_of_boundaries, //Number of boundaries for each triangle *surrogate_neighbours, //True neighbours or - if one missing - self *a, *b;                //Gradients //int N, err; int dimensions[1], N, err; //double *a, *b;  //Gradients // Convert Python arguments to C if (!PyArg_ParseTuple(args, "O", &quantity)) return NULL; domain = PyObject_GetAttrString(quantity, "domain"); if (!domain) return NULL; //Get pertinent variables centroids = (PyArrayObject*) PyObject_GetAttrString(domain, "centroid_coordinates"); if (!centroids) return NULL; centroid_values = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); if (!centroid_values) return NULL; surrogate_neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "surrogate_neighbours"); if (!surrogate_neighbours) return NULL; number_of_boundaries = (PyArrayObject*) PyObject_GetAttrString(domain, "number_of_boundaries"); if (!number_of_boundaries) return NULL; vertex_coordinates = (PyArrayObject*) PyObject_GetAttrString(domain, "vertex_coordinates"); if (!vertex_coordinates) return NULL; vertex_values = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values"); if (!vertex_values) return NULL; /* printf("In extrapolate C routine\n"); printf("d0=%d, d1=%d\n", vertex_values -> dimensions[0], vertex_values -> dimensions[1]); */ vertex_values = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values"); if (!vertex_values) return NULL; N = centroid_values -> dimensions[0]; //Release Py_DECREF(domain); //Allocate space for return vectors a and b (don't DECREF) dimensions[0] = N; a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE); //FIXME: Odd that I couldn't use normal arrays //Allocate space for return vectors a and b (don't DECREF) //a = (double*) malloc(N * sizeof(double)); //if (!a) return NULL; //b = (double*) malloc(N * sizeof(double)); //if (!b) return NULL; err = _compute_gradients(N, (double*) centroids -> data, (double*) centroid_values -> data, (int*) number_of_boundaries -> data, (int*) surrogate_neighbours -> data, (double*) a -> data, (double*) b -> data); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed"); return NULL; } err = _extrapolate(N, (double*) centroids -> data, (double*) centroid_values -> data, (double*) vertex_coordinates -> data, (double*) vertex_values -> data, (double*) a -> data, (double*) b -> data); //a, b); if (err != 0) { PyErr_SetString(PyExc_RuntimeError, "Internal function _extrapolate failed"); return NULL; } //Release Py_DECREF(centroids); Py_DECREF(centroid_values); Py_DECREF(number_of_boundaries); Py_DECREF(surrogate_neighbours); Py_DECREF(vertex_coordinates); Py_DECREF(vertex_values); Py_DECREF(a); Py_DECREF(b); //free(a); //free(b); return Py_BuildValue(""); } PyObject *limit(PyObject *self, PyObject *args) { PyObject *quantity, *domain, *Tmp; PyArrayObject *qv, //Conserved quantities at vertices *qc, //Conserved quantities at centroids *neighbours; int k, i, n, N, k3; double beta; //Safety factor double *qmin, *qmax, qn; // Convert Python arguments to C if (!PyArg_ParseTuple(args, "O", &quantity)) return NULL; domain = PyObject_GetAttrString(quantity, "domain"); if (!domain) return NULL; neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours"); //Get safety factor beta Tmp = PyObject_GetAttrString(domain, "beta"); if (!Tmp) return NULL; beta = PyFloat_AsDouble(Tmp); Py_DECREF(Tmp); Py_DECREF(domain); qc = (PyArrayObject*) PyObject_GetAttrString(quantity, "centroid_values"); qv = (PyArrayObject*) PyObject_GetAttrString(quantity, "vertex_values"); N = qc -> dimensions[0]; //Find min and max of this and neighbour's centroid values qmin = malloc(N * sizeof(double)); qmax = malloc(N * sizeof(double)); for (k=0; k data)[k]; qmax[k] = qmin[k]; for (i=0; i<3; i++) { n = ((int*) neighbours -> data)[k3+i]; if (n >= 0) { qn = ((double*) qc -> data)[n]; //Neighbour's centroid value qmin[k] = min(qmin[k], qn); qmax[k] = max(qmax[k], qn); } } } // Call underlying routine _limit(N, beta, (double*) qc -> data, (double*) qv -> data, qmin, qmax); free(qmin); free(qmax); return Py_BuildValue(""); } // Method table for python module static struct PyMethodDef MethodTable[] = { {"limit", limit, METH_VARARGS, "Print out"}, {"update", update, METH_VARARGS, "Print out"}, {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"}, {"extrapolate_second_order", extrapolate_second_order, METH_VARARGS, "Print out"}, {"interpolate_from_vertices_to_edges", interpolate_from_vertices_to_edges, METH_VARARGS, "Print out"}, {"csr_mv", csr_mv, METH_VARARGS, "Print out"}, {NULL, NULL, 0, NULL}   /* sentinel */ }; // Module initialisation void initquantity_ext(void){ Py_InitModule("quantity_ext", MethodTable); Py_InitModule("sparse_ext", MethodTable); import_array();     //Necessary for handling of NumPY structures
• ## inundation/ga/storm_surge/pyvolution/wind_example.py

 r522 #Create basic mesh (100m x 100m) N = 80 N = 100 len = 100 points, vertices, boundary = rectangular(N, N, len, len) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.default_order = 2 domain.default_order = 2 domain.store = True domain.set_name('wind_laminar') domain.set_quantity('elevation', 0.0) domain.set_quantity('level', 1.0) domain.set_quantity('friction', 0.0) domain.set_quantity('friction', 0.1) #Constant (quite extreme :-) windfield: 9000 m/s, bearing 120 degrees domain.forcing_terms.append( Wind_stress(s=9000, phi=120) ) domain.forcing_terms.append( Wind_stress(s=9000, phi=135) ) ######################
Note: See TracChangeset for help on using the changeset viewer.