source: trunk/anuga_work/development/2010-projects/anuga_1d/base/quantity_ext.c @ 7884

Last change on this file since 7884 was 7884, checked in by steve, 14 years ago

Moving 2010 project

File size: 7.3 KB
Line 
1#include "Python.h"
2#include "numpy/arrayobject.h"
3#include "math.h"
4#include <stdio.h>
5const double pi = 3.14159265358979;
6
7
8// Shared code snippets
9#include "util_ext.h"
10
11
12
13
14
15
16
17
18//=========================================================================
19// Python Glue
20//=========================================================================
21
22PyObject *limit_minmod_ext(PyObject *self, PyObject *args) {
23 
24   PyObject
25           *domain,
26           *quantity;
27
28    PyArrayObject
29        *qco,
30        *qvo,
31        *xco,
32        *xvo;
33
34    double *qc,
35            *qv,
36            *xc,
37            *xv;
38   
39  double a, b;
40  double phi, dx0, dx1;
41  int N, k, k2;
42
43   
44  // Convert Python arguments to C
45  if (!PyArg_ParseTuple(args, "O", &quantity)) {
46      PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_minmod_ext could not parse input");
47      return NULL;
48  }
49 
50
51  domain = PyObject_GetAttrString(quantity, "domain");
52
53  //printf("B = %p\n",(void*)domain);
54  if (!domain) {
55    printf("quantity_ext.c: Could not obtain python object");
56    fflush(stdout);
57    PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: Could not obtain python object domain");
58    return NULL;
59  }
60 
61
62 
63  N  = get_python_integer(quantity,"N");
64 
65  qco = get_consecutive_array(quantity, "centroid_values");
66  qvo = get_consecutive_array(quantity, "vertex_values");
67  xco = get_consecutive_array(domain, "centroids");
68  xvo = get_consecutive_array(domain, "vertices");
69   
70  qc = (double *) qco -> data;
71  qv = (double *) qvo -> data;
72  xc = (double *) xco -> data;
73  xv = (double *) xvo -> data;
74
75
76  for (k=0; k<N; k++) {
77      k2 = 2*k;
78      if (k == 0) {
79        phi = (qc[1]-qc[0])/(xc[1] - xc[0]);
80      }
81      else if (k==N-1) {
82          phi = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2]);
83      }
84      else {
85          a  = (qc[k]-qc[k-1])/(xc[k]-xc[k-1]);
86          b  = (qc[k+1]-qc[k])/(xc[k+1]-xc[k]);
87          //c  = (qc[K+1]-qc[k-1])/(xc[k+1]-xc[k-1]);
88
89          phi = 0.0;
90          if ((fabs(a) < fabs(b)) & (a*b > 0.0 )) {
91              phi = a;
92          }
93          if ((fabs(b) < fabs(a)) & (a*b > 0.0 )) {
94              phi = b;
95          }
96
97      }
98     
99      dx0 = xv[k2] - xc[k];
100      dx1 = xv[k2+1] - xc[k];
101
102
103      qv[k2]   = qc[k] + phi*dx0;
104      qv[k2+1] = qc[k] + phi*dx1;
105  }
106
107
108  Py_DECREF(qco);
109  Py_DECREF(qvo);
110  Py_DECREF(xco);
111  Py_DECREF(xvo);
112
113  // Return updated flux timestep
114  return Py_BuildValue("");
115}
116
117
118//====================================================================
119PyObject *limit_minmod_kurganov_ext(PyObject *self, PyObject *args) {
120
121   PyObject
122           *domain,
123           *quantity;
124
125    PyArrayObject
126        *qco,
127        *qvo,
128        *xco,
129        *xvo;
130
131    double *qc,
132            *qv,
133            *xc,
134            *xv;
135
136  double a, b, c;
137  double phi, dx0, dx1;
138  double theta;
139  int N, k, k2;
140
141
142  // Convert Python arguments to C
143  if (!PyArg_ParseTuple(args, "O", &quantity)) {
144      PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_ minmod_kurganov_ext could not parse input");
145      return NULL;
146  }
147
148
149  domain = PyObject_GetAttrString(quantity, "domain");
150
151  //printf("B = %p\n",(void*)domain);
152  if (!domain) {
153    printf("quantity_ext.c: Could not obtain python object");
154    fflush(stdout);
155    PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: Could not obtain python object domain");
156    return NULL;
157  }
158
159
160
161  N  = get_python_integer(quantity,"N");
162  theta = get_python_double(quantity,"beta");
163
164
165  //printf("beta = %f",theta);
166  //fflush(stdout);
167
168  qco = get_consecutive_array(quantity, "centroid_values");
169  qvo = get_consecutive_array(quantity, "vertex_values");
170  xco = get_consecutive_array(domain, "centroids");
171  xvo = get_consecutive_array(domain, "vertices");
172
173  qc = (double *) qco -> data;
174  qv = (double *) qvo -> data;
175  xc = (double *) xco -> data;
176  xv = (double *) xvo -> data;
177
178
179  for (k=0; k<N; k++) {
180      k2 = 2*k;
181      if (k == 0) {
182        phi = (qc[1]-qc[0])/(xc[1] - xc[0]);
183      }
184      else if (k==N-1) {
185          phi = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2]);
186      }
187      else {
188          a  = (qc[k]-qc[k-1])/(xc[k]-xc[k-1]);
189          b  = (qc[k+1]-qc[k])/(xc[k+1]-xc[k]);
190          c  = (qc[k+1]-qc[k-1])/(xc[k+1]-xc[k-1]);
191
192          phi = 0.0;
193          if ((sign(a)*sign(b) > 0.0) & (sign(a)*sign(c) > 0.0 )) {
194              phi = sign(a)*min(theta*min(fabs(a),fabs(b)),fabs(c));
195          }
196
197
198      }
199
200      dx0 = xv[k2] - xc[k];
201      dx1 = xv[k2+1] - xc[k];
202
203
204      qv[k2]   = qc[k] + phi*dx0;
205      qv[k2+1] = qc[k] + phi*dx1;
206  }
207
208
209  Py_DECREF(qco);
210  Py_DECREF(qvo);
211  Py_DECREF(xco);
212  Py_DECREF(xvo);
213
214  // Return updated flux timestep
215  return Py_BuildValue("");
216}
217
218 
219//====================================================================
220PyObject *limit_vanleer_ext(PyObject *self, PyObject *args) {
221
222   PyObject
223           *domain,
224           *quantity;
225
226    PyArrayObject
227        *qco,
228        *qvo,
229        *xco,
230        *xvo;
231
232    double *qc,
233            *qv,
234            *xc,
235            *xv;
236
237  double a, b;
238  double phi, dx0, dx1;
239  double theta;
240  int N, k, k2;
241
242
243  // Convert Python arguments to C
244  if (!PyArg_ParseTuple(args, "O", &quantity)) {
245      PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_vanleer_ext could not parse input");
246      return NULL;
247  }
248
249
250  domain = PyObject_GetAttrString(quantity, "domain");
251
252  //printf("B = %p\n",(void*)domain);
253  if (!domain) {
254    printf("quantity_ext.c: Could not obtain python object");
255    fflush(stdout);
256    PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: Could not obtain python object domain");
257    return NULL;
258  }
259
260
261
262  N  = get_python_integer(quantity,"N");
263  theta = get_python_double(quantity,"beta");
264
265
266  //printf("beta = %f",theta);
267  //fflush(stdout);
268
269  qco = get_consecutive_array(quantity, "centroid_values");
270  qvo = get_consecutive_array(quantity, "vertex_values");
271  xco = get_consecutive_array(domain, "centroids");
272  xvo = get_consecutive_array(domain, "vertices");
273
274  qc = (double *) qco -> data;
275  qv = (double *) qvo -> data;
276  xc = (double *) xco -> data;
277  xv = (double *) xvo -> data;
278
279
280  for (k=0; k<N; k++) {
281      k2 = 2*k;
282      if (k == 0) {
283        phi = (qc[1]-qc[0])/(xc[1] - xc[0]);
284      }
285      else if (k==N-1) {
286          phi = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2]);
287      }
288      else {
289          a  = (qc[k]-qc[k-1])/(xc[k]-xc[k-1]);
290          b  = (qc[k+1]-qc[k])/(xc[k+1]-xc[k]);
291          //c  = (qc[k+1]-qc[k-1])/(xc[k+1]-xc[k-1]);
292
293
294          phi = 0.0;
295          if ((fabs(a)+fabs(b)) > 1.0e-12) {
296              phi = (a*fabs(b)+fabs(a)*b)/(fabs(a)+fabs(b));
297          }
298          //printf("phi = %f",phi);
299          //fflush(stdout);
300
301      }
302
303      dx0 = xv[k2] - xc[k];
304      dx1 = xv[k2+1] - xc[k];
305
306
307      qv[k2]   = qc[k] + phi*dx0;
308      qv[k2+1] = qc[k] + phi*dx1;
309  }
310
311
312  Py_DECREF(qco);
313  Py_DECREF(qvo);
314  Py_DECREF(xco);
315  Py_DECREF(xvo);
316
317  // Return updated flux timestep
318  return Py_BuildValue("");
319}
320
321
322
323
324//-------------------------------
325// Method table for python module
326//-------------------------------
327
328static struct PyMethodDef MethodTable[] = {
329  {"limit_minmod_ext",          limit_minmod_ext, METH_VARARGS, "Print out"},
330  {"limit_minmod_kurganov_ext", limit_minmod_kurganov_ext, METH_VARARGS, "Print out"},
331  {"limit_vanleer_ext",         limit_vanleer_ext, METH_VARARGS, "Print out"},
332  {NULL, NULL}
333};
334
335// Module initialisation
336void initquantity_ext(void){
337  Py_InitModule("quantity_ext", MethodTable);
338  import_array();
339}
Note: See TracBrowser for help on using the repository browser.