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

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

Added in Sudi's most recent 1d codes
File size: 9.5 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//====================================================================
323PyObject *limit_vanalbada_ext(PyObject *self, PyObject *args) {
324
325   PyObject
326           *domain,
327           *quantity;
328
329    PyArrayObject
330        *qco,
331        *qvo,
332        *xco,
333        *xvo;
334
335    double *qc,
336            *qv,
337            *xc,
338            *xv;
339
340  double a, b;
341  //double c;
342  double phi, dx0, dx1;
343  double theta;
344  int N, k, k2;
345
346
347  // Convert Python arguments to C
348  if (!PyArg_ParseTuple(args, "O", &quantity)) {
349      PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: limit_vanalbada_ext could not parse input");
350      return NULL;
351  }
352
353
354  domain = PyObject_GetAttrString(quantity, "domain");
355
356  //printf("B = %p\n",(void*)domain);
357  if (!domain) {
358    printf("quantity_ext.c: Could not obtain python object");
359    fflush(stdout);
360    PyErr_SetString(PyExc_RuntimeError, "quantity_ext.c: Could not obtain python object domain");
361    return NULL;
362  }
363
364
365
366  N  = get_python_integer(quantity,"N");
367  theta = get_python_double(quantity,"beta");
368
369
370  //printf("beta = %f",theta);
371  //fflush(stdout);
372
373  qco = get_consecutive_array(quantity, "centroid_values");
374  qvo = get_consecutive_array(quantity, "vertex_values");
375  xco = get_consecutive_array(domain, "centroids");
376  xvo = get_consecutive_array(domain, "vertices");
377
378  qc = (double *) qco -> data;
379  qv = (double *) qvo -> data;
380  xc = (double *) xco -> data;
381  xv = (double *) xvo -> data;
382
383
384  for (k=0; k<N; k++) {
385      k2 = 2*k;
386      if (k == 0) {
387        phi = (qc[1]-qc[0])/(xc[1] - xc[0]);
388      }
389      else if (k==N-1) {
390          phi = (qc[N-1] - qc[N-2])/(xc[N-1] - xc[N-2]);
391      }
392      else {
393          a  = (qc[k]-qc[k-1])/(xc[k]-xc[k-1]);
394          b  = (qc[k+1]-qc[k])/(xc[k+1]-xc[k]);
395          //c  = (qc[k+1]-qc[k-1])/(xc[k+1]-xc[k-1]);
396
397
398          phi = 0.0;
399          if (a*a + b*b >= 1.0e-32) {
400              phi = (a*a*b+a*b*b)/(a*a+b*b);
401          }
402
403
404      }
405
406      dx0 = xv[k2] - xc[k];
407      dx1 = xv[k2+1] - xc[k];
408
409
410      qv[k2]   = qc[k] + phi*dx0;
411      qv[k2+1] = qc[k] + phi*dx1;
412  }
413
414
415  Py_DECREF(qco);
416  Py_DECREF(qvo);
417  Py_DECREF(xco);
418  Py_DECREF(xvo);
419
420  // Return updated flux timestep
421  return Py_BuildValue("");
422}
423
424
425
426//-------------------------------
427// Method table for python module
428//-------------------------------
429
430static struct PyMethodDef MethodTable[] = {
431  {"limit_minmod_ext",          limit_minmod_ext, METH_VARARGS, "Print out"},
432  {"limit_minmod_kurganov_ext", limit_minmod_kurganov_ext, METH_VARARGS, "Print out"},
433  {"limit_vanleer_ext",         limit_vanleer_ext, METH_VARARGS, "Print out"},
434  {"limit_vanalbada_ext",       limit_vanalbada_ext, METH_VARARGS, "Print out"},
435  {NULL, NULL}
436};
437
438// Module initialisation
439void initquantity_ext(void){
440  Py_InitModule("quantity_ext", MethodTable);
441  import_array();
442}
Note: See TracBrowser for help on using the repository browser.