 Timestamp:
 Dec 27, 2007, 6:35:08 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity_ext.c
r4886 r4897 137 137 138 138 139 int _limit_ by_vertex(int N, double beta,139 int _limit_vertices_by_all_neighbours(int N, double beta, 140 140 double* centroid_values, 141 141 double* vertex_values, … … 195 195 } 196 196 197 int _limit_ by_edge(int N, double beta,197 int _limit_edges_by_all_neighbours(int N, double beta, 198 198 double* centroid_values, 199 199 double* vertex_values, … … 240 240 241 241 //Update vertex and edge values using phi limiter 242 edge_values[k3+0] = qc + phi*dqa[0]; 243 edge_values[k3+1] = qc + phi*dqa[1]; 244 edge_values[k3+2] = qc + phi*dqa[2]; 245 246 vertex_values[k3+0] = edge_values[k3+1] + edge_values[k3+2]  edge_values[k3+0]; 247 vertex_values[k3+1] = edge_values[k3+2] + edge_values[k3+0]  edge_values[k3+1]; 248 vertex_values[k3+2] = edge_values[k3+0] + edge_values[k3+1]  edge_values[k3+2]; 249 250 } 251 252 return 0; 253 } 254 255 256 int _limit_edges_by_neighbour(int N, double beta, 257 double* centroid_values, 258 double* vertex_values, 259 double* edge_values, 260 long* neighbours) { 261 262 int i, k, k2, k3, k6; 263 long n; 264 double qmin, qmax, qn, qc; 265 double dq, dqa[3], phi, r; 266 267 for (k=0; k<N; k++){ 268 k6 = 6*k; 269 k3 = 3*k; 270 k2 = 2*k; 271 272 qc = centroid_values[k]; 273 phi = 1.0; 274 275 for (i=0; i<3; i++) { 276 dq = edge_values[k3+i]  qc; //Delta between edge and centroid values 277 dqa[i] = dq; //Save dq for use in updating vertex values 278 279 n = neighbours[k3+i]; 280 if (n >= 0) { 281 qn = centroid_values[n]; //Neighbour's centroid value 282 283 qmin = min(qc, qn); 284 qmax = max(qc, qn); 285 286 r = 1.0; 287 288 if (dq > 0.0) r = (qmax  qc)/dq; 289 if (dq < 0.0) r = (qmin  qc)/dq; 290 291 phi = min( min(r*beta, 1.0), phi); 292 } 293 } 294 295 296 //Update edge and vertex values using phi limiter 242 297 edge_values[k3+0] = qc + phi*dqa[0]; 243 298 edge_values[k3+1] = qc + phi*dqa[1]; … … 974 1029 975 1030 976 PyObject *limit_ by_vertex(PyObject *self, PyObject *args) {1031 PyObject *limit_vertices_by_all_neighbours(PyObject *self, PyObject *args) { 977 1032 //Limit slopes for each volume to eliminate artificial variance 978 1033 //introduced by e.g. second order extrapolator … … 1033 1088 N = centroid_values > dimensions[0]; 1034 1089 1035 err = _limit_ by_vertex(N, beta_w,1036 1037 1038 1039 1040 1090 err = _limit_vertices_by_all_neighbours(N, beta_w, 1091 (double*) centroid_values > data, 1092 (double*) vertex_values > data, 1093 (double*) edge_values > data, 1094 (long*) neighbours > data); 1095 1041 1096 if (err != 0) { 1042 1097 PyErr_SetString(PyExc_RuntimeError, … … 1059 1114 1060 1115 1061 PyObject *limit_ by_edge(PyObject *self, PyObject *args) {1116 PyObject *limit_edges_by_all_neighbours(PyObject *self, PyObject *args) { 1062 1117 //Limit slopes for each volume to eliminate artificial variance 1063 1118 //introduced by e.g. second order extrapolator … … 1086 1141 if (!PyArg_ParseTuple(args, "O", &quantity)) { 1087 1142 PyErr_SetString(PyExc_RuntimeError, 1088 "quantity_ext.c: limit_ by_edgecould not parse input");1143 "quantity_ext.c: limit_edges_by_all_neighbours could not parse input"); 1089 1144 return NULL; 1090 1145 } … … 1093 1148 if (!domain) { 1094 1149 PyErr_SetString(PyExc_RuntimeError, 1095 "quantity_ext.c: limit_ by_edgecould not obtain domain object from quantity");1150 "quantity_ext.c: limit_edges_by_all_neighbours could not obtain domain object from quantity"); 1096 1151 1097 1152 return NULL; … … 1102 1157 if (!Tmp) { 1103 1158 PyErr_SetString(PyExc_RuntimeError, 1104 "quantity_ext.c: limit_ by_edgecould not obtain beta_w object from domain");1159 "quantity_ext.c: limit_edges_by_all_neighbours could not obtain beta_w object from domain"); 1105 1160 1106 1161 return NULL; … … 1118 1173 N = centroid_values > dimensions[0]; 1119 1174 1120 err = _limit_ by_edge(N, beta_w,1121 1122 1123 1124 1175 err = _limit_edges_by_all_neighbours(N, beta_w, 1176 (double*) centroid_values > data, 1177 (double*) vertex_values > data, 1178 (double*) edge_values > data, 1179 (long*) neighbours > data); 1125 1180 1126 1181 if (err != 0) { … … 1143 1198 1144 1199 1200 PyObject *limit_edges_by_neighbour(PyObject *self, PyObject *args) { 1201 //Limit slopes for each volume to eliminate artificial variance 1202 //introduced by e.g. second order extrapolator 1203 1204 //This is an unsophisticated limiter as it does not take into 1205 //account dependencies among quantities. 1206 1207 //precondition: 1208 // vertex values are estimated from gradient 1209 //postcondition: 1210 // vertex and edge values are updated 1211 // 1212 1213 PyObject *quantity, *domain, *Tmp; 1214 PyArrayObject 1215 *vertex_values, //Conserved quantities at vertices 1216 *centroid_values, //Conserved quantities at centroids 1217 *edge_values, //Conserved quantities at edges 1218 *neighbours; 1219 1220 double beta_w; //Safety factor 1221 int N, err; 1222 1223 1224 // Convert Python arguments to C 1225 if (!PyArg_ParseTuple(args, "O", &quantity)) { 1226 PyErr_SetString(PyExc_RuntimeError, 1227 "quantity_ext.c: limit_edges_by_neighbour could not parse input"); 1228 return NULL; 1229 } 1230 1231 domain = PyObject_GetAttrString(quantity, "domain"); 1232 if (!domain) { 1233 PyErr_SetString(PyExc_RuntimeError, 1234 "quantity_ext.c: limit_edges_by_neighbour could not obtain domain object from quantity"); 1235 1236 return NULL; 1237 } 1238 1239 // Get safety factor beta_w 1240 Tmp = PyObject_GetAttrString(domain, "beta_w"); 1241 if (!Tmp) { 1242 PyErr_SetString(PyExc_RuntimeError, 1243 "quantity_ext.c: limit_by_vertex could not obtain beta_w object from domain"); 1244 1245 return NULL; 1246 } 1247 1248 1249 // Get pertinent variables 1250 neighbours = get_consecutive_array(domain, "neighbours"); 1251 centroid_values = get_consecutive_array(quantity, "centroid_values"); 1252 vertex_values = get_consecutive_array(quantity, "vertex_values"); 1253 edge_values = get_consecutive_array(quantity, "edge_values"); 1254 beta_w = PyFloat_AsDouble(Tmp); 1255 1256 1257 N = centroid_values > dimensions[0]; 1258 1259 err = _limit_edges_by_neighbour(N, beta_w, 1260 (double*) centroid_values > data, 1261 (double*) vertex_values > data, 1262 (double*) edge_values > data, 1263 (long*) neighbours > data); 1264 1265 if (err != 0) { 1266 PyErr_SetString(PyExc_RuntimeError, 1267 "Internal function _limit_by_vertex failed"); 1268 return NULL; 1269 } 1270 1271 1272 // Release 1273 Py_DECREF(neighbours); 1274 Py_DECREF(centroid_values); 1275 Py_DECREF(vertex_values); 1276 Py_DECREF(edge_values); 1277 Py_DECREF(Tmp); 1278 1279 1280 return Py_BuildValue(""); 1281 } 1282 1145 1283 1146 1284 // Method table for python module 1147 1285 static struct PyMethodDef MethodTable[] = { 1148 1286 {"limit_old", limit_old, METH_VARARGS, "Print out"}, 1149 {"limit_by_vertex", limit_by_vertex, METH_VARARGS, "Print out"}, 1150 {"limit_by_edge", limit_by_edge, METH_VARARGS, "Print out"}, 1287 {"limit_vertices_by_all_neighbours", limit_vertices_by_all_neighbours, METH_VARARGS, "Print out"}, 1288 {"limit_edges_by_all_neighbours", limit_edges_by_all_neighbours, METH_VARARGS, "Print out"}, 1289 {"limit_edges_by_neighbour", limit_edges_by_neighbour, METH_VARARGS, "Print out"}, 1151 1290 {"update", update, METH_VARARGS, "Print out"}, 1152 1291 {"backup_centroid_values", backup_centroid_values, METH_VARARGS, "Print out"},
Note: See TracChangeset
for help on using the changeset viewer.