source: anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c @ 4536

Last change on this file since 4536 was 4536, checked in by ole, 17 years ago

Fixed ticket:165 thanks to Duncan's excellent reduction of the problem which
had to do with ghost nodes.

File size: 16.0 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        //Semi implicit updates
190        for (k=0; k<N; k++) {
191                denominator = 1.0 - timestep*semi_implicit_update[k];
192                if (denominator == 0.0) {
193                        return -1;
194                } else {
195                        //Update conserved_quantities from semi implicit updates
196                        centroid_values[k] /= denominator;
197                }
198        }
199
200        /*  for (k=0; k<N; k++) {*/
201        /*    centroid_values[k] = exp(timestep*semi_implicit_update[k])*centroid_values[k];*/
202        /*  }*/
203
204
205        //Explicit updates
206        for (k=0; k<N; k++) {
207                centroid_values[k] += timestep*explicit_update[k];
208        }
209
210
211        //MH080605 set semi_implicit_update[k] to 0.0 here, rather than in update_conserved_quantities.py
212        for (k=0;k<N;k++){
213                semi_implicit_update[k]=0.0;
214        }
215
216        return 0;
217}
218
219
220int _average_vertex_values(int N,
221                           long* vertex_value_indices,
222                           long* number_of_triangles_per_node,
223                           double* vertex_values, 
224                           double* A) {
225  //Average vertex values to obtain one value per node
226
227  int i, index; 
228  int k = 0; //Track triangles touching each node
229  int current_node = 0;
230  double total = 0.0;
231
232  for (i=0; i<N; i++) {
233 
234    //if (current_node == N) {
235    //  printf("Current node exceeding number of nodes (%d)", N);   
236    //  return 1;
237    // }
238   
239    index = vertex_value_indices[i];
240    k += 1;
241           
242    //volume_id = index / 3
243    //vertex_id = index % 3
244    //total += self.vertex_values[volume_id, vertex_id]
245    total += vertex_values[index];
246     
247    //printf("current_node=%d, index=%d, k=%d, total=%f\n", current_node, index, k, total);
248    if (number_of_triangles_per_node[current_node] == k) {
249      A[current_node] = total/k;
250               
251      // Move on to next node
252      total = 0.0;
253      k = 0;
254      current_node += 1;
255    }
256  }
257                           
258  return 0;
259}
260
261
262/////////////////////////////////////////////////
263// Gateways to Python
264PyObject *update(PyObject *self, PyObject *args) {
265
266        PyObject *quantity;
267        PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
268
269        double timestep;
270        int N, err;
271
272
273        // Convert Python arguments to C
274        if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep)) {
275          PyErr_SetString(PyExc_RuntimeError, 
276                          "quantity_ext.c: update could not parse input");
277          return NULL;
278        }
279
280        centroid_values = get_consecutive_array(quantity, "centroid_values");
281        explicit_update = get_consecutive_array(quantity, "explicit_update");
282        semi_implicit_update = get_consecutive_array(quantity, "semi_implicit_update");
283
284        N = centroid_values -> dimensions[0];
285
286        err = _update(N, timestep,
287                      (double*) centroid_values -> data,
288                      (double*) explicit_update -> data,
289                      (double*) semi_implicit_update -> data);
290
291
292        if (err != 0) {
293          PyErr_SetString(PyExc_RuntimeError,
294                          "Zero division in semi implicit update - call Stephen :)");
295          return NULL;
296        }
297
298        //Release and return
299        Py_DECREF(centroid_values);
300        Py_DECREF(explicit_update);
301        Py_DECREF(semi_implicit_update);
302
303        return Py_BuildValue("");
304}
305
306
307PyObject *interpolate_from_vertices_to_edges(PyObject *self, PyObject *args) {
308
309        PyObject *quantity;
310        PyArrayObject *vertex_values, *edge_values;
311
312        int N, err;
313
314        // Convert Python arguments to C
315        if (!PyArg_ParseTuple(args, "O", &quantity)) {
316          PyErr_SetString(PyExc_RuntimeError, 
317                          "quantity_ext.c: interpolate_from_vertices_to_edges could not parse input");
318          return NULL;
319        }
320       
321        vertex_values = get_consecutive_array(quantity, "vertex_values");
322        edge_values = get_consecutive_array(quantity, "edge_values");
323
324        N = vertex_values -> dimensions[0];
325
326        err = _interpolate(N,
327                           (double*) vertex_values -> data,
328                           (double*) edge_values -> data);
329
330        if (err != 0) {
331          PyErr_SetString(PyExc_RuntimeError, 
332                          "Interpolate could not be computed");
333          return NULL;
334        }
335
336        //Release and return
337        Py_DECREF(vertex_values);
338        Py_DECREF(edge_values);
339
340        return Py_BuildValue("");
341}
342
343
344PyObject *average_vertex_values(PyObject *self, PyObject *args) {
345
346        PyArrayObject
347          *vertex_value_indices, 
348          *number_of_triangles_per_node,
349          *vertex_values,
350          *A;
351       
352
353        int N, err;
354
355        // Convert Python arguments to C
356        if (!PyArg_ParseTuple(args, "OOOO",
357                              &vertex_value_indices, 
358                              &number_of_triangles_per_node,
359                              &vertex_values, 
360                              &A)) {
361          PyErr_SetString(PyExc_RuntimeError, 
362                          "quantity_ext.c: average_vertex_values could not parse input");
363          return NULL;
364        }
365       
366        N = vertex_value_indices -> dimensions[0];
367        // printf("Got parameters, N=%d\n", N);
368        err = _average_vertex_values(N,
369                                     (long*) vertex_value_indices -> data,
370                                     (long*) number_of_triangles_per_node -> data,
371                                     (double*) vertex_values -> data, 
372                                     (double*) A -> data);
373
374        //printf("Error %d", err);
375        if (err != 0) {
376          PyErr_SetString(PyExc_RuntimeError, 
377                          "average_vertex_values could not be computed");
378          return NULL;
379        }
380
381        return Py_BuildValue("");
382}
383
384
385
386PyObject *compute_gradients(PyObject *self, PyObject *args) {
387
388        PyObject *quantity, *domain, *R;
389        PyArrayObject
390                *centroids,            //Coordinates at centroids
391                *centroid_values,      //Values at centroids
392                *number_of_boundaries, //Number of boundaries for each triangle
393                *surrogate_neighbours, //True neighbours or - if one missing - self
394                *a, *b;                //Return values
395
396        int dimensions[1], N, err;
397
398        // Convert Python arguments to C
399        if (!PyArg_ParseTuple(args, "O", &quantity)) {
400          PyErr_SetString(PyExc_RuntimeError, 
401                          "quantity_ext.c: compute_gradients could not parse input");   
402          return NULL;
403        }
404
405        domain = PyObject_GetAttrString(quantity, "domain");
406        if (!domain) {
407          PyErr_SetString(PyExc_RuntimeError, 
408                          "compute_gradients could not obtain domain object from quantity");
409          return NULL;
410        }
411
412        //Get pertinent variables
413
414        centroids = get_consecutive_array(domain, "centroid_coordinates");
415        centroid_values = get_consecutive_array(quantity, "centroid_values");
416        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
417        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
418
419        N = centroid_values -> dimensions[0];
420
421        //Release
422        Py_DECREF(domain);
423
424        //Allocate space for return vectors a and b (don't DECREF)
425        dimensions[0] = N;
426        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
427        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
428
429
430
431        err = _compute_gradients(N,
432                        (double*) centroids -> data,
433                        (double*) centroid_values -> data,
434                        (long*) number_of_boundaries -> data,
435                        (long*) surrogate_neighbours -> data,
436                        (double*) a -> data,
437                        (double*) b -> data);
438
439        if (err != 0) {
440          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
441          return NULL;
442        }
443
444        //Release
445        Py_DECREF(centroids);
446        Py_DECREF(centroid_values);
447        Py_DECREF(number_of_boundaries);
448        Py_DECREF(surrogate_neighbours);
449
450        //Build result, release and return
451        R = Py_BuildValue("OO", PyArray_Return(a), PyArray_Return(b));
452        Py_DECREF(a);
453        Py_DECREF(b);
454        return R;
455}
456
457
458
459PyObject *extrapolate_second_order(PyObject *self, PyObject *args) {
460
461        PyObject *quantity, *domain;
462        PyArrayObject
463                *centroids,            //Coordinates at centroids
464                *centroid_values,      //Values at centroids
465                *vertex_coordinates,   //Coordinates at vertices
466                *vertex_values,        //Values at vertices
467                *number_of_boundaries, //Number of boundaries for each triangle
468                *surrogate_neighbours, //True neighbours or - if one missing - self
469                *a, *b;                //Gradients
470
471        //int N, err;
472        int dimensions[1], N, err;
473        //double *a, *b;  //Gradients
474
475        // Convert Python arguments to C
476        if (!PyArg_ParseTuple(args, "O", &quantity)) {
477          PyErr_SetString(PyExc_RuntimeError, 
478                          "extrapolate_second_order could not parse input");   
479          return NULL;
480        }
481
482        domain = PyObject_GetAttrString(quantity, "domain");
483        if (!domain) {
484          PyErr_SetString(PyExc_RuntimeError, 
485                          "extrapolate_second_order could not obtain domain object from quantity");     
486          return NULL;
487        }
488
489        //Get pertinent variables
490        centroids = get_consecutive_array(domain, "centroid_coordinates");
491        centroid_values = get_consecutive_array(quantity, "centroid_values");
492        surrogate_neighbours = get_consecutive_array(domain, "surrogate_neighbours");
493        number_of_boundaries = get_consecutive_array(domain, "number_of_boundaries");
494        vertex_coordinates = get_consecutive_array(domain, "vertex_coordinates");
495        vertex_values = get_consecutive_array(quantity, "vertex_values");
496
497        N = centroid_values -> dimensions[0];
498
499        //Release
500        Py_DECREF(domain);
501
502        //Allocate space for return vectors a and b (don't DECREF)
503        dimensions[0] = N;
504        a = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
505        b = (PyArrayObject *) PyArray_FromDims(1, dimensions, PyArray_DOUBLE);
506
507        //FIXME: Odd that I couldn't use normal arrays
508        //Allocate space for return vectors a and b (don't DECREF)
509        //a = (double*) malloc(N * sizeof(double));
510        //if (!a) return NULL;
511        //b = (double*) malloc(N * sizeof(double));
512        //if (!b) return NULL;
513
514
515        err = _compute_gradients(N,
516                        (double*) centroids -> data,
517                        (double*) centroid_values -> data,
518                        (long*) number_of_boundaries -> data,
519                        (long*) surrogate_neighbours -> data,
520                        (double*) a -> data,
521                        (double*) b -> data);
522
523        if (err != 0) {
524          PyErr_SetString(PyExc_RuntimeError, "Gradient could not be computed");
525          return NULL;
526        }
527
528        err = _extrapolate(N,
529                        (double*) centroids -> data,
530                        (double*) centroid_values -> data,
531                        (double*) vertex_coordinates -> data,
532                        (double*) vertex_values -> data,
533                        (double*) a -> data,
534                        (double*) b -> data);
535
536
537        if (err != 0) {
538          PyErr_SetString(PyExc_RuntimeError,
539                          "Internal function _extrapolate failed");
540          return NULL;
541        }
542
543
544
545        //Release
546        Py_DECREF(centroids);
547        Py_DECREF(centroid_values);
548        Py_DECREF(number_of_boundaries);
549        Py_DECREF(surrogate_neighbours);
550        Py_DECREF(vertex_coordinates);
551        Py_DECREF(vertex_values);
552        Py_DECREF(a);
553        Py_DECREF(b);
554
555        return Py_BuildValue("");
556}
557
558
559
560PyObject *limit(PyObject *self, PyObject *args) {
561
562        PyObject *quantity, *domain, *Tmp;
563        PyArrayObject
564                *qv, //Conserved quantities at vertices
565                *qc, //Conserved quantities at centroids
566                *neighbours;
567
568        int k, i, n, N, k3;
569        double beta_w; //Safety factor
570        double *qmin, *qmax, qn;
571
572        // Convert Python arguments to C
573        if (!PyArg_ParseTuple(args, "O", &quantity)) {
574          PyErr_SetString(PyExc_RuntimeError, 
575                          "quantity_ext.c: limit could not parse input");
576          return NULL;
577        }
578
579        domain = PyObject_GetAttrString(quantity, "domain");
580        if (!domain) {
581          PyErr_SetString(PyExc_RuntimeError, 
582                          "quantity_ext.c: limit could not obtain domain object from quantity");                       
583         
584          return NULL;
585        }
586
587        //neighbours = (PyArrayObject*) PyObject_GetAttrString(domain, "neighbours");
588        neighbours = get_consecutive_array(domain, "neighbours");
589
590        //Get safety factor beta_w
591        Tmp = PyObject_GetAttrString(domain, "beta_w");
592        if (!Tmp) {
593          PyErr_SetString(PyExc_RuntimeError, 
594                          "quantity_ext.c: limit could not obtain beta_w object from domain");                 
595         
596          return NULL;
597        }       
598
599        beta_w = PyFloat_AsDouble(Tmp);
600
601        Py_DECREF(Tmp);
602        Py_DECREF(domain);
603
604
605        qc = get_consecutive_array(quantity, "centroid_values");
606        qv = get_consecutive_array(quantity, "vertex_values");
607
608
609        N = qc -> dimensions[0];
610
611        //Find min and max of this and neighbour's centroid values
612        qmin = malloc(N * sizeof(double));
613        qmax = malloc(N * sizeof(double));
614        for (k=0; k<N; k++) {
615                k3=k*3;
616
617                qmin[k] = ((double*) qc -> data)[k];
618                qmax[k] = qmin[k];
619
620                for (i=0; i<3; i++) {
621                        n = ((long*) neighbours -> data)[k3+i];
622                        if (n >= 0) {
623                                qn = ((double*) qc -> data)[n]; //Neighbour's centroid value
624
625                                qmin[k] = min(qmin[k], qn);
626                                qmax[k] = max(qmax[k], qn);
627                        }
628                        //qmin[k] = max(qmin[k],0.5*((double*) qc -> data)[k]);
629                        //qmax[k] = min(qmax[k],2.0*((double*) qc -> data)[k]);
630                }
631        }
632
633        // Call underlying routine
634        _limit(N, beta_w, (double*) qc -> data, (double*) qv -> data, qmin, qmax);
635
636        free(qmin);
637        free(qmax);
638        return Py_BuildValue("");
639}
640
641
642
643// Method table for python module
644static struct PyMethodDef MethodTable[] = {
645        {"limit", limit, METH_VARARGS, "Print out"},
646        {"update", update, METH_VARARGS, "Print out"},
647        {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},
648        {"extrapolate_second_order", extrapolate_second_order,
649                METH_VARARGS, "Print out"},
650        {"interpolate_from_vertices_to_edges",
651                interpolate_from_vertices_to_edges,
652                METH_VARARGS, "Print out"},
653        {"average_vertex_values", average_vertex_values, METH_VARARGS, "Print out"},           
654        {NULL, NULL, 0, NULL}   // sentinel
655};
656
657// Module initialisation
658void initquantity_ext(void){
659  Py_InitModule("quantity_ext", MethodTable);
660
661  import_array();     //Necessary for handling of NumPY structures
662}
Note: See TracBrowser for help on using the repository browser.