source: branches/numpy/anuga/advection/advection_ext.c @ 6883

Last change on this file since 6883 was 6304, checked in by rwilson, 16 years ago

Initial commit of numpy changes. Still a long way to go.

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