source: inundation/pyvolution/quantity_ext.c @ 1884

Last change on this file since 1884 was 1506, checked in by matthew, 19 years ago

The re-initialization of semi_implicit_update[k] to 0.0 is now done in the _update C function of quantity_ext.c, instead of the method update_conserved_quantities of domain.py. For run_profile.py, with N=128, the time spend in update_conserved_quantities is cut from 14.00 secs to 8.31 secs (averaged over three runs).

File size: 13.7 KB
Line 
1// Python - C extension for quantity module.
2//
3// To compile (Python2.3):
4//  gcc -c util_ext.c -I/usr/include/python2.3 -o util_ext.o -Wall -O
5//  gcc -shared util_ext.o  -o util_ext.so
6//
7// See the module quantity.py
8//
9//
10// Ole Nielsen, GA 2004
11
12#include "Python.h"
13#include "Numeric/arrayobject.h"
14#include "math.h"
15
16//Shared code snippets
17#include "util_ext.h"
18//#include "quantity_ext.h" //Obsolete
19
20
21
22int _compute_gradients(int N,
23                        double* centroids,
24                        double* centroid_values,
25                        long* number_of_boundaries,
26                        long* surrogate_neighbours,
27                        double* a,
28                        double* b) {
29
30  int i, k, k0, k1, k2, index3;
31  double x0, x1, x2, y0, y1, y2, q0, q1, q2; //, det;
32
33
34  for (k=0; k<N; k++) {
35    index3 = 3*k;
36
37    if (number_of_boundaries[k] < 2) {
38      //Two or three true neighbours
39
40      //Get indices of neighbours (or self when used as surrogate)
41      //k0, k1, k2 = surrogate_neighbours[k,:]
42
43      k0 = surrogate_neighbours[index3 + 0];
44      k1 = surrogate_neighbours[index3 + 1];
45      k2 = surrogate_neighbours[index3 + 2];
46
47
48      if (k0 == k1 || k1 == k2) return -1;
49
50      //Get data
51      q0 = centroid_values[k0];
52      q1 = centroid_values[k1];
53      q2 = centroid_values[k2];
54
55      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
56      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
57      x2 = centroids[k2*2]; y2 = centroids[k2*2+1];
58
59      //Gradient
60      _gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2, &a[k], &b[k]);
61
62    } else if (number_of_boundaries[k] == 2) {
63      //One true neighbour
64
65      //#Get index of the one neighbour
66      i=0; k0 = k;
67      while (i<3 && k0==k) {
68        k0 = surrogate_neighbours[index3 + i];
69        i++;
70      }
71      if (k0 == k) return -1;
72
73      k1 = k; //self
74
75      //Get data
76      q0 = centroid_values[k0];
77      q1 = centroid_values[k1];
78
79      x0 = centroids[k0*2]; y0 = centroids[k0*2+1];
80      x1 = centroids[k1*2]; y1 = centroids[k1*2+1];
81
82      //Two point gradient
83      _gradient2(x0, y0, x1, y1, q0, q1, &a[k], &b[k]);     
84     
85     
86      //Old (wrong code)
87      //det = x0*y1 - x1*y0;
88      //if (det != 0.0) {
89      //        a[k] = (y1*q0 - y0*q1)/det;
90      //        b[k] = (x0*q1 - x1*q0)/det;
91      //}
92    }
93    //    else:
94    //        #No true neighbours -
95    //        #Fall back to first order scheme
96    //        pass
97  }
98  return 0;
99}
100
101
102int _extrapolate(int N,
103                 double* centroids,
104                 double* centroid_values,
105                 double* vertex_coordinates,
106                 double* vertex_values,
107                 double* a,
108                 double* b) {
109
110  int k, k2, k3, k6;
111  double x, y, x0, y0, x1, y1, x2, y2;
112
113  for (k=0; k<N; k++) {
114    k6 = 6*k;
115    k3 = 3*k;
116    k2 = 2*k;
117
118    //Centroid coordinates
119    x = centroids[k2]; y = centroids[k2+1];
120
121    //vertex coordinates
122    //x0, y0, x1, y1, x2, y2 = X[k,:]
123    x0 = vertex_coordinates[k6 + 0];
124    y0 = vertex_coordinates[k6 + 1];
125    x1 = vertex_coordinates[k6 + 2];
126    y1 = vertex_coordinates[k6 + 3];
127    x2 = vertex_coordinates[k6 + 4];
128    y2 = vertex_coordinates[k6 + 5];
129
130    //Extrapolate
131    vertex_values[k3+0] = centroid_values[k] + a[k]*(x0-x) + b[k]*(y0-y);
132    vertex_values[k3+1] = centroid_values[k] + a[k]*(x1-x) + b[k]*(y1-y);
133    vertex_values[k3+2] = centroid_values[k] + a[k]*(x2-x) + b[k]*(y2-y);
134
135  }
136  return 0;
137}
138
139
140
141
142int _interpolate(int N,
143                 double* vertex_values,
144                 double* edge_values) {
145
146  int k, k3;
147  double q0, q1, q2;
148
149
150  for (k=0; k<N; k++) {
151    k3 = 3*k;
152
153    q0 = vertex_values[k3 + 0];
154    q1 = vertex_values[k3 + 1];
155    q2 = vertex_values[k3 + 2];
156
157    //printf("%f, %f, %f\n", q0, q1, q2);
158    edge_values[k3 + 0] = 0.5*(q1+q2);
159    edge_values[k3 + 1] = 0.5*(q0+q2);
160    edge_values[k3 + 2] = 0.5*(q0+q1);
161  }
162  return 0;
163}
164
165int _update(int N,
166            double timestep,
167            double* centroid_values,
168            double* explicit_update,
169            double* semi_implicit_update) {
170  //Update centroid values based on values stored in
171  //explicit_update and semi_implicit_update as well as given timestep
172
173
174  int k;
175  double denominator, x;
176
177
178  //Divide semi_implicit update by conserved quantity
179  for (k=0; k<N; k++) {
180    x = centroid_values[k];
181    if (x == 0.0) {
182      semi_implicit_update[k] = 0.0;
183    } else {
184      semi_implicit_update[k] /= x;
185    }
186  }
187
188
189  //Explicit updates
190  for (k=0; k<N; k++) {
191    centroid_values[k] += timestep*explicit_update[k];
192  }
193
194  //Semi implicit updates
195  for (k=0; k<N; k++) {
196    denominator = 1.0 - timestep*semi_implicit_update[k];
197    if (denominator == 0.0) {
198      return -1;
199    } else {
200      //Update conserved_quantities from semi implicit updates
201      centroid_values[k] /= denominator;
202    }
203  }
204  //MH080605 set semi_impliit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
205  for (k=0;k<N;k++)
206    semi_implicit_update[k]=0.0;
207  return 0;
208}
209
210
211/////////////////////////////////////////////////
212// Gateways to Python
213PyObject *update(PyObject *self, PyObject *args) {
214
215  PyObject *quantity;
216  PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
217
218  double timestep;
219  int N, err;
220
221
222  // Convert Python arguments to C
223  if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep))
224    return NULL;
225
226  centroid_values = get_consecutive_array(quantity, "centroid_values");
227  explicit_update = get_consecutive_array(quantity, "explicit_update");
228  semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
229
230  N = centroid_values -> dimensions[0];
231
232  err = _update(N, timestep,
233                (double*) centroid_values -> data,
234                (double*) explicit_update -> data,
235                (double*) semi_implicit_update -> data);
236
237
238  if (err != 0) {
239    PyErr_SetString(PyExc_RuntimeError,
240                    "Zero division in semi implicit update - call Stephen :)");
241    return NULL;
242  }
243
244  //Release and return
245  Py_DECREF(centroid_values);
246  Py_DECREF(explicit_update);
247  Py_DECREF(semi_implicit_update);
248
249  return Py_BuildValue("");
250}
251
252
253PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
254
255  PyObject *quantity;
256  PyArrayObject *vertex_values, *edge_values;
257
258  int N, err;
259
260  // Convert Python arguments to C
261  if (!PyArg_ParseTuple(args, "O", &quantity))
262    return NULL;
263
264  vertex_values = get_consecutive_array(quantity, "vertex_values");
265  edge_values = get_consecutive_array(quantity, "edge_values");
266
267  N = vertex_values -> dimensions[0];
268
269  err = _interpolate(N,
270                     (double*) vertex_values -> data,
271                     (double*) edge_values -> data);
272
273  if (err != 0) {
274    PyErr_SetString(PyExc_RuntimeError, "Interpolate could not be computed");
275    return NULL;
276  }
277
278  //Release and return
279  Py_DECREF(vertex_values);
280  Py_DECREF(edge_values);
281
282  return Py_BuildValue("");
283}
284
285
286PyObject *compute_gradients(PyObject *self, PyObject *args) {
287
288  PyObject *quantity, *domain, *R;
289  PyArrayObject
290    *centroids,            //Coordinates at centroids
291    *centroid_values,      //Values at centroids
292    *number_of_boundaries, //Number of boundaries for each triangle
293    *surrogate_neighbours, //True neighbours or - if one missing - self
294    *a, *b;                //Return values
295
296  int dimensions[1], N, err;
297
298  // Convert Python arguments to C
299  if (!PyArg_ParseTuple(args, "O", &quantity))
300    return NULL;
301
302  domain = PyObject_GetAttrString(quantity, "domain");
303  if (!domain)
304    return NULL;
305
306  //Get pertinent variables
307
308  centroids = get_consecutive_array(domain, "centroid_coordinates");
309  centroid_values = get_consecutive_array(quantity, "centroid_values");
310  surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
311  number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
312
313  N = centroid_values -> dimensions[0];
314
315  //Release
316  Py_DECREF(domain);
317
318  //Allocate space for return vectors a and b (don't DECREF)
319  dimensions[0] = N;
320  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
321  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
322
323
324
325  err = _compute_gradients(N,
326                           (double*) centroids -> data,
327                           (double*) centroid_values -> data,
328                           (long*) number_of_boundaries -> data,
329                           (long*) surrogate_neighbours -> data,
330                           (double*) a -> data,
331                           (double*) b -> data);
332
333  if (err != 0) {
334    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
335    return NULL;
336  }
337
338  //Release
339  Py_DECREF(centroids);
340  Py_DECREF(centroid_values);
341  Py_DECREF(number_of_boundaries);
342  Py_DECREF(surrogate_neighbours);
343
344  //Build result, release and return
345  R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
346  Py_DECREF(a);
347  Py_DECREF(b);
348  return R;
349}
350
351
352
353PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
354
355  PyObject *quantity, *domain;
356  PyArrayObject
357    *centroids,            //Coordinates at centroids
358    *centroid_values,      //Values at centroids
359    *vertex_coordinates,   //Coordinates at vertices
360    *vertex_values,        //Values at vertices
361    *number_of_boundaries, //Number of boundaries for each triangle
362    *surrogate_neighbours, //True neighbours or - if one missing - self
363    *a, *b;                //Gradients
364
365  //int N, err;
366  int dimensions[1], N, err;
367  //double *a, *b;  //Gradients
368
369  // Convert Python arguments to C
370  if (!PyArg_ParseTuple(args, "O", &quantity))
371    return NULL;
372
373  domain = PyObject_GetAttrString(quantity, "domain");
374  if (!domain)
375    return NULL;
376
377  //Get pertinent variables
378  centroids = get_consecutive_array(domain, "centroid_coordinates");
379  centroid_values = get_consecutive_array(quantity, "centroid_values");
380  surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
381  number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
382  vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
383  vertex_values = get_consecutive_array(quantity, "vertex_values");
384
385  N = centroid_values -> dimensions[0];
386
387  //Release
388  Py_DECREF(domain);
389
390  //Allocate space for return vectors a and b (don't DECREF)
391  dimensions[0] = N;
392  a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
393  b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
394
395  //FIXME: Odd that I couldn't use normal arrays
396  //Allocate space for return vectors a and b (don't DECREF)
397  //a = (double*) malloc(N * sizeof(double));
398  //if (!a) return NULL;
399  //b = (double*) malloc(N * sizeof(double));
400  //if (!b) return NULL;
401
402
403  err = _compute_gradients(N,
404                           (double*) centroids -> data,
405                           (double*) centroid_values -> data,
406                           (long*) number_of_boundaries -> data,
407                           (long*) surrogate_neighbours -> data,
408                           (double*) a -> data,
409                           (double*) b -> data);
410
411  if (err != 0) {
412    PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
413    return NULL;
414  }
415
416  err = _extrapolate(N,
417                     (double*) centroids -> data,
418                     (double*) centroid_values -> data,
419                     (double*) vertex_coordinates -> data,
420                     (double*) vertex_values -> data,
421                     (double*) a -> data,
422                     (double*) b -> data);
423
424
425  if (err != 0) {
426    PyErr_SetString(PyExc_RuntimeError,
427                    "Internal function _extrapolate failed");
428    return NULL;
429  }
430
431
432
433  //Release
434  Py_DECREF(centroids);
435  Py_DECREF(centroid_values);
436  Py_DECREF(number_of_boundaries);
437  Py_DECREF(surrogate_neighbours);
438  Py_DECREF(vertex_coordinates);
439  Py_DECREF(vertex_values);
440  Py_DECREF(a);
441  Py_DECREF(b);
442
443  return Py_BuildValue("");
444}
445
446
447
448PyObject *limit(PyObject *self, PyObject *args) {
449
450  PyObject *quantity, *domain, *Tmp;
451  PyArrayObject
452    *qv, //Conserved quantities at vertices
453    *qc, //Conserved quantities at centroids
454    *neighbours;
455
456  int k, i, n, N, k3;
457  double beta_w; //Safety factor
458  double *qmin, *qmax, qn;
459
460  // Convert Python arguments to C
461  if (!PyArg_ParseTuple(args, "O", &quantity))
462    return NULL;
463
464  domain = PyObject_GetAttrString(quantity, "domain");
465  if (!domain)
466    return NULL;
467
468  //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
469  neighbours = get_consecutive_array(domain, "neighbours");
470
471  //Get safety factor beta_w
472  Tmp = PyObject_GetAttrString(domain, "beta_w");
473  if (!Tmp)
474    return NULL;
475
476  beta_w = PyFloat_AsDouble(Tmp);
477
478  Py_DECREF(Tmp);
479  Py_DECREF(domain);
480
481
482  qc = get_consecutive_array(quantity, "centroid_values");
483  qv = get_consecutive_array(quantity, "vertex_values");
484
485
486  N = qc -> dimensions[0];
487
488  //Find min and max of this and neighbour's centroid values
489  qmin = malloc(N * sizeof(double));
490  qmax = malloc(N * sizeof(double));
491  for (k=0; k<N; k++) {
492    k3=k*3;
493
494    qmin[k] = ((double*) qc -> data)[k];
495    qmax[k] = qmin[k];
496
497    for (i=0; i<3; i++) {
498      n = ((long*) neighbours -> data)[k3+i];
499      if (n >= 0) {
500        qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
501
502        qmin[k] = min(qmin[k], qn);
503        qmax[k] = max(qmax[k], qn);
504      }
505    }
506  }
507
508  // Call underlying routine
509  _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
510
511  free(qmin);
512  free(qmax);
513  return Py_BuildValue("");
514}
515
516
517
518// Method table for python module
519static struct PyMethodDef MethodTable[] = {
520  {"limit", limit, METH_VARARGS, "Print out"},
521  {"update", update, METH_VARARGS, "Print out"},
522  {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
523  {"extrapolate_second_order", extrapolate_second_order,
524   METH_VARARGS, "Print out"},
525  {"interpolate_from_vertices_to_edges",
526   interpolate_from_vertices_to_edges,
527   METH_VARARGS, "Print out"},
528  {NULL, NULL, 0, NULL}   // sentinel
529};
530
531// Module initialisation
532void initquantity_ext(void){
533  Py_InitModule("quantity_ext", MethodTable);
534
535  import_array();     //Necessary for handling of NumPY structures
536}
Note: See TracBrowser for help on using the repository browser.