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

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

Updated some methods for quantity. Looks like we can use old
limiting system with larger values of beta.

File size: 7.1 KB
Line 
1// Python - C extension module for advection.py
2//
3// To compile (Python2.X):
4// use python ../utilities/compile.py to
5// compile all C files within a directory
6//
7//
8// Steve Roberts, ANU 2008
9
10
11#include "Python.h"
12#include "Numeric/arrayobject.h"
13#include "math.h"
14#include "stdio.h"
15
16// Shared code snippets
17#include "util_ext.h"
18
19
20//-------------------------------------------
21// Low level routines (called from wrappers)
22//------------------------------------------
23
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,
36                    double  huge_timestep,
37                    double  max_timestep,
38                    int ntri,
39                    int nbdry){
40
41 
42        //Local Variables
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;
51        int  k_i,n_m,k_i_j;
52        int k,i,j,n,m;
53        int k3;
54       
55        //Loop through triangles
56
57        timestep = max_timestep;
58
59        for (k=0; k<ntri; k++){
60            optimal_timestep = huge_timestep;
61            flux = 0.0;
62            k3 = 3*k;
63            for (i=0; i<3; i++){
64                k_i = k3+i;
65                //Quantities inside triangle facing neighbour i
66                ql = quantity_edge[k_i];
67
68
69                //Quantities at neighbour on nearest face
70                n = domain_neighbours[k_i];
71                if (n < 0) {
72                    m = -n-1; //Convert neg flag to index
73                    qr = quantity_bdry[m];
74                } else {
75                    m = domain_neighbour_edges[k_i];
76                    n_m = 3*n+m;
77                    qr = quantity_edge[n_m];
78                }
79
80
81                //Outward pointing normal vector
82                for (j=0; j<2; j++){
83                    k_i_j = 6*k+2*i+j;
84                    normal[j] = domain_normals[k_i_j];
85                }
86
87
88                //Flux computation using provided function
89                normal_velocity = domain_velocity[0]*normal[0] + domain_velocity[1]*normal[1];
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);
98                flux = flux - edgeflux * domain_edgelengths[k_i];
99
100                //Update optimal_timestep
101                if (domain_tri_full_flag[k] == 1) {
102                    if (max_speed != 0.0) {
103                        optimal_timestep = (optimal_timestep>domain_radii[k]/max_speed) ? 
104                          domain_radii[k]/max_speed : optimal_timestep;
105                    }
106                }
107
108            }
109
110            //Normalise by area and store for when all conserved
111            //quantities get updated
112            quantity_update[k] = flux/domain_areas[k];
113
114            timestep = (timestep>optimal_timestep) ? optimal_timestep : timestep;
115
116        }
117
118        return timestep;
119}
120
121
122//-----------------------------------------------------
123// Python method Wrappers
124//-----------------------------------------------------
125
126
127
128PyObject *compute_fluxes(PyObject *self, PyObject *args) {
129  /*Compute all fluxes and the timestep suitable for all volumes
130    in domain.
131
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.
135
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;
150 
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);
227}
228
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.