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 "numpy/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 | |
---|
24 | double _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 | |
---|
128 | PyObject *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 | //------------------------------- |
---|
234 | static 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 |
---|
245 | void initadvection_ext(void){ |
---|
246 | Py_InitModule("advection_ext", MethodTable); |
---|
247 | import_array(); // Necessary for handling of NumPY structures |
---|
248 | } |
---|