[7884] | 1 | #include "Python.h" |
---|
| 2 | #include "numpy/arrayobject.h" |
---|
| 3 | #include "math.h" |
---|
| 4 | #include <stdio.h> |
---|
| 5 | const 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 | |
---|
| 22 | PyObject *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 | //==================================================================== |
---|
| 119 | PyObject *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 | //==================================================================== |
---|
| 220 | PyObject *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 | |
---|
| 328 | static 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 |
---|
| 336 | void initquantity_ext(void){ |
---|
| 337 | Py_InitModule("quantity_ext", MethodTable); |
---|
| 338 | import_array(); |
---|
| 339 | } |
---|