source: anuga_core/source/anuga/advection/advection_ext.c @ 4978

Last change on this file since 4978 was 4978, checked in by steve, 16 years ago

Added a C extension for the advection directory (gave up on using f2py!)

File size: 6.9 KB
Line 
1// Python - C extension module for advection.py
2//
3// To compile (Python2.X):
4// use python ../utilities/compile.py
5//
6// See the module advection.py
7//
8//
9// Steve Roberts, ANU 2008
10
11
12#include "Python.h"
13#include "Numeric/arrayobject.h"
14#include "math.h"
15#include "stdio.h"
16
17// Shared code snippets
18#include "util_ext.h"
19
20
21
22double _compute_fluxes(
23                    double* quantity_update,
24                    double* quantity_edge,
25                    double* quantity_bdry,
26                    long*   domain_neighbours,
27                    long*   domain_neighbour_edges,
28                    double* domain_normals,
29                    double* domain_areas,
30                    double* domain_radii,
31                    double* domain_edgelengths,
32                    long*   domain_tri_full_flag,
33                    double* domain_velocity,
34                    double  huge_timestep,
35                    double  max_timestep,
36                    int ntri,
37                    int nbdry){
38
39 
40        //Loop
41
42        double qr,ql;
43        double normal[2];
44        double normal_velocity;
45        double flux, edgeflux;
46        double max_speed;
47        double optimal_timestep;
48        double timestep;
49        int  k_i,n_m,k_i_j;
50        int k,i,j,n,m;
51        int k3;
52       
53        timestep = max_timestep;
54
55        for (k=0; k<ntri; k++){
56            optimal_timestep = huge_timestep;
57            flux = 0.0;
58            k3 = 3*k;
59            for (i=0; i<3; i++){
60                k_i = k3+i;
61                //Quantities inside volume facing neighbour i
62                ql = quantity_edge[k_i];
63
64
65                //Quantities at neighbour on nearest face
66                n = domain_neighbours[k_i];
67                if (n < 0) {
68                    m = -n-1; //Convert neg flag to index
69                    qr = quantity_bdry[m];
70                } else {
71                    m = domain_neighbour_edges[k_i];
72                    n_m = 3*n+m;
73                    qr = quantity_edge[n_m];
74                }
75
76
77                //Outward pointing normal vector
78                for (j=0; j<2; j++){
79                    k_i_j = 6*k+2*i+j;
80                    normal[j] = domain_normals[k_i_j];
81                }
82
83
84                //Flux computation using provided function
85                normal_velocity = domain_velocity[0]*normal[0] + domain_velocity[1]*normal[1];
86
87                if (normal_velocity < 0) {
88                    edgeflux = qr * normal_velocity;
89                } else {
90                    edgeflux = ql * normal_velocity;
91                }
92
93                max_speed = fabs(normal_velocity);
94                flux = flux - edgeflux * domain_edgelengths[k_i];
95
96                //Update optimal_timestep
97                if (domain_tri_full_flag[k] == 1) {
98                    if (max_speed != 0.0) {
99                        optimal_timestep = (optimal_timestep>domain_radii[k]/max_speed) ? 
100                          domain_radii[k]/max_speed : optimal_timestep;
101                    }
102                }
103
104            }
105
106            //Normalise by area and store for when all conserved
107            //quantities get updated
108            quantity_update[k] = flux/domain_areas[k];
109
110            timestep = (timestep>optimal_timestep) ? optimal_timestep : timestep;
111
112        }
113
114        return timestep;
115}
116
117
118//-----------------------------------------------------
119// Python method Wrappers
120//-----------------------------------------------------
121
122
123
124PyObject *compute_fluxes(PyObject *self, PyObject *args) {
125  /*Compute all fluxes and the timestep suitable for all volumes
126    in domain.
127
128    Fluxes across each edge are scaled by edgelengths and summed up
129    Resulting flux is then scaled by area and stored in
130    explicit_update for the conserved quantity stage.
131
132    The maximal allowable speed computed for each volume
133    is converted to a timestep that must not be exceeded. The minimum of
134    those is computed as the next overall timestep.
135
136    Python call:
137    timestep = advection_ext.compute_fluxes(domain,quantity,huge_timestep,max_timestep)
138
139    Post conditions:
140      domain.explicit_update is reset to computed flux values
141      domain.timestep is set to the largest step satisfying all volumes.
142
143  */
144
145  PyObject *domain, *quantity;
146  PyArrayObject
147    * quantity_update,
148    * quantity_edge,
149    * quantity_bdry,
150    * domain_neighbours,
151    * domain_neighbour_edges,
152    * domain_normals,
153    * domain_areas,
154    * domain_radii,
155    * domain_edgelengths,
156    * domain_tri_full_flag,
157    * domain_velocity;
158                 
159
160  // Local variables
161  int ntri, nbdry;
162  double huge_timestep, max_timestep;
163  double timestep;
164
165
166  // Convert Python arguments to C
167  if (!PyArg_ParseTuple(args, "OOdd",
168                        &domain,
169                        &quantity,
170                        &huge_timestep,
171                        &max_timestep)) {
172    PyErr_SetString(PyExc_RuntimeError, "advection_ext.c: compute_fluxes could not parse input");
173    return NULL;
174  }
175
176  quantity_edge          = get_consecutive_array(quantity, "edge_values");
177  quantity_bdry          = get_consecutive_array(quantity, "boundary_values");
178  quantity_update        = get_consecutive_array(quantity, "explicit_update");
179  domain_neighbours      = get_consecutive_array(domain,   "neighbours");
180  domain_neighbour_edges = get_consecutive_array(domain,   "neighbour_edges");
181  domain_normals         = get_consecutive_array(domain,   "normals");
182  domain_areas           = get_consecutive_array(domain,   "areas");
183  domain_radii           = get_consecutive_array(domain,   "radii");
184  domain_edgelengths     = get_consecutive_array(domain,   "edgelengths");
185  domain_tri_full_flag   = get_consecutive_array(domain,   "tri_full_flag");
186  domain_velocity        = get_consecutive_array(domain,   "velocity"); 
187   
188  ntri  = quantity_edge -> dimensions[0];
189  nbdry = quantity_bdry -> dimensions[0]; 
190
191  timestep = _compute_fluxes((double*) quantity_update -> data,
192                             (double*) quantity_edge -> data,
193                             (double*) quantity_bdry -> data,
194                             (long*)   domain_neighbours -> data,
195                             (long*)   domain_neighbour_edges -> data,
196                             (double*) domain_normals -> data,
197                             (double*) domain_areas ->data,
198                             (double*) domain_radii -> data,
199                             (double*) domain_edgelengths -> data,
200                             (long*)   domain_tri_full_flag -> data,
201                             (double*) domain_velocity -> data,
202                             huge_timestep,
203                             max_timestep,
204                             ntri,
205                             nbdry);
206
207  // Release and return
208  Py_DECREF(quantity_update);
209  Py_DECREF(quantity_edge);
210  Py_DECREF(quantity_bdry);
211  Py_DECREF(domain_neighbours);
212  Py_DECREF(domain_neighbour_edges);
213  Py_DECREF(domain_normals);
214  Py_DECREF(domain_areas);
215  Py_DECREF(domain_radii);
216  Py_DECREF(domain_edgelengths);
217  Py_DECREF(domain_tri_full_flag);
218  Py_DECREF(domain_velocity);
219
220
221  return Py_BuildValue("d", timestep);
222}
223
224
225
226//-------------------------------
227// Method table for python module
228//-------------------------------
229static struct PyMethodDef MethodTable[] = {
230  /* The cast of the function is necessary since PyCFunction values
231   * only take two PyObject* parameters, and rotate() takes
232   * three.
233   */
234
235  {"compute_fluxes", compute_fluxes, METH_VARARGS, "Print out"},
236  {NULL, NULL}
237};
238
239// Module initialisation
240void initadvection_ext(void){
241  Py_InitModule("advection_ext", MethodTable);
242  import_array(); // Necessary for handling of NumPY structures
243}
Note: See TracBrowser for help on using the repository browser.